貝葉斯 A/B 測試公式沒有任何意義

我正在使用貝葉斯 ab 測試中的公式，以便使用貝葉斯方法計算 AB 測試的結果。

在哪裡

• 一加 A 的成功次數
• 一加 A 的失敗次數
• 一加 B 的成功次數
• 一加 B 的失敗次數
• 是Beta 函數

示例數據：

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

一個標準的非貝葉斯道具測試給了我顯著的結果（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 

而我對貝葉斯公式的實現（使用鏈接中的解釋）給了我非常奇怪的結果：

# 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

這意味著是，考慮到這些數據，這沒有任何意義。

有人可以澄清一下嗎？

在您引用的網站上有一個通知

beta 函數產生非常大的數字，因此如果您在程序中獲得無限值，請務必使用對數，如上面的代碼中所示。您的標準庫的 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 具有更大的比例。因此，它不是**p值，而是 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)))

引用自：https://stats.stackexchange.com/questions/141141