基于后验密度积分求解90%概率上限的鲁棒方法问询

针对你遇到的网格法鲁棒性不足问题，推荐采用自适应数值积分+根求解的组合方案，完全避免手动设置区间和网格大小的困扰，同时保持一维问题的计算效率。

核心思路 #

• 利用R内置的integrate()函数处理后验密度的自适应积分，无需预先划定固定区间（支持无穷上限）
• 通过uniroot()函数求解CDF等于目标概率（90%）的根，自动定位所需分位数
• 可选：通过优化方法快速定位后验模态，辅助确定根求解的初始区间，进一步提升效率

完整实现代码 #

library(numDeriv) # 用于估计后验方差（可选）

# 生成样本
y <- rgamma(22, 2, 0.5)
n <- length(y)
tau <- 0.9
target_prob <- 0.9 # 对应1-alpha=0.9

# 定义未归一化后验密度函数
unnorm_post <- function(q) {
  mu <- mean((y < q) - tau)
  sigma <- max(var((y < q) - tau), 1e-8)
  exp(-0.5 * log(sigma) - (n / 2) * (mu^2 / sigma))
}

# 计算后验密度的归一化常数（积分从0到无穷）
Z <- integrate(unnorm_post, lower = 0, upper = Inf)$value

# 定义目标函数：CDF(q) - target_prob，我们需要找到该函数的根
target_func <- function(q) {
  if (q <= 0) return(-target_prob)
  integrate(unnorm_post, lower = 0, upper = q)$value / Z - target_prob
}

# 方法1：自动扩展区间找到根求解的上下界
q_upper <- max(y)
while (target_func(q_upper) < 0) {
  q_upper <- q_upper * 1.5 # 逐步扩大区间直到CDF超过目标概率
}

# 方法2：利用后验模态和标准差确定区间（更高效，需numDeriv包）
# mode_est <- optim(par = mean(y), fn = function(q) -log_posterior(q), 
#                   method = "Brent", lower = 0, upper = max(y)*2)$par
# hess <- hessian(function(q) -log_posterior(q), x = mode_est)
# post_sd <- sqrt(1/hess)
# q_upper <- mode_est + 3 * post_sd

# 求解根
result <- uniroot(target_func, lower = 0, upper = q_upper)
upper_quantile <- result$root

# 输出结果
cat("90%后验分位数：", upper_quantile, "\n")

方案优势 #

• 鲁棒性强：无需手动设置区间上限b，integrate()自动处理无穷区间，自适应调整积分精度
• 效率更高：自适应积分仅在密度较高的区域分配更多计算资源，避免网格法在低密度区域的无效计算
• 精度可控：integrate()和uniroot()都支持设置精度参数（如rel.tol），满足不同场景需求
• 通用性好：适用于任何一维后验分布，无需针对特定分布调整逻辑

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