Korzystam z formuły z testu ab Bayesian , aby obliczyć wyniki testu AB z wykorzystaniem metody Bayesa.
gdzie
- w jednym plus liczba sukcesów dla A
- w jednym plus liczba awarii dla A.
- w jednym plus liczba sukcesów dla B
- w jednym plus liczba awarii dla 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 wynosi , co nie ma sensu, biorąc pod uwagę te dane.
Czy ktoś może to wyjaśnić?
p-value
tagiem? Myślałem, że Bayesianie nie chcą mieć nic wspólnego z wartościami p.p-value
tag, ponieważ nie jest on powiązany.Odpowiedzi:
Na stronie, którą zacytowałeś, jest powiadomienie
więc Twoja implementacja jest błędna. Poniżej podaję poprawiony kod:
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:
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ć
vapply
czystszego i nieco szybszego kodu:źródło
vapply
nie jest bardziej wektoryzowany niżfor
pętla, wręcz przeciwnie, są one zasadniczo takie same. Ładna odpowiedź.for
Pętle C / C ++ / Fortan == wektoryzowane;for
Pętla R ! = Wektoryzacja. Jest to w zasadzie definicja wektoryzacji.vapply
w przyszłości może być konieczne ponowne przemyślenie mojego podejścia . Mam nadzieję, że raz na zawsze otrzymam miłą odpowiedź.vapply
naprawdę jest. Zobacz moją odpowiedź tutaj