有序Logit模型对数似然优化报错：vmmin初始值非有限值

问题情况 #

• 构建有序Logit模型时，用BFGS、L-BFGS-B、Newton-Raphson等方法优化对数似然函数，反复出现报错：Error in optim(par = ..., fn = ..., method = ..., ...) : L-BFGS-B needs finite values of 'fn'
• 试过多种初始参数（c(0,0,0,0)、c(1,-0.5,0.5,-0.5)），问题没解决
• 排查过数据集，没发现异常；尝试二分法和Newton-Raphson时，还出现对数概率含NaN的警告

问题代码 #

# 定义逻辑分布的对数似然函数
loglikelihood <- function(par, z) {
  # 从参数向量中提取参数
  beta <- par[1]
  alpha <- par[2:4]
  
  # 从数据框中提取变量
  poss <- z[, 1]
  quartile <- z[, 2]
  
  # 计算每个分位数的概率
  prob1 <- exp(alpha[1] - poss * beta) / (1 + exp(alpha[1] - poss * beta))
  prob2 <- exp(alpha[2] - poss * beta) / (1 + exp(alpha[2] - poss * beta)) - prob1
  prob3 <- exp(alpha[3] - poss * beta) / (1 + exp(alpha[3] - poss * beta)) - exp(alpha[2] - poss * beta) / (1 + exp(alpha[2] - poss * beta))
  prob4 <- 1 - exp(alpha[3] - poss * beta) / (1 + exp(alpha[3] - poss * beta))
  
  # 根据观测到的分位数组合概率
  probabilities <- ifelse(quartile == 1, prob1,
                          ifelse(quartile == 2, prob2,
                                 ifelse(quartile == 3, prob3, prob4)))
  
  # 通过求和概率的对数计算对数似然
  ll <- sum(log(probabilities))
  
  # 返回负对数似然（用于最小化）
  return(-ll)
}

# 用初始参数值测试对数似然函数
ll <- loglikelihood(par = c(0, 0, 0, 0), z = table)
print(ll)


p0 <- c(0, 0, 0, 0)
MLE <- optim(par = p0, fn = loglikelihood, method = "BFGS", hessian = TRUE, control = list(fnscale = -1), z = table)
print(MLE)

问题根源 #

1. 数值溢出问题：当alpha[i] - poss * beta的绝对值很大时，exp()会溢出成Inf，导致概率计算出现0、Inf或NaN，取对数后就会产生-Inf或NaN，让似然函数值变成非有限值
2. 阈值参数无约束：有序Logit要求阈值alpha必须满足alpha[1] < alpha[2] < alpha[3]，但当前代码没加这个约束，优化时阈值顺序可能颠倒，导致prob2、prob3变成负数，取对数后出现NaN
3. 似然函数逻辑冲突：代码返回的是负对数似然（要最小化），但调用optim时又加了fnscale = -1，相当于把目标改成最大化对数似然，逻辑混乱容易引发数值问题

解决步骤 #

1. 用稳定函数计算概率 #

用R内置的plogis()代替手动计算逻辑分布CDF，它内部做了数值稳定处理，能避免exp溢出：

prob1 <- plogis(alpha[1] - poss * beta)
prob2 <- plogis(alpha[2] - poss * beta) - prob1
prob3 <- plogis(alpha[3] - poss * beta) - plogis(alpha[2] - poss * beta)
prob4 <- 1 - plogis(alpha[3] - poss * beta)

2. 强制阈值参数的顺序约束 #

通过参数变换保证alpha[1] < alpha[2] < alpha[3]，比如用指数变换让增量为正：

loglikelihood <- function(par, z) {
  beta <- par[1]
  alpha1 <- par[2]
  # 用exp保证增量为正，确保alpha1 < alpha2 < alpha3
  delta1 <- exp(par[3])
  delta2 <- exp(par[4])
  alpha2 <- alpha1 + delta1
  alpha3 <- alpha2 + delta2
  
  poss <- z[, 1]
  quartile <- z[, 2]
  
  prob1 <- plogis(alpha1 - poss * beta)
  prob2 <- plogis(alpha2 - poss * beta) - prob1
  prob3 <- plogis(alpha3 - poss * beta) - plogis(alpha2 - poss * beta)
  prob4 <- 1 - plogis(alpha3 - poss * beta)
  
  # 用case_when替代嵌套ifelse，可读性更好
  probabilities <- dplyr::case_when(
    quartile == 1 ~ prob1,
    quartile == 2 ~ prob2,
    quartile == 3 ~ prob3,
    quartile == 4 ~ prob4
  )
  
  # 替换极小概率，避免log(0)产生-Inf
  probabilities <- pmax(probabilities, 1e-10)
  ll <- sum(log(probabilities))
  # 返回负对数似然，让optim最小化它
  return(-ll)
}

调用optim时去掉fnscale = -1，保持逻辑一致：

p0 <- c(0, 0, 0, 0)  # 初始参数对应alpha1=0, alpha2=1, alpha3=2
MLE <- optim(par = p0, fn = loglikelihood, method = "BFGS", hessian = TRUE, z = table)
print(MLE)

3. 处理极端概率值 #

在计算对数前，把概率值替换成极小正数（比如1e-10），避免log(0)的问题：

probabilities <- pmax(probabilities, 1e-10)

4. 先验证初始参数的似然值 #

优化前先单独运行初始参数的似然函数，确认结果是有限值：

ll <- loglikelihood(par = p0, z = table)
print(ll)

如果输出是有限数，再执行优化。

内容的提问来源于stack exchange，提问作者Mirela Cașcaval