如何在R语言中求解给定方程中的sigma参数值？

已知条件与方程 #

给定核心方程：

prev = P12/(P11 + P12)

其中各变量的计算公式如下：

P11 <- exp(-(sigma+alpha_h)*t-(beta1_dot_h/beta2_h)*(exp(beta2_h*(x+t))-exp(beta2_h*x)))

# P12拆分为三部分计算
part1 <- exp((beta1_dot_h/beta2_h)*exp(beta2_h*x)*(1-(1+gama)*exp(beta2_h*t)) - (alpha_dd + (1 + gama)*alpha_h)*t - (beta1_dot_dd/beta2_dd)*exp(beta2_dd*(x+t)))

part2 <- exp(gama*(beta1_dot_h/beta2_h)*exp(beta2_h*(x+t))*(1 - beta2_h*(t/2)) + (beta1_dot_dd/beta2_dd)*exp(beta2_dd*(x+(t/2)))*(1 - beta2_dd*(t/2))) * sigma

part3 <- (exp(((gama*alpha_h) + alpha_dd + (gama*beta1_dot_h)*exp(beta2_h*(x+(t/2))) + beta1_dot_dd*exp(beta2_dd*(x+(t/2))))*t) - 1) / ((gama*alpha_h) + alpha_dd - sigma + (gama*beta1_dot_h) * exp(beta2_h*(x+(t/2))) + (beta1_dot_dd*exp(beta2_dd*(x+t/2))))

P12 <- part1 * part2 * part3

已知参数取值（除sigma外）：

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

求解思路 #

将原方程变形为：
prev*(P11 + P12) = P12
进一步整理得到等式：
prev*P11 = P12*(1 - prev)
构造目标函数：
f(sigma) = prev*P11 - P12*(1 - prev)
我们需要找到使得f(sigma) = 0的sigma值，这属于单变量方程求根问题，可使用R语言内置的uniroot()函数求解。

具体R代码实现 #

# 1. 定义已知参数
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

# 2. 定义目标函数：输入sigma，返回f(sigma)的计算值
target_func <- function(sigma) {
  # 计算P11
  P11 <- exp(-(sigma + alpha_h)*t - (beta1_dot_h/beta2_h)*(exp(beta2_h*(x+t)) - exp(beta2_h*x)))
  
  # 计算P12的三个组成部分
  part1 <- exp((beta1_dot_h/beta2_h)*exp(beta2_h*x)*(1 - (1+gama)*exp(beta2_h*t)) - (alpha_dd + (1+gama)*alpha_h)*t - (beta1_dot_dd/beta2_dd)*exp(beta2_dd*(x+t)))
  
  part2 <- exp(gama*(beta1_dot_h/beta2_h)*exp(beta2_h*(x+t))*(1 - beta2_h*(t/2)) + (beta1_dot_dd/beta2_dd)*exp(beta2_dd*(x+(t/2)))*(1 - beta2_dd*(t/2))) * sigma
  
  part3_numerator <- exp(((gama*alpha_h) + alpha_dd + (gama*beta1_dot_h)*exp(beta2_h*(x+(t/2))) + beta1_dot_dd*exp(beta2_dd*(x+(t/2))))*t) - 1
  part3_denominator <- (gama*alpha_h) + alpha_dd - sigma + (gama*beta1_dot_h)*exp(beta2_h*(x+(t/2))) + beta1_dot_dd*exp(beta2_dd*(x+t/2))
  part3 <- part3_numerator / part3_denominator
  
  P12 <- part1 * part2 * part3
  
  # 返回目标函数值
  return(prev*P11 - P12*(1 - prev))
}

# 3. 使用uniroot求解，指定sigma的搜索区间（可根据实际情况调整）
# 先尝试区间c(0, 1)，若报错可扩大或调整区间范围
result <- uniroot(target_func, interval = c(0, 1))

# 输出求解结果
cat("求解得到的sigma值为：", result$root, "\n")

注意事项 #

• 如果uniroot()报错提示区间内函数值符号不异号，说明需要调整搜索区间，比如尝试c(-1, 1)或其他合理范围。
• 由于gama=0，代码中含gama的项会自动简化为0，不影响计算结果。

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