NaN ошибки с bbmle

Question

Этот вопрос относится к моему предыдущему вопросу здесь и к набору данных, представленному в статье Новое обобщение линейного экспоненциального распределения: теория и применение .Для этих данных, адаптируя код, предложенный Беном Болкером, у нас есть

library(stats4)
library(bbmle)

x <- scan(textConnection("115 181 255 418 441 461 516 739 743 789
807 865 924 983 1024 1062 1063 1165 1191 1222 1222 1251 1277 1290 1357 1369 1408 1455 1478 1549 1578 1578 1599 1603 1605 1696 1735 1799 1815 1852"))

dd  <- data.frame(x)

dLE <- function(x,lambda,theta,log=TRUE){
    r <- log(lambda+theta*x)-(lambda*x+(theta/2)*x^2)
    if (log) return(r) else return(exp(r))
}

svec <- list(lambda=0.0009499,theta=0.000002)
m1 <- mle2( x ~ dLE(lambda,theta),
      data=dd,
      start=svec,
      control=list(parscale=unlist(svec)))
coef(m1)

, который возвращает несколько ошибок (произведенных NaN) и значения для mles, которые весьма отличаются от приведенных в таблице 2 статьи.Почему это так и как это можно исправить?

Answer 1

После некоторого исследования мое мнение таково, что статья просто дает неверные результаты. Результаты, полученные от optim(), дают результаты, которые выглядят намного лучше, чем те, о которых сообщается в статье. Я всегда мог что-то упустить; Я бы предложил вам связаться с соответствующим автором.

(Предупреждения не обязательно являются проблемой - они означают, что оптимизатор пробовал некоторые комбинации, которые приводят к получению журналов отрицательных чисел по пути, что не означает, что конечный результат неверен - но я согласен, что это всегда хорошая идея для устранения предупреждений в случае, если они каким-то образом портят результат.)

отборочные

library(bbmle)
## load data, in a format as similar to original table
## as possible (looking for typos)
x <- scan(textConnection("115 181 255 418 441 461 516 739 743 789
                          807 865 924 983 1024 1062 1063 1165 1191 1222
                         1222 1251 1277 1290 1357 1369 1408 1455 1478 1549
                         1578 1578 1599 1603 1605 1696 1735 1799 1815 1852"))
dd  <- data.frame(x)  
## parameters listed in table 2
svec <- list(lambda=9.499e-4,theta=2e-6)

Функции

## PDF (as above)
dLE <- function(x,lambda,theta,log=TRUE){
    r <- log(lambda+theta*x)-(lambda*x+(theta/2)*x^2)
    if (log) return(r) else return(exp(r))
}
## CDF (for checking)    
pLE <- function(x,lambda,theta) {
    1-exp(-(lambda*x+(theta/2)*x^2))
}

подходит модель

Я использовал method="L-BFGS-B", потому что это облегчает установку нижних границ параметров (что позволяет избежать предупреждений).

m1 <- mle2( x ~ dLE(lambda,theta),
      data=dd,
      start=svec,
      control=list(parscale=unlist(svec)),
      method="L-BFGS-B",
      lower=c(0,0))

Результаты

coef(m1)
##      lambda        theta 
## 0.000000e+00 1.316733e-06 
-logLik(m1)  ## 305.99 (much better than 335, reported in the paper)

* * 1 022 графа

Давайте еще раз проверим, можем ли мы воспроизвести этот рисунок из бумаги:

png("SO55032275.png")
par(las=1)
plot(ecdf(dd$x),col="red")
with(svec,curve(pLE(x,lambda,theta),add=TRUE,col=1))
with(as.list(coef(m1)),curve(pLE(x,lambda,theta),add=TRUE,col=3,lty=2))
legend("topleft",col=c(2,1,3),lty=c(NA,1,3),pch=c(16,NA,NA),
       c("ecdf","paper (lam=9e-4, th=2e-6)","ours (lam=0, th=1.3e-6)"))
dev.off()

Ecdf и CDF, нарисованные с параметрами из бумажного соответствия; CDF, нарисованный с параметрами, оцененными здесь, намного лучше (фактически он выглядит лучше и имеет меньшую логарифмическую вероятность, чем соответствие KLE, о котором сообщалось в статье). Я пришел к выводу, что что-то не так с припадками на бумаге.