如何定位导致optimr报错的初始值并规避L-BFGS-B有限值问题

问题描述 #

使用optimr最小化选择行为负对数似然函数时，出现如下错误：

Error in optim(par = par, fn = efn, gr = egr, lower = lower, upper = upper,  : 
  L-BFGS-B needs finite values of 'fn'

难以定位触发错误的初始值，也不清楚如何通过设置边界规避问题。

错误原因分析 #

1. 选择概率为0导致log(0)：当cp = p[c]为0时，log(cp)返回-Inf，求和后NegLL变成Inf，不符合L-BFGS-B对有限函数值的要求。
2. 参数边界设置不合理：phi_param = invlogit(params[2])的范围是(0,1)，但optimr中设置lower=c(-Inf,0)，当params[2]趋近于-Inf时，phi_param趋近于0，可能导致概率计算异常。
3. 全局变量依赖：函数中使用的N_blocks、num_trials_per_block、beta_test未在函数内定义或作为参数传入，可能引发运行时错误。
4. 示例数据语法错误：displayed_stim、a、r、tf、Block向量的拼接缺少逗号，导致无法生成正确的数据框。

修复步骤 #

1. 避免log(0) #

给选择概率添加极小值（如1e-10），确保log()输入为有限值：

cp <- p[c]
cp <- max(cp, 1e-10)  # 防止概率为0
choiceProb <- c(choiceProb, cp)

2. 合理设置参数边界 #

由于invlogit(x)的输出范围是(0,1)，给原始参数设置合理边界，确保转换后的参数在有效范围内：

• 学习率alpha_param应在(0,1)，对应params[1]的边界设为(-10, 10)（invlogit(10)≈1，invlogit(-10)≈0）
• 注意力权重phi_param应在(0,1)，对应params[2]的边界设为(-10, 10)
  调用optimr时：

lower = c(-10, -10), upper = c(10, 10)

3. 消除全局变量依赖 #

将N_blocks、num_trials_per_block、beta_test作为参数传入函数，或在函数内计算：

new_loglikelihood_ACL <- function(params, rep, val, beta_test) {
  N_blocks <- length(val)
  num_trials_per_block <- nrow(rep %>% filter(Block == 1))
  # ... 其余代码
}

4. 修复示例数据语法错误 #

在向量拼接处补全逗号，比如：

displayed_stim <- c("2,6,7","2,4,9","3,5,7","1,5,9","3,4,8","2,4,9","2,6,7","3,5,7","1,5,9","2,4,9","1,5,9","3,5,7","3,4,8",
"2,6,7","1,6,8","3,4,8","1,6,8","1,6,8","3,4,8","1,5,9","3,5,7","1,6,8","3,5,7","1,5,9","1,5,9","3,4,8",
# ... 其余元素
)

同理修复a、r、tf、Block向量的逗号问题。

5. 添加参数有效性检查 #

在函数开头检查转换后的参数是否在合理范围，避免无效计算：

alpha_param <- invlogit(params[1])
phi_param <- invlogit(params[2])

# 确保参数在有效范围内
if (alpha_param <= 1e-10 | alpha_param >= 1 - 1e-10 | phi_param <= 1e-10 | phi_param >= 1 - 1e-10) {
  return(1e10)  # 返回极大值，让优化器避开这些参数
}

完整修正代码 #

修正后的似然函数 #

new_loglikelihood_ACL <- function(params, rep, val, beta_test) {
  alpha_param <- invlogit(params[1])
  phi_param <- invlogit(params[2])
  
  # 参数有效性检查
  if (alpha_param <= 1e-10 | alpha_param >= 1 - 1e-10 | phi_param <= 1e-10 | phi_param >= 1 - 1e-10) {
    return(1e10)
  }
  
  choiceProb_Full <- c()
  N_blocks <- length(val)
  
  p3 <- 1 - phi_param
  
  for(Block_Num in 1:N_blocks){
    choiceProb <- c()
    Val_DF <- data.frame(val[[Block_Num]])
    
    Rep_DF <- rep %>% filter(Block == Block_Num)
    bs <- c(Rep_DF$displayed_stim)
    choice <- c(Rep_DF$a)
    outcome <- c(Rep_DF$r)
    tf <- c(Rep_DF$tf)
    num_trials_per_block <- nrow(Rep_DF)
    
    dimension_hint <- ifelse(tf[1] %in% c('yellow', 'blue', 'orange'), 'color', 'shape') 
    
    phi_color <- ifelse(dimension_hint == 'color', phi_param, p3)
    phi_shape <- ifelse(dimension_hint == 'shape', phi_param, p3)
    
    for(Trial_Num in 1:num_trials_per_block){
      
      b <- as.numeric(unlist(strsplit(bs[Trial_Num], ",")))
      bandits_df <- subset(Val_DF, stimuli_num %in% b)
      
      Q_bandits <- (bandits_df$v_shape * phi_shape) + (bandits_df$v_color * phi_color)
      
      p_num <- exp(beta_test * Q_bandits) 
      p <- p_num / sum(p_num) 
      c <- choice[Trial_Num]
      cp <- p[c]
      cp <- max(cp, 1e-10)  # 避免log(0)
      choiceProb <- c(choiceProb, cp)
      
      color_num  <- subset(Val_DF, features_color == bandits_df$features_color[c])$stimuli_num
      shape_num <- subset(Val_DF, features_shape == bandits_df$features_shape[c])$stimuli_num

      # 更新值
      Val_DF$v_shape[shape_num] <- bandits_df$v_shape[c] + alpha_param * (outcome[Trial_Num] - Q_bandits[c]) * phi_shape
      Val_DF$v_color[color_num] <- bandits_df$v_color[c] + alpha_param * (outcome[Trial_Num] - Q_bandits[c]) * phi_color
    }
    choiceProb_Full <- c(choiceProb_Full, choiceProb)
  }
  NegLL <- -1 * sum(log(choiceProb_Full))
  return(NegLL)
}

修正后的示例数据与optimr调用 #

create_val_df <- function(){
  Q_shapes = c(0.1666667,0.1666667,0.1666667)
  v_shape = rep(Q_shapes, each=3)
  Q_color = c(0.1666667,0.1666667,0.1666667)
  
  Q_stimuli = as.data.frame(v_shape)
  Q_stimuli$v_color = rep(Q_color, each=3)  # 修正长度匹配问题
  Q_stimuli$features = c("yellow, circle","yellow, oval","yellow, square",  
                         "blue, circle","blue, oval","blue, square","orange, circle",  
                         "orange, oval","orange, square")  

  Q_stimuli$features_color <- c("yellow","yellow","yellow","blue",
                                "blue","blue","orange","orange","orange")

  Q_stimuli$features_shape <-  c("circle","oval","square","circle",
                                 "oval","square","circle","oval","square")

  Q_stimuli$stimuli_num <- c(1:9)

  finalList <- list(Q_stimuli,Q_stimuli,Q_stimuli,Q_stimuli,Q_stimuli,Q_stimuli)
  for(i in 1:6) finalList[[i]]$Block <- i
  
  return(finalList)
}

values_data = create_val_df()

# 修复向量逗号问题
displayed_stim <- c("2,6,7","2,4,9","3,5,7","1,5,9","3,4,8","2,4,9","2,6,7","3,5,7","1,5,9","2,4,9","1,5,9","3,5,7","3,4,8",
"2,6,7","1,6,8","3,4,8","1,6,8","1,6,8","3,4,8","1,5,9","3,5,7","1,6,8","3,5,7","1,5,9","1,5,9","3,4,8",
"2,4,9","2,4,9","1,6,8","1,6,8","2,6,7","3,5,7","3,4,8","2,6,7","2,6,7","2,4,9","1,5,9","1,5,9","3,4,8",
"2,6,7","3,5,7","1,6,8","3,5,7","1,6,8","1,6,8","2,6,7","3,4,8","2,6,7","3,4,8","2,4,9","3,5,7","2,4,9",
"1,5,9","2,4,9","1,6,8","3,4,8","2,6,7","1,6,8","2,4,9","1,5,9","1,6,8","2,4,9","3,4,8","2,6,7","2,4,9",
"2,6,7","3,5,7","1,5,9","1,5,9","3,5,7","3,5,7","3,4,8","2,4,9","2,6,7","2,6,7","1,5,9","1,5,9","3,5,7",
"1,6,8","2,6,7","1,6,8","1,6,8","3,5,7","3,4,8","3,5,7","3,4,8","2,4,9","2,4,9","1,5,9","3,4,8","3,5,7",
"2,4,9","2,4,9","3,5,7","2,4,9","1,5,9","3,4,8","3,5,7","1,5,9","3,4,8","2,6,7","1,6,8","1,6,8","2,6,7",
"2,6,7","3,4,8","1,5,9","1,6,8")
a <- c(3,1,2,2,3,1,2,2,1,1,1,1,1,2,2,1,2,2,1,1,1,1,1,1,1,1,1,3,1,1,1,1,1,1,1,1,1,2,1,1,3,1,3,3,3,3,2,3,2,2,3,2,
1,2,1,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,3,1,1,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,3,3,1,2,1,2,3,2,2,3,1,3,3,1,
1,3,2,3)
r <- c(0,1,0,1,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,0,1,0,0,0,1,1,1,1,1,1,1,0,0,0,0,1,0,0,0,0,1,1,1,1,1,0,1,
1,1,0,1,0,1,0,1,0,0,1,1,1,1,1,1,0,1,1,1,0,1,0,1,0,1,1,0,1,1,0,0,0,1,1,1,1,1,0,0,1,1,1,0,1,1,1,1,0,0,1,1,
1,0,1,1)
tf <- c("square","square","square","square","square","square","square","square","square","square","square",
"square","square","square","square","square","square","square","yellow","yellow","yellow","yellow",
"yellow","yellow","yellow","yellow","yellow","yellow","yellow","yellow","yellow","yellow","yellow",
"yellow","yellow","yellow","circle","circle","circle","circle","circle","circle","circle","circle",
"circle","circle","circle","circle","circle","circle","circle","circle","circle","circle","orange",
"orange","orange","orange","orange","orange","orange","orange","orange","orange","orange","orange",
"orange","orange","orange","orange","orange","orange","blue","blue","blue","blue","blue",
"blue","blue","blue","blue","blue","blue","blue","blue","blue","blue","blue",
"blue","blue","oval","oval","oval","oval","oval","oval","oval","oval","oval",
"oval","oval","oval","oval","oval","oval","oval","oval","oval")
Block <- c(rep(1,18), rep(2,18), rep(3,18), rep(4,18), rep(5,18), rep(6,18))  # 简化Block生成，确保长度匹配

rep_data <- data.frame(displayed_stim,a,r,tf,Block)

# 设置beta_test（需根据实际情况调整）
beta_test <- 1

starting_alpha <- logit(runif(1,0.1,0.9))
starting_phi <- logit(runif(1,0.6,0.9))

optim_output <- optimx::optimr(
  par = c(starting_alpha, starting_phi),
  fn = new_loglikelihood_ACL,
  method = "L-BFGS-B",
  lower = c(-10, -10), 
  upper = c(10, 10),
  val = values_data,
  rep = rep_data,
  beta_test = beta_test,
  control = list(ndeps = c(1e-5, 1e-5), maxit = 10000, trace = 1)  # trace=1可追踪迭代过程
)

调试技巧 #

• 使用control=list(trace=1)查看迭代过程中的参数和函数值，定位触发错误的参数
• 在似然函数中添加print(params)和print(NegLL)，实时输出计算结果
• 测试固定参数值，验证函数是否返回有限值

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