我试图实现最大似然估计以下2参数贝塔泊松模型
通过研究StackOverflow中的其他解决方案，我得出了一个结论（可能是错误的），即“optim”函数可能是最好的解决方案。因此，我尝试编写了以下代码：
Log_Lik_BP <- function(alpha_beta, D, Y, N){
alpha <- alpha_beta[1]
beta <- alpha_beta[2]
sum(Y * log(1 - (1 + (D/beta))^(-alpha))) - alpha*sum(N-Y)*log(1+(D/beta))
}
optim(par = c(0,0), fn = Log_Lik_BP, D = df$D, N = df$N, Y = df$Y)但是我没有非常成功，我收到了以下错误消息，其中“长度22”指的是 Dataframe 中的22个数据点
Error in optim(par = c(0, 0), fn = Log_Lik_BP, D = df$D, N = df$N, Y = df$Y) :
objective function in optim evaluates to length 22 not 1我应用此MLE估计的数据集如下：
D <- c(4.2,8,10.6,11.9,12.2,12.6,13,13.2,28.3,28.7,68.2,69,71.1,83,91.7,95.5,106.5,136.5,171.2,186.9,309.3,1557.1)
N <- c(1,1,1,1,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)
Y <- c(0,0,1,1,2,1,0,1,0,1,1,0,1,1,1,1,1,1,1,1,1,1)
df <- data.frame(D, N, Y)提前感谢您对我的代码的建议和更正
1条答案
按热度按时间gpnt7bae1#
这里有一堆问题，我将介绍调试过程。
在解决优化问题时，首先在起始参数（例如
Log_Lik_BP(c(0,0), D, N, Y)）处评估函数以消除一层复杂性总是很有帮助的。此外，它有助于（1）将debug()或debugonce()应用于函数以逐步完成它，或者（2）在函数范围之外运行代码以查看发生了什么。1.运行
Log_Lik_BP(c(0,0), D, N, Y))，结果：一个22个NaN值的向量。我们需要做的第一件事是看看为什么它的长度是22而不是1（我们稍后会回到NaN问题）。1.设置
alpha <- beta <- 0并评估单个组件sum(Y * log(1 - (1 + (D/beta))^(-alpha)))和alpha*sum(N-Y)*log(1+(D/beta))。结果：第一个组成部分是长度-1，正如它应该的那样，第二个组成部分是长度22。啊哈！公式中的分组不正确，我们需要括号来包括log(1+(D/beta))项，即alpha*sum((N-Y)*log(1+(D/beta)))。1.相应地固定
Log_Lik_BP。结果：如果我们在c(0,0)处对函数求值，我们再次得到NaN;如果我们尝试optim()，我们得到“函数不能在初始参数处求值”。1.从公式来看，
c(0,0)似乎是一个有问题的起点，因为我们需要除以beta。如果我们从c(1,1)开始，会发生什么？结果：optim()运行时没有错误，但我们得到了看起来很有趣的答案-默认的Nelder-Mead算法用完了迭代，看起来beta正在收敛到零。1.有点跑题了，但我尝试了另一种方法，允许设置边界（
method = "L-BFGS-B"）并设置alpha和beta（lower = c(0.001, 0.001)）的下限（当beta=0停止时，将下限设置为0会遇到麻烦，因为“L-BFGS-B需要'fn'的有限值”）。结果：这很快就得到了一个解，但是结果看起来仍然不稳定（alpha很大，而beta达到了边界）。（这一步可能是不必要的，请参阅下一步以获得剩余问题的真实的解决方案。）1.要知道
optim()默认情况下会最小化目标函数，并且您已经定义了一个应该被最大化的对数似然函数。将control = list(fnscale = -1)改为最大化。结果：合理的参数估计和负对数似然。您仍然应该检查结果在上下文中是否有意义！