La formule du test bayésien A / B n'a aucun sens

J'utilise la formule de test ab bayésien afin de calculer les résultats du test AB en utilisant la méthodologie bayésienne.

où

• $\alpha_A$ en un plus le nombre de succès pour A
• $\beta_A$ en un plus le nombre d'échecs pour A
• $\alpha_B$ en un plus le nombre de succès pour B
• $\beta_B$ en un plus le nombre d'échecs pour B
• $B$ est la fonction bêta

Exemples de données:

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

Un test d'hélice non bayésien standard me donne des résultats significatifs (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 

alors que ma mise en œuvre de la formule de Bayes (en utilisant les explications du lien) m'a donné des résultats très étranges:

# 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

cela signifie que $P(TEST > CONTROL)$ est $0$ , ce qui n'a aucun sens compte tenu de ces données.

Quelqu'un pourrait-il clarifier?

Sur le site que vous citez il y a un avis

La fonction bêta produit de très grands nombres, donc si vous obtenez des valeurs infinies dans votre programme, assurez-vous de travailler avec des logarithmes, comme dans le code ci-dessus. La fonction log-beta de votre bibliothèque standard sera utile ici.

donc votre implémentation est fausse. Ci-dessous, je fournis le code corrigé:

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

}

Il génère un total = 0,9576921, c'est-à-dire "les chances que B batte A à long terme" (en citant votre lien) ce qui semble valide puisque B dans votre exemple a une plus grande proportion. Donc, ce n'est pas une valeur p mais plutôt une probabilité que B soit supérieur à A (vous ne vous attendez pas à ce qu'il soit <0,05).

Vous pouvez exécuter des simulations simples pour vérifier les résultats:

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

Dans les deux cas, la réponse est oui.

En ce qui concerne le code, notez que pour la boucle n'est pas nécessaire et généralement ils ralentissent les choses dans R, vous pouvez donc utiliser alternativement vapplypour un code plus propre et un peu plus rapide:

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