在bbmle包中实现MLE估计的问题(针对R)

Question

我正在试图验证Zubair等人在使用数据集1时，在题为$\alpha$ -Lomax分布的文章中为$\beta$、$\beta$和$\lambda$获得的MLEs。本文使用了以下代码(见附录B)：

library(bbmle)
x <- c(66, 117, 132, 111, 107, 85, 89, 79, 91, 97, 138, 103, 111, 86, 78, 96, 93, 101, 102, 110, 95, 96, 88, 122, 115, 92, 137, 91, 84, 96, 97, 100, 105, 104, 137, 80, 104, 104, 106, 84, 92, 86, 104, 132, 94, 99, 102, 101, 104, 107, 99, 85, 95, 89, 102, 100, 98, 97, 104, 114, 111, 98, 99, 102, 91, 95, 111, 104, 97, 98, 102, 109, 88, 91, 103, 94, 105, 103, 96, 100, 101, 98, 97, 97, 101, 102, 98, 94, 100, 98, 99, 92, 102, 87 , 99, 62, 92, 100, 96, 98) 
n <- length(x)
ll_LLx <- function(alpha, beta, lambda){
n*log(lambda*alpha/beta)-sum(log(1+x/beta))
-(lambda+1)*sum(log(log((1+x/beta)^alpha)))
-2*sum(log(1+(log((1+x/beta)^alpha))^(-lambda)))
}
mle.res <- mle2(ll_LLx, start=list(alpha=alpha, beta=beta, lambda=lambda),
hessian.opt=TRUE)
summary(mle.res)

本文使用此代码得到了该数据集的$\hat{\alpha} = 0.5302，\hat{\β}= 17.6331，\hat{λ}= 35.6364$的值。我的问题很简单:如何在R中实现这段代码而不产生错误？这段代码似乎列出了参数'alpha'，'beta‘和'lambda'，但没有给它们分配数值起始值。因此，我试图在代码之前为这些参数设置合理的起始值，如下所示：

alpha=0.5
beta=17
lambda=35

然而，这产生了意想不到的错误：

Error in optim(par = c(alpha = 0.5, beta = 17, lambda = 35), fn = function (p)  : 
  non-finite finite-difference value [1]
In addition: There were 50 or more warnings (use warnings() to see the first 50)

这里发生了什么，我该怎么解决呢？

Answer 1

有两个问题。

第一

> alpha=0.53
> beta=17.6
> lambda=35.6
> ll_fragment<-function(alpha,beta,gamma) -2*sum(log(1+(log((1+x/beta)^alpha))^(-lambda)))
> 
> ll_LLx(alpha,beta,gamma)
[1] -205.132
> ll_fragment(alpha,beta,gamma)
[1] -205.132

也就是说，在打印引入的换行代码时，当您从PDF中复制代码时，您将得到一系列的三个表达式。函数返回的值就是最后一个表达式的值。

第二，如果将代码与等式中定义的逻辑似然进行比较.在方程6.1之前的未编号方程中定义的方程以$n\log\frac{\lambda\alpha}{\β}$开始。代码有n*log(lambda*alpha/beta)。这些看起来一样，但mle2是一个最小的，所以他们应该有相反的标志。逻辑似然的方程与方程2.2中的pdf相匹配，所以我假设它是正确的，并且给出的代码将试图最小化对数可能性。

如果我修理断线和标牌，我就会得到

> mle.res <- mle2(ll1a, start=list(alpha=alpha, beta=beta, lambda=lambda),hessian=TRUE)
Warning messages:
1: In log(lambda * alpha/beta) : NaNs produced
2: In log(log((1 + x/beta)^alpha)) : NaNs produced
3: In log(lambda * alpha/beta) : NaNs produced
4: In log(log((1 + x/beta)^alpha)) : NaNs produced
5: In log(lambda * alpha/beta) : NaNs produced
6: In log(log((1 + x/beta)^alpha)) : NaNs produced
> summary(mle.res)
Maximum likelihood estimation

Call:
mle2(minuslogl = ll1a, start = list(alpha = alpha, beta = beta, 
    lambda = lambda), hessian.opts = TRUE)

Coefficients:
        Estimate Std. Error z value     Pr(z)    
alpha   0.530499   0.018522  28.641 < 2.2e-16 ***
beta   17.649226   1.357944  12.997 < 2.2e-16 ***
lambda 35.636033   2.260395  15.765 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

-2 log L: 771.2541 

与报纸完全一致。

这就是为什么应该需要存储库中的代码，而不是PDF中的代码，但对于这类论文来说，即使是获取代码也是向前迈进了一大步。