与 bbmle 的 NaN 错误

Question

这个问题与我之前的问题和论文A New Generalization of Linear Exponential Distribution中提出的数据集有关：理论与应用.对于此数据，采用 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)

其中 returns 几个错误（产生 NaN）和 mles 的值与论文 Table 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 拟合更好，并且具有更低的对数似然）。我的结论是论文中的拟合存在严重错误。

与 bbmle 的 NaN 错误

NaN errors with bbmle

r

nan

mle

预赛

函数

适合模特

结果

图表