使用R中optim()估计gamma_estimate的MLE遇问题求解决

我用R的optim()函数估计gamma_estimate参数的MLE，设置参数下界为1、上界为1.9，但遇到以下问题：

• 使用"Brent"和"L-BFGS-B"方法时，无论初始值和边界怎么设，结果始终是下界值；
• 无边界的"Nelder-Mead"方法得到0.07711332，显然不符合预设的边界要求，结果不可靠；
• 使用"SANN"方法时直接报错。

附上原自定义目标函数和optim()调用代码：

原自定义函数代码 #

gamma_mle <- function(gamma_estimates) {

  # Initialization

  Mu <-test[nrow(test), 2:4] %>% as.double()
  var <- test[nrow(test), 5:7] %>% as.double()
  pi <- test[nrow(test), 8:10] %>% as.double()
  alpha <- test[nrow(test), 11:13] %>% as.double()
  beta <- test[nrow(test), 14:16] %>% as.double()
  loglike <- c(0, test[nrow(test), 17]) %>% as.double()
  read_cover <- 100
  n = 2
  k = 3

  
  ## E step
  tau <- matrix("0", k) %>% as.list()
  total <- matrix(0, nrow = nrow(data))
  
  for (i in 1:k) {
    tau[[i]]<- pi[i] * dbeta(data, alpha[i], beta[i])}
  
  for (i in 1:k) {
    total <- total + tau[[i]]}
  
  for (i in 1:k) {
    tau[[i]]<- tau[[i]] %>% as.matrix() / total %>% as.matrix()
  }
  
  # Update Rule
  # Update mixing coefficients
  for (i in 1:k) {
    pi[[i]] <- sum(tau[[i]]) / length(data)
  }
  
  # Update Mu
  for (i in 1:k) {
    Mu[[i]] <- sum(tau[[i]] * data)/sum(tau[[i]])
  }
  
  # Update variance
  for (i in 1:k) {
    var[[i]] <- ((Mu[[i]] * (1 - Mu[[i]])) * gamma_estimates) / read_cover
  }
  
  # Convert Mu, var to updated alpha, beta 
  updated <- matrix("0", k) %>% as.list()
  for (i in 1:k) {
    updated[[i]] <- estBetaParams(Mu[[i]], var[[i]]) %>% unlist()
    alpha[[i]] <- c(updated[[i]][1])
    beta[[i]] <- c(updated[[i]][2])
  }
  
  ## M step
  #Maximize the loglikelihood
  for (i in 1:k) {
    loglike <-sum(log(pi[i] * dbeta(data, alpha[i], beta[i])))
  }
  return(loglike)
}

原optim()调用代码 #

gamma_estimates <- optim(fn = gamma_mle, 
                         par = c(1.50), # parameter = k, gamma_estimates
                         lower = c(1.000), #Lower bound on parameters
                         upper = c(1.990), #Upper bound on parameters
                         hessian = FALSE, #Return Hessian for SEs
                         method = "Brent", 
                         control = list(maxit = 10000, parscale = 0.00001)
)

gamma_estimates <- optim(fn = gamma_mle,
                         par = c(1.30), # parameter = k, gamma_estimates
                         lower = c(1.00), #Lower bound on parameters
                         upper = c(1.6), #Upper bound on parameters
                         hessian = FALSE, #Return Hessian for SEs
                         method = "L-BFGS-B"
)

gamma_estimates <- optim(fn = gamma_mle,
                         par = c(1.5), # parameter = k, gamma_estimates
                         hessian = FALSE, #Return Hessian for SEs
                         method = "Nelder-Mead",
                         control = list(maxit = 10000, parscale = 0.00001)
)

gamma_estimates <- optim(fn = gamma_mle, 
                         par = c(1.01), # parameter = k, gamma_estimates
                         #lower = c(-Inf), #Lower bound on parameters
                         #upper = c(Inf), #Upper bound on parameters
                         hessian = TRUE, #Return Hessian for SEs
                         method = "SANN"
)

核心问题排查 #

