使用optim函数实现双参数Beta-Poisson模型MLE的参数问题排查

问题描述 #

我尝试为双参数Beta-Poisson模型实现极大似然估计（MLE），参考资料后使用optim函数编写代码，但运行时出现错误：目标函数返回长度为22而非1，恳请提供修正建议。

原代码 #

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)

错误信息 #

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

数据集 #

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)

修正建议 #

错误原因 #

optim要求目标函数返回单个标量值（对数似然的总和），原代码中第二部分- alpha*sum(N-Y)*log(1+(D/beta))未对(N-Y)*log(1+(D/beta))整体求和，导致返回长度为22的向量，不符合要求。此外初始值c(0,0)会导致分母为0的计算错误（beta=0时D/beta无意义），且optim默认最小化目标函数，而我们需要最大化对数似然。

修正后的代码 #

方式1：直接最大化对数似然

Log_Lik_BP <- function(alpha_beta, D, Y, N){
  alpha <- alpha_beta[1]
  beta <- alpha_beta[2]
  # 限制参数为正，避免无效计算
  if(alpha <= 0 || beta <= 0) return(-Inf)
  # 修正求和逻辑，确保返回单个标量
  sum(Y * log(1 - (1 + (D/beta))^(-alpha))) - alpha*sum( (N-Y)*log(1+(D/beta)) )
}

# 使用合理初始值，通过fnscale=-1实现最大化
optim(par = c(1, 10), fn = Log_Lik_BP, D = df$D, N = df$N, Y = df$Y, 
      control = list(fnscale = -1))

方式2：最小化负对数似然（更常用）

Neg_Log_Lik_BP <- function(alpha_beta, D, Y, N){
  alpha <- alpha_beta[1]
  beta <- alpha_beta[2]
  if(alpha <= 0 || beta <= 0) return(Inf)
  -sum(Y * log(1 - (1 + (D/beta))^(-alpha))) + alpha*sum( (N-Y)*log(1+(D/beta)) )
}

optim(par = c(1, 10), fn = Neg_Log_Lik_BP, D = df$D, N = df$N, Y = df$Y)

关键修正点 #

• 求和逻辑修正：对(N-Y)*log(1+(D/beta))整体求和，确保目标函数返回单个数值。
• 参数范围限制：添加alpha和beta的正性判断，避免计算中出现NaN或无穷大，提升优化稳定性。
• 合理初始值：替换c(0,0)为c(1,10)这类有效初始参数，避免分母为0的错误。
• 优化方向调整：通过fnscale=-1或负对数似然，适配optim的默认最小化行为，实现对数似然最大化。

内容的提问来源于stack exchange，提问作者user2085797