如何解决MLE调用中出现的'Error in optim: non-finite finite-difference value [1]'错误

问题重现 #

运行以下最小可复现代码：

library(SuppDists)
library(stats4)

mle(function(x) -log(dPearson(x, N=21, rho=0.999)), 0.5)

触发错误：

Error in optim(start, f, method = method, hessian = TRUE, ...) :
non-finite finite-difference value [1]

将rho改为0.99时无错误，但尝试指定method = "L-BFGS-B", lower = 0.00001, upper = 0.99999时，出现新错误：

L-BFGS-B needs finite values of 'fn'

实际使用场景代码：

library(SuppDists)
library(stats4)

test_rMZ_2rDZ <- function(rMZ, nMZ, rDZ, nDZ) {
  # Run MLE
  H1_MZ_MLE = mle(function(x) -log(dPearson(x, N=nMZ, rho=rMZ)), 0.5)
  H1_DZ_MLE = mle(function(x) -log(dPearson(x, N=nDZ, rho=rDZ)), 0.5)
  
  H0_MLE = mle(function(x) -(log(dPearson(x, N=nMZ, rho=rMZ)) + log(dPearson(x/2, N=nDZ, rho=rDZ))), 0.5)
  
  # Extract log-likelihoods
  H0_logLik = logLik(H0_MLE)
  H1_logLik = logLik(H1_MZ_MLE) + logLik(H1_DZ_MLE)
  
  # Run likelihood ratio test 
  test_stat = 2 * (H1_logLik - H0_logLik)
  pchisq(test_stat, df = 1, lower.tail = FALSE)
}

test_rMZ_2rDZ(rMZ = 0.999, nMZ = 21, rDZ = 0.98, nDZ = 52)

错误原因 #

当rho接近1时，dPearson计算出的概率密度值可能趋近于0，导致log(dPearson)返回-Inf，取负号后变成Inf（非有限值），优化算法无法处理。即使添加参数上下限，优化过程中仍可能采样到使目标函数非有限的x值，触发L-BFGS-B的报错。

解决方案 #

1. 目标函数添加数值稳定性处理 #

在计算-log(dPearson)前，先检查dPearson的结果，若过小则替换为极小正数，避免log(0)的情况。

2. 选择更合理的初始值 #

当rho接近1时，Pearson相关系数的MLE应接近rho，将初始值设为rho附近而非固定0.5，减少优化过程中遇到极端值的概率。

3. 严格限制参数范围并使用鲁棒优化方法 #

指定L-BFGS-B方法，同时设置参数上下限为远离0和1的极小/极大值（如1e-6和1-1e-6），避免参数接近边界导致的数值奇异。

修改后的实际场景代码 #

library(SuppDists)
library(stats4)

test_rMZ_2rDZ <- function(rMZ, nMZ, rDZ, nDZ) {
  # 定义带数值稳定处理的负对数似然函数
  neg_log_lik <- function(x, N, rho) {
    dp <- dPearson(x, N = N, rho = rho)
    # 避免log(0)，将极小值替换为1e-16
    dp <- ifelse(dp < 1e-16, 1e-16, dp)
    -log(dp)
  }
  
  # 定义H0的负对数似然函数
  neg_log_lik_H0 <- function(x, nMZ, rMZ, nDZ, rDZ) {
    dp_mz <- dPearson(x, N = nMZ, rho = rMZ)
    dp_dz <- dPearson(x/2, N = nDZ, rho = rDZ)
    dp_mz <- ifelse(dp_mz < 1e-16, 1e-16, dp_mz)
    dp_dz <- ifelse(dp_dz < 1e-16, 1e-16, dp_dz)
    -(log(dp_mz) + log(dp_dz))
  }
  
  # Run MLE，使用更合理的初始值和鲁棒优化方法
  H1_MZ_MLE <- mle(neg_log_lik, start = list(x = rMZ), 
                   method = "L-BFGS-B", lower = 1e-6, upper = 1 - 1e-6,
                   N = nMZ, rho = rMZ)
  H1_DZ_MLE <- mle(neg_log_lik, start = list(x = rDZ), 
                   method = "L-BFGS-B", lower = 1e-6, upper = 1 - 1e-6,
                   N = nDZ, rho = rDZ)
  
  H0_MLE <- mle(neg_log_lik_H0, start = list(x = (rMZ + rDZ)/2), 
                method = "L-BFGS-B", lower = 1e-6, upper = 2*(1 - 1e-6),
                nMZ = nMZ, rMZ = rMZ, nDZ = nDZ, rDZ = rDZ)
  
  # Extract log-likelihoods
  H0_logLik <- logLik(H0_MLE)
  H1_logLik <- logLik(H1_MZ_MLE) + logLik(H1_DZ_MLE)
  
  # Run likelihood ratio test 
  test_stat <- 2 * (H1_logLik - H0_logLik)
  pchisq(test_stat, df = 1, lower.tail = FALSE)
}

# 测试修改后的函数
test_rMZ_2rDZ(rMZ = 0.999, nMZ = 21, rDZ = 0.98, nDZ = 52)

说明 #

• 数值稳定处理通过ifelse将极小的dPearson值替换为1e-16，确保log(dp)返回有限值。
• 初始值设置为rho（或rho的均值），让优化起点更接近真实MLE，减少迭代过程中遇到极端值的概率。
• H0的参数上限设为2*(1-1e-6)，因为x/2需要小于1，所以x最大接近2。

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