机器算法验证 - 贝叶斯 A/B 测试公式没有任何意义 - 吾爱随笔录

贝叶斯 A/B 测试公式没有任何意义

机器算法验证 r 贝叶斯测试

2022-03-02 17:04:15

我正在使用贝叶斯 ab 测试中的公式，以便使用贝叶斯方法计算 AB 测试的结果。

Pr (p_{B} > p_{A}) = \sum_{i = 0}^{α_{B} - 1} \frac{B (α_{A} + i, β_{B} + β_{A})}{(β_{B} + i) B (1 + i, β_{B}) B (α_{A}, β_{A})}

$\Pr(p_B > p_A) = \sum^{\alpha_B-1}_{i=0} \frac{B(\alpha_A+i,\beta_B+\beta_A)}{(\beta_B+i)B(1+i,\beta_B)B(\alpha_A, \beta_A)}$

在哪里

$\alpha_A$ 一加 A 的成功次数
$\beta_A$ 一加 A 的失败次数
$\alpha_B$ 一加 B 的成功次数
$\beta_B$ 一加 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

这意味着 $P(TEST > CONTROL)$ 是 $0$ ，考虑到这些数据，这没有任何意义。

有人可以澄清一下吗？

1个回答

在您引用的网站上有一个通知

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

其它你可能感兴趣的问题

上一篇是否应该在学生 T 检验中校正标准偏差？下一篇如何评估连续变量的比例风险假设