在 bbmle 包中实现 MLE 估计的问题(对于 R)
Problem Implementing MLE Estimation in bbmle Package (for R)
我正在尝试验证为 $\alpha$、$\beta$ 和 $ 获得的 MLE \lambda$ 用于 Zubair 等人在使用数据集 1 时标题为 A Study of Logistic-Lomax Distribution 的论文中的 Logistic-Lomax 分布。该论文使用以下代码来执行此操作(请参阅附录 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{\beta} = 17.6331, \hat{\lambda} = 35.6364$使用此代码的此数据集的 MLE。我的问题很简单:如何在 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)
这里发生了什么,我该如何解决?
有两个问题。
第一个
> 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}{\beta}$。代码有
n*log(lambda*alpha/beta)
。这些看起来一样,但是 mle2
是一个 minimiser,所以它们应该有相反的符号。对数似然的方程与方程 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 中的代码,但即使获得代码也是此类论文的一大进步
我正在尝试验证为 $\alpha$、$\beta$ 和 $ 获得的 MLE \lambda$ 用于 Zubair 等人在使用数据集 1 时标题为 A Study of Logistic-Lomax Distribution 的论文中的 Logistic-Lomax 分布。该论文使用以下代码来执行此操作(请参阅附录 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{\beta} = 17.6331, \hat{\lambda} = 35.6364$使用此代码的此数据集的 MLE。我的问题很简单:如何在 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)
这里发生了什么,我该如何解决?
有两个问题。
第一个
> 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}{\beta}$。代码有
n*log(lambda*alpha/beta)
。这些看起来一样,但是 mle2
是一个 minimiser,所以它们应该有相反的符号。对数似然的方程与方程 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 中的代码,但即使获得代码也是此类论文的一大进步