r - bbmle 的 NaN 错误

Question

这个问题与我之前的问题和论文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 中给出的完全不同。为什么会这样，如何纠正？

标签： rnanmle

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 拟合要低）。我的结论是论文中的拟合存在严重错误。