带随机目标函数与等式约束的R语言参数优化问题求助

一、先解决alabama的调用错误 #

你遇到的递归求值错误，是因为alabama默认会把...中的所有参数同时传给目标函数fn和约束函数heq，但两个函数需要的额外参数不匹配（compare_distributions需ref_dist，opt_constraint需theta），导致参数传递冲突。

修正调用方式，用alabama的args.fn和args.heq分别传递各自的额外参数：

library(alabama)

# 修正后的调用
alabama(par = c(0, 1, 0, 1),
        fn = compare_distributions,
        heq = opt_constraint,
        args.fn = list(ref_dist = household_income),
        args.heq = list(theta = 0.5))

不过即使修正了调用，由于目标函数是随机的（每次调用rlnorm生成不同样本，KS统计量会波动），这类确定性优化器（alabama、nloptr的大部分算法）很难收敛——优化器依赖目标函数值的一致性判断搜索方向，随机波动会干扰收敛逻辑。

二、适合随机目标函数的R优化包推荐 #

针对带噪声（随机）的黑箱优化问题，推荐以下鲁棒的启发式优化包：

1. GenSA（模拟退火） #

模拟退火算法对非凸、噪声目标函数兼容性好，不需要目标函数的梯度，适合你的场景。
示例代码：

library(GenSA)

# 把等式约束转化为惩罚项，加到目标函数中
obj_fun_with_penalty <- function(x, ref_dist, theta, penalty_weight = 100) {
  ks_stat <- compare_distributions(x, ref_dist)
  constraint_val <- opt_constraint(x, theta)
  # 用平方项惩罚约束违反
  return(ks_stat + penalty_weight * (constraint_val)^2)
}

# 设置参数边界（避免sigma为负）
lower <- c(-Inf, 0.01, -Inf, 0.01)
upper <- c(Inf, 5, Inf, 5)

# 运行优化（设置种子减少随机性）
set.seed(123)
result <- GenSA(par = c(0,1,0,1),
                fn = obj_fun_with_penalty,
                lower = lower,
                upper = upper,
                ref_dist = household_income,
                theta = 0.5)

# 查看结果
print(result$par)
print(result$value)

2. DEoptim（差分进化） #

基于种群的进化算法，适合黑箱优化，对噪声鲁棒，支持通过惩罚项处理约束。
示例代码：

library(DEoptim)

# 带惩罚项的目标函数
obj_fun_penalty <- function(x) {
  ks_stat <- compare_distributions(x, household_income)
  constraint_val <- opt_constraint(x, 0.5)
  return(ks_stat + 100 * (constraint_val)^2)
}

# 参数边界
lower <- c(-5, 0.01, -5, 0.01)
upper <- c(5, 5, 5, 5)

# 运行优化
set.seed(456)
result <- DEoptim(obj_fun_penalty, lower, upper,
                  control = DEoptim.control(itermax = 100, popsize = 50))

# 查看结果
print(result$optim$bestmem)
print(result$optim$bestval)

3. cmaes（协方差矩阵自适应进化策略） #

适合高维、非凸、噪声目标，自动调整搜索策略，无需手动设置过多参数。
示例代码：

library(cmaes)

# 带惩罚项的目标函数
obj_fun <- function(x) {
  ks_stat <- compare_distributions(x, household_income)
  constraint_val <- opt_constraint(x, 0.5)
  return(ks_stat + 100 * (constraint_val)^2)
}

# 运行优化
set.seed(789)
result <- cmaes(c(0,1,0,1), obj_fun,
                lower = c(-Inf, 0.01, -Inf, 0.01),
                upper = c(Inf, 5, Inf, 5))

# 查看结果
print(result$par)
print(result$value)

三、额外建议 #

• 目标函数的随机性会影响优化稳定性，建议增加模拟样本量，或多次模拟取KS统计量的平均值减少波动：

  compare_distributions <- function(x, ref_dist, n_sim = 200000, n_rep = 3) {
    ks_vals <- replicate(n_rep, {
      E <- rlnorm(n = n_sim, meanlog = x[1], sdlog = x[2])
      L <- rlnorm(n = n_sim, meanlog = x[3], sdlog = x[4])
      Y <- E + L
      unname(ks.test(Y, ref_dist)$statistic)
    })
    return(mean(ks_vals))
  }
  

• 你的约束函数推导正确：对数正态分布均值为exp(mu + sigma²/2)，E[L]/(E[E]+E[L])=θ等价于(mu_L + sigma_L²/2) - (mu_E + sigma_E²/2) = log(θ/(1-θ))。

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