Formuła do testowania Bayesian A / B nie ma żadnego sensu

10

Korzystam z formuły z testu ab Bayesian , aby obliczyć wyniki testu AB z wykorzystaniem metody Bayesa.

Pr(pB>pA)=i=0αB1B(αA+i,βB+βA)(βB+i)B(1+i,βB)B(αA,βA)

gdzie

  • αA w jednym plus liczba sukcesów dla A
  • βA w jednym plus liczba awarii dla A.
  • αB w jednym plus liczba sukcesów dla B
  • βB w jednym plus liczba awarii dla B
  • B to funkcja Beta

Przykładowe dane:

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

Standardowy test podpory bayesowskiej daje mi znaczące wyniki (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 

podczas gdy moja implementacja formuły Bayesa (używając wyjaśnień w linku) dała mi bardzo dziwne wyniki:

# 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

oznacza to, że P(TEST>CONTROL) wynosi 0 , co nie ma sensu, biorąc pod uwagę te dane.

Czy ktoś może to wyjaśnić?

Yehoshaphat Schellekens
źródło
Pytanie z analizy bayesowskiej z p-valuetagiem? Myślałem, że Bayesianie nie chcą mieć nic wspólnego z wartościami p.
Dilip Sarwate,
Masz rację! właśnie pomyślałem, że przyciągnie to więcej uwagi!
Yehoshaphat Schellekens
@ YehoshaphatSchellekens, jeśli to był prawdziwy powód, dla którego usunęłem p-valuetag, ponieważ nie jest on powiązany.
Tim
Pewnie nie ma problemu.
Yehoshaphat Schellekens

Odpowiedzi:

17

Na stronie, którą zacytowałeś, jest powiadomienie

Funkcja beta generuje bardzo duże liczby, więc jeśli otrzymujesz nieskończone wartości w swoim programie, pamiętaj o pracy z logarytmami, jak w powyższym kodzie. Przyda się tutaj funkcja log-beta twojej biblioteki standardowej.

więc Twoja implementacja jest błędna. Poniżej podaję poprawiony kod:

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

}

Daje to wynik całkowity = 0,9576921, czyli „szanse, że B w dłuższej perspektywie pokona A” (cytując link), co brzmi poprawnie, ponieważ B w twoim przykładzie ma większą proporcję. Tak, to nie p -value lecz prawdopodobieństwo, że B jest większa niż A (robisz nie oczekiwać, że będzie <0,05).

Możesz uruchomić proste symulacje, aby sprawdzić wyniki:

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

W obu przypadkach odpowiedź brzmi „tak”.


Jeśli chodzi o kod, zauważ, że pętla jest niepotrzebna i generalnie spowalniają działanie w R, więc możesz alternatywnie użyć vapplyczystszego i nieco szybszego kodu:

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)))
Tim
źródło
Hmm ... Zastanawiam się, czy faktycznie przetestowałeś prędkość, ponieważ vapplynie jest bardziej wektoryzowany niż forpętla, wręcz przeciwnie, są one zasadniczo takie same. Ładna odpowiedź.
David Arenburg,
1
forPętle C / C ++ / Fortan == wektoryzowane; forPętla R ! = Wektoryzacja. Jest to w zasadzie definicja wektoryzacji.
David Arenburg,
1
@ YehoshaphatSchellekens punkt z logami nie dotyczy określonego oprogramowania, ale ogólnych obliczeń statystycznych. W podanym przykładzie na stronie podany jest kod Julii - Julia jest również bardzo dobrym językiem do programowania statystycznego i nadal używane są dzienniki.
Tim
2
Właściwie właśnie zadałem pytanie dotyczące tej dokładnej dyskusji, którą przeprowadziliśmy na temat SO, vapplyw przyszłości może być konieczne ponowne przemyślenie mojego podejścia . Mam nadzieję, że raz na zawsze otrzymam miłą odpowiedź.
David Arenburg,
2
Ok, po długim zastanowieniu i podsumowaniu komentarzy i odpowiedzi, które otrzymałem na moje pytanie, wydaje mi się, że doszedłem do ogólnego zrozumienia tego, co vapplynaprawdę jest. Zobacz moją odpowiedź tutaj
David Arenburg,