La fórmula para las pruebas bayesianas A / B no tiene ningún sentido

Estoy usando la fórmula de la prueba ab bayesiana para calcular los resultados de la prueba AB utilizando la metodología bayesiana.

dónde

• $\alpha_A$ en uno más el número de éxitos para A
• $\beta_A$ en uno más el número de fallas para A
• $\alpha_B$ en uno más el número de éxitos para B
• $\beta_B$ en uno más el número de fallas para B
• $B$ es la función Beta

Datos de ejemplo:

control: 1000 trials with 78 successes
test: 1000 trials with 100 successes

Una prueba estándar de apoyo no bayesiano me da resultados significativos (p <10%):

prop.test(n=c(1000,1000), x=c(100,78), correct=F)

#   2-sample test for equality of proportions without continuity correction
# 
# data:  c(100, 78) out of c(1000, 1000)
# X-squared = 2.9847, df = 1, p-value = 0.08405
# alternative hypothesis: two.sided
# 95 percent confidence interval:
#  -0.0029398  0.0469398
# sample estimates:
# prop 1 prop 2 
#  0.100  0.078 

mientras que mi implementación de la fórmula de Bayes (usando las explicaciones en el enlace) me dio resultados muy extraños:

# success control+1
a_control <- 78+1
# failures control+1
b_control <- 1000-78+1
# success control+1
a_test <- 100+1
# failures control+1
b_test <- 1000-100+1

is_control_better <- 0
for (i in 0:(a_test-1) ) {
  is_control_better <- is_control_better+beta(a_control+i,b_control+b_test) / 
                       (b_test+i)*beta(1+i,b_test)*beta(a_control,b_control)

}

round(is_control_better, 4)
# [1] 0

eso significa que $P(TEST > CONTROL)$ es $0$ , lo que no tiene ningún sentido dados estos datos.

¿Alguien podría aclarar?

En el sitio que cita hay un aviso

La función beta produce números muy grandes, por lo que si obtiene valores infinitos en su programa, asegúrese de trabajar con logaritmos, como en el código anterior. La función log-beta de su biblioteca estándar será útil aquí.

entonces su implementación es incorrecta. A continuación proporciono el código corregido:

a_A <- 78+1
b_A <- 1000-78+1
a_B <- 100+1
b_B <- 1000-100+1

total <- 0

for (i in 0:(a_B-1) ) {
  total <- total + exp(lbeta(a_A+i, b_B+b_A)
                       - log(b_B+i)
                       - lbeta(1+i, b_B)
                       - lbeta(a_A, b_A))

}

Produce un total = 0.9576921, es decir "probabilidades de que B venza a A a largo plazo" (citando su enlace), lo que parece válido ya que B en su ejemplo tiene una mayor proporción. Por lo tanto, es no un p -valor, sino más bien una probabilidad de que B es mayor que A (usted no esperar que sea <0,05).

Puede ejecutar las simulaciones simples para verificar los resultados:

set.seed(123)

# does Binomial distributions with proportions
# from your data give similar estimates?

mean(rbinom(n, 1000, a_B/1000)>rbinom(n, 1000, a_A/1000))

# and does values simulated in a similar fashion to
# the model yield similar results?

fun2 <- function(n=1000) {
  pA <- rbeta(1, a_A, b_A)
  pB <- rbeta(1, a_B, b_B)
  mean(rbinom(n, 1000, pB) > rbinom(n, 1000, pA))
}

summary(replicate(1000, fun2(1000)))

En ambos casos la respuesta es sí.

En cuanto al código, tenga en cuenta que for loop es innecesario y generalmente hacen que las cosas sean más lentas en R, por lo que puede usar alternativamente vapplyun código más limpio y un poco más rápido:

fun <- function(i) exp(lbeta(a_A+i, b_B+b_A)
             - log(b_B+i)
             - lbeta(1+i, b_B)
             - lbeta(a_A, b_A))

sum(vapply(0:(a_B-1), fun, numeric(1)))