我有两个比例（P1，P2），我想表征P1-P2。假设每一个都有beta，则每个beta的后验分布也将是beta。我最初的问题是关于两个beta分布的差异的分布。我发现了一个溶液（蛮力）

Question

## 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)

这效果很好，这是我对Pham-Gia和Turkkan（1993）论文的实施：

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()

我在那里设置的数据（8/10 vs 6/10），一切都很好，两种方法都匹配了。但是当我将其更改为：

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

我重复所有内容，在此步骤中我会有一个错误：

## getting the pdf of difference:
y <- func_PD_vec(x)

i使用for循环重复了此步骤，以查找发生此错误的第一个X，并发现第一个错误在x = -0.02时发生，这是错误，我得到了：”集成中的错误（a1.simple，0，1，a = a，b = b，b.prime = b.prime，c = c，：：： 非有限函数值“

有一种方法可以解决吗？我以为是由F1函数引起的，它计算了Appell的F1超几何函数，但不确定如何修复它。

我试图寻找其他（也许更稳定或高效）的实现，但找不到。
thanks！

我也与您使用过的来源有类似的问题。

INSTEAD，我发现“ phase1b”软件包中包含的功能更稳定，请访问

https：//github.com/genentech/phase1b

功能DBETADIFF和PBETADIFF是您所需要的。

 #devtools::install_github("https://github.com/Genentech/phase1b/", force = TRUE)
library(phase1b)
library(ggplot2)
## prior for both groups
a_pri <- 1
b_pri <- 1
## observed success and N for each group
g1_succ <- 80
g1_N <- 100
g2_succ <- 60
g2_N <- 100
## 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

x <- seq(-0.99,0.99,0.01)
df <- data.frame(p=x,pdf=dbetadiff(x,parX=c(a_g2,b_g2),parY =c(a_g1,b_g1) ))
ggplot(df,aes(p,pdf))+geom_line()

for (i in 1:nrow(df)){
   df$cdf[i] <- pbetadiff(x[i],parX=c(a_g2,b_g2),parY =c(a_g1,b_g1) )
 }
ggplot(df,aes(p,cdf))+geom_line()

Answer 1