R中优化伪收敛：求解MLE时积分处理方案咨询

我明白你在拟合带有复杂积分项的负对数似然时遇到了伪收敛，核心问题是当前的Gauss-Laguerre积分精度不足，导致似然函数的梯度/曲率计算不稳定，让优化器误判为收敛。下面是几个针对性的改进方案，你可以逐一尝试：

1. 优化Gauss-Laguerre积分的参数与实现 #

你当前的积分实现没有充分匹配被积函数的形式，通过变量替换和调整积分参数可以大幅提升精度：

变量替换适配积分形式 #

观察你的被积函数：
$$\frac{y^{{a-\mu-1}}{(y+t)}{x+a}} e^{-\theta/y}$$
做变量替换 $z = \theta/y$（即 $y = \theta/z$），积分可转化为更适配Gauss-Laguerre权重（$e^{-z}$）的形式：
$$\theta^{a-\mu} \int_0^\infty \frac{z^{x+\mu - a}}{(\theta + t z)^{x+a}} e^{-z} dz$$

改用statmod包的高效积分实现 #

statmod的Gauss-Laguerre实现更稳定，结合变量替换后的代码如下：

library(statmod)

loglik <- function(para, data) {
  n.quad <- 128  # 适当增加节点数提升精度
  a <- para[1]
  mu <- para[2]
  the <- para[3]
  
  # 生成适配替换后形式的Gauss-Laguerre节点与权重
  quad <- gauss.quad(n.quad, kind = "laguerre", alpha = mu - a)
  nodes <- quad$nodes
  weights <- quad$weights
  
  int.ing <- function(z, x, t) {
    # 变量替换后的被积函数
    (the^(a - mu) * z^(x + mu - a)) / (the + t * z)^(x + a)
  }
  
  int.f <- function(x, t) {
    sum(weights * int.ing(nodes, x, t))
  }
  
  v.int.f <- Vectorize(int.f)
  int <- v.int.f(data$count, data$time)
  
  loglik.value <- lgamma(a + data$count) - lgamma(a) + mu*log(the) - lgamma(mu) + log(int)
  log.sum <- sum(loglik.value)
  return(-log.sum)
}

2. 用自适应积分替代固定节点积分 #

R的cubature包提供自适应单/多变量积分，能自动调整精度直到满足阈值，避免固定节点数导致的精度不足：

library(cubature)

loglik <- function(para, data) {
  a <- para[1]
  mu <- para[2]
  the <- para[3]
  
  int.f <- function(x, t) {
    integrand <- function(y) {
      (y^(a - mu - 1) / (y + t)^(x + a)) * exp(-the / y)
    }
    # 设置严格精度阈值，平衡精度与速度
    adaptIntegrate(integrand, lowerLimit = 0, upperLimit = Inf, tol = 1e-8)$integral
  }
  
  v.int.f <- Vectorize(int.f)
  int <- v.int.f(data$count, data$time)
  
  loglik.value <- lgamma(a + data$count) - lgamma(a) + mu*log(the) - lgamma(mu) + log(int)
  log.sum <- sum(loglik.value)
  return(-log.sum)
}

3. 优化优化器设置与初始值 #

伪收敛也可能源于初始值偏离最优值或优化器约束过松：

• 调整初始值：模拟数据的真实theta为480，将初始值设为c(2, 0.01, 400)更接近真实值；
• 放宽优化器限制：增加迭代次数、缩小精度阈值，同时给参数设置极小的下界（避免趋近于0时的数值问题）：

start.value <- c(2, 0.01, 400)
fit <- nlminb(start = start.value, loglik, data=data, lower = c(1e-6, 1e-6, 1e-6), 
              control = list(trace = T, iter.max = 200, rel.tol = 1e-10))

4. 提升积分的数值稳定性 #

当参数a、mu接近时，被积函数可能出现数值溢出现象，可通过对数转换避免：

int.ing.log <- function(y, x, t, a, mu, the) {
  (a - mu - 1)*log(y) - (x + a)*log(y + t) - the/y
}

int.f <- function(x, t, a, mu, the) {
  integrand <- function(y) {
    exp(int.ing.log(y, x, t, a, mu, the))
  }
  adaptIntegrate(integrand, 0, Inf)$integral
}

尝试以上方法后，应该能解决积分收敛慢的问题，让优化器找到真正的收敛点。

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