两个Beta随机变量差的分布实现问题

我需要分析两个比例p1、p2的差分（p1-p2）分布特性，两者均采用Beta先验，后验也为Beta分布。已实现两种解法：

• 暴力抽样法：小样本（8/10 vs 6/10）和大样本（80/100 vs 60/100）下均能正常运行
• 基于Pham-Gia和Turkkan(1993)论文的闭式解法：小样本结果与抽样法匹配，但大样本时计算差分PDF（y <- func_PD_vec(x)）报错，首个报错x值为-0.02，错误信息：

Error in integrate(A1.simple, 0, 1, a = a, b = b, b.prime = b.prime, c = c, : non-finite function value

推测是tolerance::F1函数计算Appell's F1超几何函数时的数值稳定性问题，求修复方案。

暴力抽样代码 #

## prior for both groups
a_pri <- 1
b_pri <- 1
## observed success and N for each group
g1_succ <- 8
g1_N <- 10
g2_succ <- 6
g2_N <- 10
## posterior parameters
a_g1 <- a_pri + g1_succ
b_g1 <- b_pri + g1_N - g1_succ

a_g2 <- a_pri + g2_succ
b_g2 <- b_pri + g2_N - g2_succ
## sampling from each group's posterior
samp_g1 <- rbeta(n = 10000,shape1 = a_g1,shape2 = b_g1)
samp_g2 <- rbeta(n = 10000,shape1 = a_g2,shape2 = b_g2)
## brute force the difference and visualization:
samp_diff <- samp_g1 - samp_g2

df_BF <- data.frame(samp_g1,samp_g2,s_diff= samp_diff) 
ggplot(df_BF)+geom_density(aes(x=s_diff))+xlim(-1,1)

闭式解法代码 #

AA1 <- a_g1;BB1 <- b_g1;AA2 <- a_g2;BB2 <- b_g2
## writing the function, only function of PP (the difference)
func_PD <- function(PP){  
  A <- beta(AA1, BB1) * beta(AA2, BB2)
  if(PP > 0 & PP <= 1){
    f_p <- beta(AA2,BB1)*(PP^(BB1+BB2-1))*((1-PP)^(AA2+BB1-1))*
      tolerance::F1(a = BB1, b = AA1+BB1+AA2+BB2-2,
                    b.prime = 1-AA1, c = BB1+AA2, x = 1-PP, y = 1-PP^2)/A
  }
  if(PP >= -1 & PP < 0){
    f_p <- beta(AA1,BB2)*((-PP)^(BB1+BB2-1))*((1+PP)^(AA1+BB2-1))*
      tolerance::F1(a = BB2, b = 1-AA2,
                    b.prime = AA1+BB1+AA2+BB2-2, c = AA1+BB2, x = 1-PP^2, y = 1+PP)/A
  }
  if(PP==0){
    if(AA1+AA2>1 & BB1+BB2>0){
      f_p <- beta(AA1+AA2-1,BB1+BB2-1)/A
    }
    if(!(AA1+AA2>1 & BB1+BB2>0)){
      f_p <- NA
    }
  }
  return(f_p)
}
## vectorizing it:
func_PD_vec <- Vectorize(func_PD)
## setting x values varying from -1 to 1:
x <- seq(-0.99,0.99,0.01)
## getting the pdf of difference:
y <- func_PD_vec(x)
## getting cdf by integrating:
df <- data.frame(p=x,pdf=y)
df$cdf <- NA
A_ch <- 0
for (i in 1:nrow(df)){
  if(A_ch > 0.995) {break}
  KK <- integrate(func_PD_vec, lower=-1, upper=df$p[i])
  df$cdf[i] <- KK$value
  A_ch <- df$cdf[i]
  # df$cdf[i] <- trapzfun(func_PD, a=-1,b= df$p[i])$value
}
df$cdf[which(is.na(df$cdf))] <- 1
## visualizing:
ggplot(df,aes(p,pdf))+geom_point()

大样本测试代码 #

g1_succ <- 80
g1_N <- 100
g2_succ <- 60
g2_N <- 100

1. 替换不稳定的Appell F1实现 #

tolerance::F1在大参数下数值稳定性差，改用hypergeo::AppellF1函数，其数值计算更鲁棒：

• 先安装加载包：install.packages("hypergeo"); library(hypergeo)
• 修改闭式解法中F1的调用，注意参数映射：hypergeo::AppellF1(a, b1, b2, c, x, y)对应原tolerance::F1的a、b、b.prime、c、x、y

修改后的核心函数示例：

func_PD <- function(PP){  
  A <- beta(AA1, BB1) * beta(AA2, BB2)
  if(PP > 0 & PP <= 1){
    f_p <- beta(AA2,BB1)*exp((BB1+BB2-1)*log(PP))*exp((AA2+BB1-1)*log(1-PP))*
      hypergeo::AppellF1(a = BB1, b1 = AA1+BB1+AA2+BB2-2,
                    b2 = 1-AA1, c = BB1+AA2, x = 1-PP, y = 1-PP^2)/A
  }
  if(PP >= -1 & PP < 0){
    f_p <- beta(AA1,BB2)*exp((BB1+BB2-1)*log(-PP))*exp((AA1+BB2-1)*log(1+PP))*
      hypergeo::AppellF1(a = BB2, b1 = 1-AA2,
                    b2 = AA1+BB1+AA2+BB2-2, c = AA1+BB2, x = 1-PP^2, y = 1+PP)/A
  }
  if(PP==0){
    f_p <- ifelse(AA1+AA2>1 & BB1+BB2>0, beta(AA1+AA2-1,BB1+BB2-1)/A, NA)
  }
  return(f_p)
}

2. 数值计算优化 #

• 用exp()+log()替换直接幂运算，避免大参数下的数值溢出（已在上述代码中实现）
• 对极端PP值（接近±1）增加边界处理，比如设置极小的epsilon阈值，避免log(0)或幂运算错误

3. 大样本替代方案：正态近似 #

当样本量足够大时，Beta后验可近似为正态分布，差分分布也可近似：

• 计算两个Beta后验的均值和方差：

  mu1 <- a_g1/(a_g1 + b_g1)
  var1 <- (a_g1*b_g1)/((a_g1+b_g1)^2*(a_g1+b_g1+1))
  mu2 <- a_g2/(a_g2 + b_g2)
  var2 <- (a_g2*b_g2)/((a_g2+b_g2)^2*(a_g2+b_g2+1))
  

• 差分的近似正态分布：N(mu1 - mu2, var1 + var2)
  该方法计算速度快，大样本下精度足够。

内容的提问来源于stack exchange，提问作者Farzin Shamloo