求Beta分布95% HPD区间算法及dbeta逆函数实现方法

一、计算Beta(7,25)的95% HPD区间 #

方法1：使用现成工具包 #

最便捷的方式是借助HDInterval包的hdi()函数，它直接支持Beta分布的HPD区间计算：

# 安装并加载包
install.packages("HDInterval")
library(HDInterval)

# 计算95% HPD区间
hdi(dbeta, shape1 = 7, shape2 = 25, credMass = 0.95)

方法2：手动实现（适配你的代码思路） #

如果想手动实现，核心逻辑是找到一个密度阈值c，让密度≥c的所有p值对应的累积概率为95%，且这个区间是最短的（HPD的核心特性）。可以通过数值优化实现：

# 定义目标函数：给定阈值c，计算密度≥c区间的累积概率，让其等于0.95
target <- function(c, shape1, shape2, credMass) {
  # 找到密度等于c的两个临界点（Beta(7,25)是单峰分布，存在两个解）
  lower <- uniroot(function(x) dbeta(x, shape1, shape2) - c, interval = c(0, qbeta(0.5, shape1, shape2)))$root
  upper <- uniroot(function(x) dbeta(x, shape1, shape2) - c, interval = c(qbeta(0.5, shape1, shape2), 1))$root
  # 计算区间内的累积概率
  prob <- pbeta(upper, shape1, shape2) - pbeta(lower, shape1, shape2)
  return(prob - credMass)
}

# 求解最优阈值c
c_opt <- uniroot(target, interval = c(0, dbeta(qbeta(0.5,7,25),7,25)), shape1=7, shape2=25, credMass=0.95)$root

# 计算HPD区间的上下限
hpd_lower <- uniroot(function(x) dbeta(x,7,25)-c_opt, interval = c(0, qbeta(0.5,7,25)))$root
hpd_upper <- uniroot(function(x) dbeta(x,7,25)-c_opt, interval = c(qbeta(0.5,7,25),1))$root

# 输出结果
cat("95% HPD区间：(", round(hpd_lower,4), ", ", round(hpd_upper,4), ")\n", sep="")

二、根据密度值反推p值（dbeta的"逆函数"实现） #

R中没有直接的dbeta逆函数（因为密度函数不是单调的），但对于Beta(7,25)这种单峰分布，给定一个≤峰值的密度值，会对应两个p值，可以用uniroot()分别求解：

# 定义函数：输入密度值，返回对应的两个p值
dbeta_inv <- function(d, shape1, shape2) {
  peak_pos <- qbeta(0.5, shape1, shape2) # 分布峰值的位置
  max_density <- dbeta(peak_pos, shape1, shape2) # 分布的最大密度值
  
  if(d > max_density) stop("输入的密度值超过该Beta分布的峰值，无解")
  if(d == max_density) return(c(peak_pos, peak_pos))
  
  # 求解左侧小于峰值的p值
  lower_p <- uniroot(function(x) dbeta(x, shape1, shape2) - d, interval = c(0, peak_pos))$root
  # 求解右侧大于峰值的p值
  upper_p <- uniroot(function(x) dbeta(x, shape1, shape2) - d, interval = c(peak_pos, 1))$root
  
  return(c(lower_p, upper_p))
}

# 示例：查找密度值为2时对应的p值
dbeta_inv(2, 7, 25)

注意：如果输入的密度值大于分布峰值，函数会报错，因为不存在对应的p值。

内容的提问来源于stack exchange，提问作者Juan Kim