R语言Newton Raphson法做MLE报初始值越界、生成NaN问题求解

• 似然函数推导错误：你当前编写的对数似然完全不符合删失数据的极大似然估计结构，且存在大量数值溢出、负数取对数的问题：
  1. 冗余嵌套的exp()/log()计算极易触发数值上溢，比如exp(t^beta)会快速变为无穷大，导致后续计算出现负数，取对数后直接生成NaN
  2. (data$t_i)^(beta - 1)在beta<1时，若t_i取值极小会直接变成无穷大，和后续项相乘后要么溢出要么变负，触发log计算报错
• 未加参数约束：lambda和beta均为必须大于0的参数，未给优化器加约束的情况下，迭代过程很容易跑出负参数值，代入计算直接生成NaN
• 无数值保护逻辑：所有涉及exp、log的计算没有做截断处理，极易出现上溢/下溢问题

1. 重新推导正确的对数似然 #

根据你的数据生成逻辑，对应生存函数为S(t) = exp(-lambda * (exp(t^beta) - 1))，删失数据的对数似然规则为：事件发生的观测加log(f(t))，删失观测加log(S(t))，原代码中的似然结构完全错误，这是核心问题。

2. 增加数值稳定性保护 #

对所有可能出现0、负数、溢出的计算节点做截断处理，避免非法值输入到log/exp函数中。

3. 增加参数约束 #

强制优化过程中lambda和beta始终为正，避免参数越界。

修复后可运行代码 #

lambda <- 0.02
beta <- 0.5
n <- 100
N <- 1000
lambda_hat <- beta_hat <- cp <- NULL

library(survival); library(maxLik)

set.seed(20)
for (i in 1:N) {
  u <- runif(n)
  c_i <- rexp(n, 0.00001)
  # 生成数据时做截断，避免t_i异常过大
  t_i <- pmin((log(pmax(1 - (1/lambda)*log(1 - u), 1e-10)))^(1/beta), 1e6)
  s_i <- 1*(t_i < c_i)
  t <- pmin(t_i, c_i)
  data <- data.frame(t = t, s_i = s_i)
  
  LLF <- function(para) {
    lambda <- para[1]
    beta <- para[2]
    # 参数越界直接返回负无穷，优化器会自动放弃该迭代方向
    if(lambda <= 0 || beta <= 0) return(-Inf)
    # 截断t^beta避免exp上溢，exp(20)已足够大无需计算更大值
    t_beta <- pmin(t^beta, 20)
    exp_tbeta <- exp(t_beta)
    # 事件观测的对数似然
    ll_event <- s_i * (log(lambda) + log(beta) + (beta - 1)*log(pmax(t, 1e-10)) + t_beta)
    # 删失观测的对数似然
    ll_cens <- (1 - s_i) * (-lambda*(exp_tbeta - 1))
    return(sum(ll_event + ll_cens))
  }
  # 调用maxLik，加参数正约束，换用更稳健的BFGS优化方法
  mle <- maxLik(LLF, start = c(lambda = 0.02, beta = 0.5), 
                constraints = list(ineqA = matrix(c(1,0, 0,1), nrow=2, byrow=T), 
                                   ineqB = c(1e-5, 1e-5)),
                method = "BFGS")
  lambda_hat[i] <- mle$estimate[1]
  beta_hat[i] <- mle$estimate[2]
}

如果需要换回Newton-Raphson方法，可以先用上述BFGS方法跑出收敛的参数值作为初始值，再代入Newton-Raphson迭代即可。

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