与bbmle的NaN错误

Question

这个问题与我之前的问题here和A New Generalization of Linear Exponential Distribution: Theory and Application论文中提供的数据集有关。对于这些数据，我们有适应Ben Bolker提出的代码

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)

图形

让我们仔细看看是否可以从论文中复制这个数字：

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拟合更低的对数似然性）。我得出结论，论文的适用性存在严重问题。