Ошибки NaN с bbmle

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

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 статьи. Почему это так и как это можно исправить?

r nan mle

Will 06.03.2019 источник

Ответы (1)

arrow_upward
2
arrow_downward

После некоторых исследований я пришел к выводу, что статья просто имеет неверные результаты. Результаты, которые я получаю от 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)

график

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

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-подгонка, указанная в статье). Я прихожу к выводу, что с вписыванием в бумагу что-то не так.

Ben Bolker 07.03.2019

comment

Большое спасибо. Да, я был озадачен их результатами, особенно отрицательным результатом, который я получил. Я еще не проверял параметры других подогнанных моделей. - Will; 08.03.2019