Формула для байесовского А / Б тестирования не имеет никакого смысла

10

Я использую формулу из байесовского ab-тестирования , чтобы вычислить результаты теста AB, используя байесовскую методологию.

Pr(пВ>пA)знак равноΣязнак равно0αВ-1В(αA+я,βВ+βA)(βВ+я)В(1+я,βВ)В(αA,βA)

где

  • αA в один плюс количество успехов для A
  • βA в один плюс количество сбоев для A
  • αВ в один плюс количество успехов для B
  • βВ в один плюс количество сбоев для B
  • В - бета-функция

Пример данных:

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

Стандартный не байесовский тест дает мне существенные результаты (р <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 

в то время как моя реализация формулы Байеса (используя объяснения в ссылке) дала мне очень странные результаты:

# 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

это означает, что п(TЕST>СОNTрОL) равно 0 , что не имеет никакого смысла, учитывая эти данные.

Может кто-нибудь уточнить?

Yehoshaphat Schellekens
источник
Байесовский анализ с p-valueтегом? Я думал, что байесовцы отказались иметь какое-либо отношение к р-значениям.
Дилип Сарвэйт
Твое право! Просто думал, что это привлечет больше внимания!
Yehoshaphat Schellekens
@YehoshaphatSchellekens, если это реальная причина, по которой я удаляю p-valueтег, поскольку он не связан.
Тим
Конечно, нет проблем.
Yehoshaphat Schellekens

Ответы:

17

На сайте вы цитируете есть уведомление

Бета-функция выдает очень большие числа, поэтому, если вы получаете бесконечные значения в вашей программе, обязательно работайте с логарифмами, как в коде выше. Здесь вам пригодится функция log-beta вашей стандартной библиотеки.

поэтому ваша реализация неверна. Ниже приведен исправленный код:

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

}

Он выводит итоговое значение = 0,9576921, то есть «вероятность того, что B превзойдет A в долгосрочной перспективе» (цитируя вашу ссылку), что звучит правильно, так как B в вашем примере имеет большую пропорцию. Таким образом, это не р -значение, а вероятность того, что B больше , то A (вы не ожидаете , что это будет <0,05).

Вы можете запустить простое моделирование, чтобы проверить результаты:

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

В обоих случаях ответ - да.


Что касается кода, обратите внимание, что цикл for не нужен, и, как правило, он замедляет работу в R, так что вы можете альтернативно использовать его vapplyдля более чистого и немного более быстрого кода:

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)))
Тим
источник
Хм ... Интересно, проверяли ли вы на скорость, потому что vapplyне более вектозированы, чем forпетли, напротив, они в основном одинаковы. Хороший ответ, хотя.
Дэвид Аренбург
1
forЦиклы C / C ++ / Fortan == векторизация; R forпетля! = Векторизация. Это в основном определение векторизации.
Дэвид Аренбург
1
@YehoshaphatSchellekens: с журналами дело не в определенном программном обеспечении, а в общих статистических вычислениях. В примере на сайте, который вы цитируете, указан код julia - julia также является очень хорошим языком для статистического программирования, и все еще используются журналы.
Тим
2
На самом деле, я только что задал вопрос относительно этой конкретной дискуссии о SO, мне может понадобиться переосмыслить свой подход vapplyв будущем. Я надеюсь, что я получу хороший ответ раз и навсегда.
Дэвид Аренбург
2
Хорошо, после долгих размышлений и обобщения комментариев и ответов, которые я получил на мой вопрос, я думаю, что пришел к некоторому общему пониманию того, что на vapplyсамом деле. Смотрите мой ответ здесь
Дэвид Аренбург