1. 目标函数方向错误：optim()默认是最小化目标函数，但你需要最大化对数似然。当前函数直接返回对数似然，优化算法会寻找最小值，这也是边界方法始终返回下界的核心原因（如果对数似然随gamma增大而递减，最小值会落在左边界）。
2. 对数似然计算错误：M步骤的循环中，每次迭代都会覆盖loglike变量，最终只返回第k个混合分量的对数似然，而非整个混合模型的对数似然总和。
3. 参数更新逻辑失效：函数中Mu、pi、alpha等参数每次调用都从test最后一行读取初始值，没有在E-M步骤中迭代更新，相当于每次优化迭代都重置了模型参数，完全违背E-M算法的迭代逻辑。
4. 潜在辅助函数问题：estBetaParams函数的实现是否正确，直接影响alpha和beta的更新，进而干扰对数似然计算。

修正后的代码实现 #

1. 修正目标函数

# 先确保estBetaParams函数正确实现（Beta分布参数转换）
estBetaParams <- function(mu, var) {
  if (mu <= 0 || mu >= 1 || var <= 0) stop("无效的均值或方差，需满足0<mu<1且var>0")
  phi <- (mu * (1 - mu)) / var - 1
  alpha <- mu * phi
  beta <- (1 - mu) * phi
  return(c(alpha = alpha, beta = beta))
}

gamma_mle <- function(gamma_est) {
  # 读取初始模型参数
  init_params <- test[nrow(test), ]
  Mu <- as.double(init_params[2:4])
  pi <- as.double(init_params[8:10])
  alpha <- as.double(init_params[11:13])
  beta <- as.double(init_params[14:16])
  read_cover <- 100
  k <- 3
  
  ## E step：计算后验概率tau
  tau <- vector("list", k)
  total <- numeric(nrow(data))
  for (i in 1:k) {
    tau[[i]] <- pi[i] * dbeta(data, alpha[i], beta[i])
  }
  total <- rowSums(sapply(tau, cbind))
  tau <- lapply(tau, function(x) x / total)
  
  ## M step：更新模型参数
  # 更新混合系数pi
  pi_new <- sapply(tau, sum) / length(data)
  # 更新Mu
  Mu_new <- sapply(1:k, function(i) sum(tau[[i]] * data) / sum(tau[[i]]))
  # 更新方差（依赖gamma_est）
  var_new <- (Mu_new * (1 - Mu_new) * gamma_est) / read_cover
  # 更新alpha和beta
  ab_new <- lapply(1:k, function(i) estBetaParams(Mu_new[i], var_new[i]))
  alpha_new <- sapply(ab_new, `[`, 1)
  beta_new <- sapply(ab_new, `[`, 2)
  
  # 计算整个混合模型的对数似然
  loglike <- sum(log(rowSums(sapply(1:k, function(i) {
    pi_new[i] * dbeta(data, alpha_new[i], beta_new[i])
  }))))
  
  # 返回负对数似然，适配optim的最小化逻辑
  return(-loglike)
}

2. 修正optim调用

# 使用L-BFGS-B方法，最小化负对数似然等价于最大化对数似然
gamma_estimates <- optim(
  fn = gamma_mle,
  par = 1.5,  # 初始值
  lower = 1.0,
  upper = 1.9,
  method = "L-BFGS-B",
  control = list(maxit = 10000)  # 可根据收敛情况调整maxit
)

# 查看结果
print(gamma_estimates)

额外注意事项 #

• 参数传入优化：建议将初始模型参数（Mu、pi等）作为gamma_mle的参数传入，避免依赖全局变量，降低代码耦合风险。
• 数据有效性检查：确保data中的值在(0,1)范围内（符合Beta分布取值），否则dbeta会返回0或NaN，导致对数似然计算错误。
• 收敛性验证：查看gamma_estimates$convergence值，0表示收敛成功；若未收敛，可调整control中的parscale（参数缩放）或maxit（最大迭代次数），或尝试多组初始值。
• SANN方法适配：模拟退火方法对目标函数平滑性要求高，需额外调整温度参数，不建议优先使用；若要尝试，需确保目标函数无NaN值且连续。

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