R语言二元正态分布数值积分求解满足概率条件的t1=t2阈值

因为约束为t1=t2，我们直接令公共值为t，将问题转化为单变量求根问题：找到t使得 P(X1>t) + P(X1<t, X2>t) = 0.05，具体实现如下：

1. 依赖包准备 #

我们需要用到mvtnorm包计算二元正态分布的联合概率，先安装加载：

# 安装包（仅需运行一次）
install.packages("mvtnorm")
# 加载包
library(mvtnorm)

2. 定义参数和目标函数 #

# 定义分布参数
mu <- c(0, 0)
cov_mat <- matrix(c(1, 1/2, 1/2, 1), nrow = 2)

# 定义目标函数：输入t，返回概率和与0.05的差值
target_fun <- function(t) {
  # 计算P(X1 > t)
  p1 <- 1 - pnorm(t, mean = 0, sd = 1)
  # 计算P(X1 < t, X2 > t)
  p2 <- pmvnorm(
    lower = c(-Inf, t),
    upper = c(t, Inf),
    mean = mu,
    sigma = cov_mat
  )
  # 返回差值，求根时要让这个值等于0
  return(p1 + p2 - 0.05)
}

3. 单变量求根求解t #

用R内置的uniroot函数在合理区间内求根，这里我们选区间c(1, 3)（常规单变量0.05上侧分位数约1.64，加入第二项后t会更大，区间覆盖合理范围即可）：

# 求根
solve_res <- uniroot(target_fun, interval = c(1, 3))
# 输出求解得到的t值（即t1=t2的取值）
t_result <- solve_res$root
print(t_result)

结果验证 #

可以代入计算验证结果是否符合要求：

# 验证概率和
p1_test <- 1 - pnorm(t_result)
p2_test <- pmvnorm(lower = c(-Inf, t_result), upper = c(t_result, Inf), mean = mu, sigma = cov_mat)
cat("概率和为:", p1_test + p2_test, "\n")

正常运行后得到的t约为1.8，代入计算的概率和会非常接近0.05。

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