如何在R语言中求解含sigma的自定义方程？（uniroot/optimise无效）

已知自定义R函数sigmagm(sigma)，除sigma外所有变量参数均已给定，目标是找到使该函数值为0的sigma值。以下是具体解决步骤：

第一步：代入已知参数简化函数 #

给定参数中gama=0，可直接消去所有含gama的项，大幅简化计算逻辑，避免不必要的复杂运算和潜在错误。简化后的函数及参数定义如下：

# 定义所有已知参数
alpha_h <- -0.001188
beta1_h <- -5.85382
beta2_h <- 0.43633
alpha_dd <- -0.002006
beta1_dd <- -6.24563
beta2_dd <- 0.01495
beta1_dot_dd <- exp(beta1_dd)
beta1_dot_h <- exp(beta1_h)
prev <- 0.057115867
t <- 5
x <- 15
gama <- 0

# 简化后的sigmagm函数
sigmagm_simplified <- function(sigma) {
  # 预计算固定项，减少重复运算
  term1 <- (-sigma * t) - (alpha_h * t)
  term2 <- -(beta1_dot_h / beta2_h) * (exp(beta2_h*(x+t)) - exp(beta2_h*x))
  term3 <- -(beta1_dot_h / beta2_h) * exp(beta2_h*x) * (1 - exp(beta2_h*t))
  term4 <- (alpha_dd + alpha_h) * t
  term5 <- (beta1_dot_dd / beta2_dd) * exp(beta2_dd*(x+t))
  term6 <- -(beta1_dot_dd / beta2_dd) * exp(beta2_dd*(x + t/2)) * (1 - beta2_dd*(t/2))
  term7 <- -log(sigma)
  
  # 计算对数项内的表达式
  exp_arg <- alpha_dd - sigma + beta1_dot_dd * exp(beta2_dd*(x + t/2)) * t
  term8 <- -log(exp(exp_arg) - 1)
  term9 <- -log(alpha_dd - sigma + beta1_dot_dd * exp(beta2_dd*(x + t/2)))
  term10 <- -log((1 - prev)/prev)
  
  total <- term1 + term2 + term3 + term4 + term5 + term6 + term7 + term8 + term9 + term10
  return(total)
}

第二步：确定sigma的有效定义域 #

函数包含多个对数项，sigma必须满足：

• sigma > 0（log(sigma)要求输入为正）
• exp(exp_arg) - 1 > 0 → exp_arg > 0
• alpha_dd - sigma + beta1_dot_dd * exp(beta2_dd*(x + t/2)) > 0

先计算约束上限，确定sigma的取值范围：

# 计算两个约束的上限值
constraint1 <- alpha_dd + beta1_dot_dd * exp(beta2_dd*(x + t/2)) * t
constraint2 <- alpha_dd + beta1_dot_dd * exp(beta2_dd*(x + t/2))
sigma_upper <- min(constraint1, constraint2)
cat("sigma的有效区间：(0,", sigma_upper, ")\n")

第三步：使用uniroot求解方程根 #

uniroot要求输入区间两端的函数值符号相反，先测试区间内的点确认符号变化：

# 测试区间内的点，确认符号变化
cat("sigmagm_simplified(0.001):", sigmagm_simplified(0.001), "\n")
cat("sigmagm_simplified(0.0005):", sigmagm_simplified(0.0005), "\n")

根据测试结果选择合适区间，执行uniroot求解：

# 执行根求解
result <- uniroot(sigmagm_simplified, interval = c(0.0001, 0.001))
cat("求解得到的sigma值：", result$root, "\n")
cat("对应函数值：", result$f.root, "\n")

第四步：备选方案——使用optimise最小化函数绝对值 #

如果uniroot因区间选择问题失败，可尝试最小化|sigmagm_simplified(sigma)|，找到使函数最接近0的sigma：

# 最小化函数绝对值的平方，提升求解稳定性
opt_result <- optimise(function(s) abs(sigmagm_simplified(s))^2, interval = c(0.0001, sigma_upper))
cat("优化得到的sigma值：", opt_result$minimum, "\n")
cat("对应函数值平方：", opt_result$objective, "\n")

常见问题排查 #

• 未简化原函数：gama=0时大量项可直接消去，复杂函数易出现计算溢出或精度问题
• 区间选择错误：uniroot必须保证区间两端函数符号相反，需先测试多个点确认
• 定义域违规：sigma必须满足所有对数项输入为正，否则函数会返回NaN或Inf

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