Оценка максимального правдоподобия с равномерным распределением по R приводит к абсурдному результату

Question

Я хочу использовать функцию mle, чтобы получить оценки a и b в распределении Unif(a,b). Но я получаю абсурдные оценки нигде близко к 1 и 3.

library(stats4) 
set.seed(20161208) 

N <- 100 
c <- runif(N, 1, 3) 
LL <- function(min, max) { 
    R <- runif(100, min, max) 
    suppressWarnings((-sum(log(R)))) 
    } 
mle(minuslogl = LL, start = list(min = 1, max = 3), method = "BFGS", 
    lower = c(-Inf, 0), upper = c(Inf, Inf)) 

я получил:

Call: 
mle(minuslogl = LL, start = list(min = 1, max = 3), method = "BFGS") 

Coefficients: 
    min  max 
150.8114 503.6586 

Любые идеи, что происходит? Заранее спасибо!

Answer 1

Я бы хотел указать, где ваш код является неправильным.

1. Вам нужно dunif не runif. Вы можете определить:

   LL <- function (a, b) -sum(dunif(x, a, b, log.p = TRUE)) 
   

   В моем коде ниже я не использовал dunif, так как плотность только 1/(b - a), так что я написал это непосредственно.

2. Вы генерируете образцы внутри целевой функции. Для U[a,b] это нормально, так как его плотность не содержит x. Но для других распределений объектная функция изменяется на каждой итерации.
3. С короткими ограничениями вам нужно method = "L-BFGS-B", а не обычное "BFGS". И вы не используете правильные ограничения.

---

Теперь более подробно ...

Для длина- вектора n образца x из U[a, b], вероятности является (b - a)^(-n) и отрицательные лог-правдоподобия n * log(b - a). Очевидно, что MLE - a = min(x) и b = max(x).

Численная оптимизация совершенно не нужна и на самом деле невозможна без ограничений. Посмотрите на вектор градиента:

(n/(a - b), n/(b - a)) 

Частная производная w.r.t. a/b всегда негатив/позитив и не может быть 0.

Численный подход становится возможным при наложении ограничений коробки: -Inf < a <= min(x) и max(x) <= b < Inf. Мы точно знаем, что итерация заканчивается на границе.

Мой код ниже использует как optim, так и mle. Примечание mle потерпит неудачу, когда он переворачивает матрицу Гессе, как в единственном числе:

-(b - a)^2 (b - a)^2 
(b - a)^2 -(b - a)^2 

Код:

## 100 samples 
set.seed(20161208); x <- runif(100, 1, 3) 
# range(x) 
# [1] 1.026776 2.984544 

## using `optim` 
nll <- function (par) log(par[2] - par[1]) ## objective function 
gr_nll <- function (par) c(-1, 1)/diff(par) ## gradient function 
optim(par = c(0,4), fn = nll, gr = gr_nll, method = "L-BFGS-B", 
     lower = c(-Inf, max(x)), upper = c(min(x), Inf), hessian = TRUE) 
#$par 
#[1] 1.026776 2.984544 ## <- reaches boundary! 
# 
# ... 
# 
#$hessian ## <- indeed singular!! 
#   [,1]  [,2] 
#[1,] -0.2609022 0.2609022 
#[2,] 0.2609022 -0.2609022 

## using `stats4::mle` 
library(stats4) 
nll. <- function (a, b) log(b - a) 
mle(minuslogl = nll., start = list(a = 0, b = 4), method = "L-BFGS-B", 
    lower = c(-Inf, max(x)), upper = c(min(x), Inf)) 
#Error in solve.default(oout$hessian) : 
# Lapack routine dgesv: system is exactly singular: U[2,2] = 0