Formula untuk Bayesian A / B Testing tidak masuk akal

10

Saya menggunakan rumus dari pengujian ab Bayesian untuk menghitung hasil tes AB menggunakan metodologi Bayesian.

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

dimana

  • αA dalam satu ditambah jumlah keberhasilan untuk A
  • βA dalam satu ditambah jumlah kegagalan untuk A
  • αB dalam satu ditambah jumlah keberhasilan untuk B
  • βB dalam satu ditambah jumlah kegagalan untuk B
  • B adalah fungsi Beta

Contoh data:

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

Tes prop non Bayesian standar memberi saya hasil yang signifikan (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 

sementara implementasi formula Bayes saya (menggunakan penjelasan di tautan) memberi saya hasil yang sangat aneh:

# 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

itu berarti bahwa adalah , yang tidak masuk akal mengingat data ini.P(TEST>CONTROL)0

Bisakah seseorang mengklarifikasi?

Yehoshaphat Schellekens
sumber
A questiob analisis Bayesian dengan p-valuetag? Saya pikir orang Bayesia menolak untuk melakukan apa pun dengan nilai-p.
Dilip Sarwate
Kanan Anda! hanya berpikir bahwa itu akan menarik lebih banyak perhatian!
Yehoshaphat Schellekens
@YehoshaphatSchellekens jika itu alasan sebenarnya saya menghapus p-valuetag karena tidak terkait.
Tim
Tentu tidak masalah.
Yehoshaphat Schellekens

Jawaban:

17

Di situs yang Anda kutip ada pemberitahuan

Fungsi beta menghasilkan angka yang sangat besar, jadi jika Anda mendapatkan nilai tak terbatas dalam program Anda, pastikan untuk bekerja dengan logaritma, seperti pada kode di atas. Fungsi log-beta perpustakaan standar Anda akan berguna di sini.

jadi implementasi Anda salah. Di bawah ini saya berikan kode yang diperbaiki:

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

}

Ini menghasilkan total = 0,9576921, yaitu "peluang bahwa B akan mengalahkan A dalam jangka panjang" (mengutip tautan Anda) yang terdengar valid karena B dalam contoh Anda memiliki proporsi yang lebih besar. Jadi, ini bukan nilai- p melainkan probabilitas bahwa B lebih besar dari A (Anda tidak berharap itu menjadi <0,05).

Anda dapat menjalankan simulasi sederhana untuk memeriksa hasilnya:

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

Dalam kedua kasus jawabannya adalah ya.


Seperti tentang kode, perhatikan bahwa untuk loop tidak perlu dan umumnya mereka membuat hal-hal lebih lambat di R, sehingga Anda dapat menggunakannya sebagai alternatif vapplyuntuk kode pembersih dan sedikit lebih cepat:

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
sumber
Hmm ... Saya ingin tahu apakah Anda benar-benar diuji untuk kecepatan, karena vapplytidak lebih vectorozied daripada forloop, sebaliknya, mereka pada dasarnya sama. Jawaban yang bagus.
David Arenburg
1
C / C ++ / Fortan forloop == vektor; R forloop! = Vektor. Ini pada dasarnya definisi vektor.
David Arenburg
1
@YehoshaphatSchellekens poin dengan log bukan tentang perangkat lunak tertentu tetapi komputasi statistik umum. Dalam contoh di situs Anda mengutip kode julia disediakan - julia juga bahasa yang sangat baik untuk pemrograman statistik dan log masih digunakan.
Tim
2
Sebenarnya, saya baru saja mengajukan pertanyaan mengenai diskusi tepat yang kami lakukan di SO, saya mungkin perlu memikirkan kembali pendekatan saya ke vapplydepan. Saya harap saya akan mendapatkan jawaban yang bagus sekali dan untuk selamanya.
David Arenburg
2
Ok, setelah lama berpikir dan meringkas komentar dan jawaban yang saya terima pada pertanyaan saya, saya pikir saya datang dengan pemahaman umum tentang apa yang vapplysebenarnya. Lihat jawaban saya di sini
David Arenburg