R语言optimize函数未找到最优值问题排查

使用R语言的optimize函数时出现异常：两个高度相似的数据集t1和t2，针对t1，optimize能正常返回最优参数与极小的目标函数值；但针对t2，函数返回区间[0,10]的端点值9.999995，目标函数值为预设的高惩罚值1e10，而将t1的最优参数代入t2可得到更优的目标函数值。

待最小化函数与数据集 #

my_beta_fc <- function(data, par){
  
  # data is a dataframe of the ndays distribution with 2 columns:
  # x = number of days ranging from 0 to x+
  # f_x = reach (absolute)
  
  a <- par

  data <- data[,1:2]
  names(data) <- c("x", "f_x")
  
  n <- max(data$x)
  
  p0 <- data %>% 
    mutate(reach_pc = f_x / sum(f_x)) %>% 
    filter(x == 0) %>% 
    pull(reach_pc)
    
  xbar <- data %>% 
    summarise(mean = sum(f_x*x)/sum(f_x)) %>% 
    pull(mean)
  
  b <- a*(n-xbar)/xbar
    
  # we want to minimise the difference between observed and calculated p0:
  result <- abs(p0 - (gamma(n+b)*gamma(a+b))/(gamma(a+n+b)*gamma(b)))
  
  # if the evaluated function returns NaN, replace with a high penalty
  # to steer the optimization away from regions where the function returns NaN:
  if(is.nan(result)|is.na(result)){
    return(1e10)
  }
  
  return(result)
  
}

t1 <- data.frame(
  n_days = 0:7,
  reach = c(40979971, 2110778, 1126387, 729457, 541512, 346607, 236263, 198262)
) 

t2 <- data.frame(
  n_days = 0:7,
  reach = c(41233610, 2017354, 1063684, 694576, 518144, 330215, 223006, 188648)
) 

t1数据集的运行代码及结果 #

param_solution <- optimize(
  f = function(param) my_beta_fc(data = t1, param),
  interval = c(0,10),
  tol = 0.00001
)

结果：

> param_solution
$minimum
[1] 0.05495449

$objective
[1] 0.0000003038929

t2数据集的运行代码及结果 #

param_solution <- optimize(
  f = function(param) my_beta_fc(data = t2, param),
  interval = c(0,10),
  tol = 0.00001
)

结果：

> param_solution
$minimum
[1] 9.999995

$objective
[1] 10000000000

异常原因分析 #

1. Gamma函数数值溢出触发惩罚机制 #

t2的均值xbar约为0.173，略小于t1的0.179。当参数a增大时，b = a*(n-xbar)/xbar的增长速度比t1更快（t2的(n-xbar)/xbar≈39.45，t1≈38.09），导致Gamma函数的输入参数迅速变得极大，超出R的数值计算范围，产生NaN，触发代码中返回高惩罚值1e10的逻辑。

2. 黄金分割搜索的局限性 #

optimize采用黄金分割搜索法，它需要在区间内找到函数值变化的趋势。如果搜索过程中，大部分区间内的函数值都被替换为1e10，函数会误将区间端点判定为极小值点——因为在它的搜索路径中，端点的惩罚值和其他区域一致，甚至被当作“最优”选择。

3. 优化区间设置不合理 #

合适的a值应该接近t1的最优解（0.05左右），但当前优化区间是[0,10]，范围过大。当a偏离极小值区域时，很快触发Gamma溢出，optimize无法有效探索到真正的极小值区域，最终走到区间端点。

验证 #

将t1的最优参数代入t2的目标函数，可得到远小于1e10的结果：

my_beta_fc(t2, 0.05495449)
# 输出约为0.00000028（具体值可自行计算）

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