R语言SFA模型最大似然估计函数优化及参数校准问题问询

一、适配的优化方法 #

• 更换无梯度优化算法：BFGS依赖梯度计算，SFA似然函数包含pnorm项，参数边缘区域梯度波动大容易导致收敛到局部最优，优先选择optim默认的Nelder-Mead方法，无需计算梯度，对非平滑似然函数适配性更好。
• 带边界约束的梯度优化：如果要使用梯度类方法，选择L-BFGS-B方法，强制给方差参数设置正下界，避免优化过程中出现负方差导致似然计算报错。
• 先用专业SFA包做基准验证：可以调用frontier包的sfa()函数得到官方实现的估计结果，既可以验证自己编写的似然函数正确性，也可以将官方结果作为自定义似然优化的初始值。

二、合理初始值设置方案 #

• 系数初始值从OLS回归获取：先跑线性回归lm(Y ~ X1 + X2 + X3)，得到的截距和自变量系数直接作为b0、b1、b2、b3的初始值，远优于随机设置的初始值。
• 方差参数从OLS残差推导：计算OLS回归的残差标准差，按照7:3的比例拆分给SigmaV和Gamma作为初始值，符合你提供的真实参数比例，避免初始值偏离可行域过远。
• 网格搜索多组初始值：将6个参数的合理范围切分为2-3个水平，组合成多组初始值分别跑优化，选择最终似然值最大的那组结果，规避局部最优问题。

三、optim参数调整技巧 #

• 添加参数边界约束：使用L-BFGS-B方法时，添加lower = c(-Inf, -Inf, -Inf, -Inf, 1e-5, 1e-5)参数，保证SigmaV和Gamma始终为正，避免似然计算出现NaN。
• 调整收敛控制参数：在control列表中添加reltol = 1e-12缩小收敛容忍阈值，添加maxit = 10000提高最大迭代次数，让优化过程充分迭代直到真正收敛。
• 优化似然函数数值稳定性：你当前编写的for循环效率低且容易出现数值下溢，推荐修改为向量化版本，直接用dnorm、pnorm的对数参数提升稳定性：

LL <- function(theta){
  b0 <- theta[1]
  b1 <- theta[2]
  b2 <- theta[3]
  b3 <- theta[4]
  SigmaV <- theta[5]
  Gamma <- theta[6]
  # 非法参数直接返回负无穷，避免优化走偏
  if(SigmaV <= 0 | Gamma <= 0) return(-Inf)
  sigma_total <- sqrt(SigmaV^2 + Gamma^2)
  res <- Y - b0 - b1 * X1 - b2 * X2 - b3 * X3
  sum(
    log(2) - log(sigma_total) +
    dnorm(res / sigma_total, log = TRUE) +
    pnorm(-res * Gamma / (SigmaV * sigma_total), log.p = TRUE)
  )
}

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