R语言求解给定方程 估算待估参数pi_1的取值

已知参数定义 #

所有提前给定的参数取值如下（R代码形式）：

k_0 = 0.21
k_1 = 0.21
m = 52

alpha = 0.05
beta  = 0.2
pi_0 = 0.669
# pi_1 为待估参数
power <- 1-beta
cz <- 20
z_alpha <- qnorm(p= alpha/2, lower.tail=FALSE)
Z_beta <- qnorm(p= beta, lower.tail=FALSE)   

待解方程 #

其余参数固定时，待求解的方程为：

cz  <- 1 + ((z_alpha + Z_beta)^2)*((pi_0*(1-pi_0)/m+ pi_1*(1- pi_1)/m + (((k_0)^2)*((pi_0)^2) + ((k_1)^2)*((pi_1)^2)))/((pi_0 - pi_1)^2))

求解方法 #

这是典型的一元非线性方程求根问题，不需要手动推导复杂的解析解，直接用R内置的uniroot()单变量求根函数即可得到精确解，操作步骤如下：

• 首先把方程做移项处理，构造目标函数：函数输入为pi_1，输出为「方程计算得到的cz值」和「目标cz值20」的差，我们要找的就是输出值为0时对应的输入值
• 确定pi_1的合法搜索范围：pi_1作为概率参数取值范围在(0,1)之间，注意方程分母存在(pi_0 - pi_1)^2项，需要避开pi_1=pi_0=0.669的无意义点，因此把搜索区间拆为(0, 0.668)和(0.67, 1)两段分别查找
• 将目标函数和搜索区间传入uniroot()，设置合适的计算容差即可得到结果

可直接运行的完整求解代码如下：

# 赋值固定参数
k_0 = 0.21
k_1 = 0.21
m = 52
alpha = 0.05
beta  = 0.2
pi_0 = 0.669
cz_target <- 20
z_alpha <- qnorm(p= alpha/2, lower.tail=FALSE)
Z_beta <- qnorm(p= beta, lower.tail=FALSE)

# 构造求根目标函数
cal_cz_diff <- function(pi_1) {
  cz_calc <- 1 + ((z_alpha + Z_beta)^2) * (
    (
      pi_0*(1-pi_0)/m + pi_1*(1-pi_1)/m + (k_0^2 * pi_0^2 + k_1^2 * pi_1^2)
    ) / (pi_0 - pi_1)^2
  )
  return(cz_calc - cz_target)
}

# 分区间求根
res_lower <- uniroot(cal_cz_diff, interval = c(0.001, 0.668), tol = 1e-7)
res_higher <- uniroot(cal_cz_diff, interval = c(0.67, 0.999), tol = 1e-7)

# 打印结果
cat("小于pi_0的pi_1解：", res_lower$root, "\n")
cat("大于pi_0的pi_1解：", res_higher$root, "\n")

如果只需要近似估算，可以在0到1的区间内按固定步长（比如0.001）生成一系列pi_1的候选值，逐个代入方程计算cz值，找到计算结果最接近20的候选值就是近似解，步长越小估算精度越高。

内容的提问来源于stack exchange，提问作者Ahir Bhairav Orai