導入
社内の勉強会で、ベイジアンABテストについて発表があり、気になったのでABテストとシミュレーション実験で比較してみます。
普通のABテストで新たな施策が勝ち傾向はあるけど有意差が出なかったから採用しないともったいない気持ちになりますよね。
ベイジアンABテストでは各群の分布が求まるため、分布を見て判断したり、勝率を計算して判断することができます。
たとえば、分布を見た結果「勝率80%なら採用しよう！」ということですね。
しかし、「これって有意水準を緩めていることと同じでは？」と疑問を持ったので、簡単なシミュレーション実験で、p値と勝率を並べてみようと思います。
すなわち「普通のABテストで勝ち傾向ではあるが有意差は出なかったとき、ベイジアンABテストはどうなるのか？」です。
ABテストについての統計知識は以前まとめたのでこちらをご覧ください：
実験コードや一部数式は折りたたんであるので、気になる方は展開してください。
ベイジアンABテストとは？
ベイジアンABテストは、ベイズ主義の立場で行うABテストです。
事前分布と得られたデータの尤度を組み合わせて、事後分布を計算します。
各群の事後分布が得られるため、「A群に対してB群が勝つ確率」「0.2から0.3に真の値が含まれる確率」など自由な計算ができます。
| ABテスト | ベイジアンABテスト | |
|---|---|---|
| 主義 | 頻度論 | ベイズ（分布） |
| 評価軸 | 実験の正当性の評価 | 勝つ確率の評価 |
| 得られる区間 | 信頼区間 | 信用区間 |
ベイジアンABテストでは事後分布を見て、施策を良しとするかの判断を行いますが、わかりやすいのは勝率なので、今回は勝率を確認します。
参考文献
-
19-3. 95％信頼区間のもつ意味
- 信頼区間の意味合いから「実験の正当性」を見ている感覚がわかりやすいです
定式化とシミュレーション実験
以下の2つの仮定でそれぞれ実験を行います。
| 事前分布 | データの生成分布 | 事後分布 |
|---|---|---|
| 正規分布 | 正規分布 | 正規分布 |
| ベータ分布 | ベルヌーイ分布 | ベータ分布 |
普通のABテストの有意水準は0.05くらいと思っておくことにします。
正規分布
定式化
サンプルが正規分布から生成されていると仮定：
p(x|\mu)=\mathrm{Norm}(x|\mu,\sigma^2)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left\{-\frac{(x-\mu)^2}{2\sigma^2}\right\}
- $\mu$: 平均値（未知）
- $\sigma^2$: 分散（既知とする）
パラメータの事前分布：
p(\mu)=\mathrm{Norm}(\mu | \mu_0,\sigma_0^2)=\frac{1}{\sqrt{2\pi\sigma_0^2}}\exp\left\{-\frac{(\mu-\mu_0)^2}{2\sigma_0^2}\right\}
- $\mu_0, \sigma_0^2$: 事前知識
事後分布（分散既知）：
$$
p(\mu|D) = \mathrm{Norm}\left(\frac{n\bar{x}/\sigma^2+\mu_0/\sigma_0^2}{n/\sigma^2+1/\sigma_0^2}, \frac{1}{n/\sigma^2+1/\sigma_0^2}\right)
$$
- $\mu_0, \sigma_0^2$: 事前分布のパラメータ
- $\sigma^2$: データの従う分散（既知設定）
- $\bar x$: データの平均
- $n$: データ数
事後分布の計算
\begin{align*}
p(\mu|D)&\propto p(\mu)p(D|\mu)\\
&=p(\mu)\prod_{x_i\in D} p(x_i|\mu)\\
&\propto
\exp\left\{-\frac{(\mu-\mu_0)^2}{2\sigma_0^2}
-\sum_i\frac{(x_i-\mu)^2}{2\sigma^2}\right\}\\
&\propto \exp\left\{-\frac{1}{2}\left((n/\sigma^2+1/\sigma_0^2)\mu^2-2(\sum_i x_i/\sigma^2+\mu_0/\sigma_0^2)\mu+\mathrm{const.}\right)\right\}\\
&\propto \exp\left\{-\frac{1}{2}(n/\sigma^2+1/\sigma_0^2)\left(\mu-\frac{\sum_i x_i/\sigma^2+\mu_0/\sigma_0^2}{n/\sigma^2+1/\sigma_0^2}\right)^2\right\}
\end{align*}
勝率の計算：
$x_A \sim N(\mu_A, \sigma_A^2), x_B \sim N(\mu_B, \sigma_B^2)$のとき、
$$
x_A - x_B \sim N(\mu_A - \mu_B, \sigma_A^2 + \sigma_B^2)
$$
$P(x_A\le x_B)=P(x_A-x_B\le 0)=\mathrm{cdf}(0)$
この分布のCDF（累積密度関数）の値の$x=0$地点の値（$-\inftyから0までの積分値$）を見れば、Bの勝率がわかります（1から引けばAの勝率）
実験
実験設定
| A群 | B群 | |
|---|---|---|
| サンプル数 | 200 | 190 |
| 真の平均値（未知） | 2 | 3 |
| 標準偏差（既知） | 4 | 4 |
| 事前分布の平均値（決める） | 2.5 | 2.5 |
| 事前分布の標準偏差（決める） | 4 | 4 |
この実験を何回か実施したときに出てくるパターンを4つ抽出して記載します。
当然、等確率で登場するわけではなく、実行結果3,4は起きにくいです。
コード
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, ttest_ind
# 実験設定
num_a = 200
num_b = 190
true_mu_a = 2
true_mu_b = 3
sigma = 4 # 標準偏差(既知)
data_a=np.random.normal(true_mu_a, sigma, size=num_a)
data_b=np.random.normal(true_mu_b, sigma, size=num_b)
mean_a = data_a.mean()
mean_b = data_b.mean()
print(f"{mean_a=}")
print(f"{mean_b=}")
# ABテスト
test=ttest_ind(data_a, data_b, equal_var=False, alternative="two-sided")
print(f"{test.pvalue=}(two-sided)")
test=ttest_ind(data_a, data_b, equal_var=False, alternative="less")
print(f"{test.pvalue=}(less)")
# ベイジアンAB
sigma0 = sigma
mu0 = 2.5
x = np.linspace(0, 5, 1000)
mu_a = (num_a*mean_a/sigma**2 + mu0/sigma0**2)/(num_a/sigma**2 + 1/sigma0**2)
var_a = 1/(num_a/sigma0**2 + 1/sigma0**2)
sigma_a = np.sqrt(var_a)
mu_b = (num_b*mean_b/sigma**2 + mu0/sigma0**2)/(num_b/sigma**2 + 1/sigma0**2)
var_b = 1/(num_b/sigma0**2 + 1/sigma0**2)
sigma_b = np.sqrt(var_b)
y_a = norm.pdf(x, mu_a, sigma_a)
y_b = norm.pdf(x, mu_b, sigma_b)
plt.plot(x, y_a, label="A")
plt.plot(x, y_b, "--", label="B")
plt.legend()
plt.savefig("normal.png")
# 勝率計算
mu_ab=mu_a-mu_b
sigma_ab=np.sqrt(var_a+var_b)
a_less_than_b_prob=norm.cdf(0, mu_ab, sigma_ab)
print(f"p(A<B) = {a_less_than_b_prob}")
# サンプリングによる勝率計算
simulation_size=100000
x_a = np.random.normal(mu_a, sigma_a, size=simulation_size)
x_b = np.random.normal(mu_b, sigma_b, size=simulation_size)
print(f"p(A<B) = {np.mean(x_a<x_b)}")
-
実行結果1（有意差出て正しく採択されるケース）
- サンプルの平均値
- A: 1.95946
- B: 3.31938
- p値
- 片側検定（A<B）：0.00044
- 両側検定（A≠B）：0.00089
- 事後分布
p(A<B) = 0.99959
- サンプルの平均値
-
実行結果2（勝ち傾向はあるが有意差出ないケース）
- サンプルの平均値
- A: 2.20669
- B: 2.72540
- p値
- 片側検定（A<B）：0.10507
- 両側検定（A≠B）：0.21014
- 事後分布
p(A<B) = 0.89916
- サンプルの平均値
-
実行結果3（負け傾向になったとき）
- サンプルの平均値
- A: 2.61990
- B: 2.42435
- p値
- 片側検定（A<B）：0.69309
- 両側検定（A≠B）：0.61382
- 事後分布
p(A<B) = 0.31513
- サンプルの平均値
-
実行結果4（大きく負けたとき）
- サンプルの平均値
- A: 2.73743
- B: 1.58086
- p値
- 片側検定（A<B）：0.99839
- 両側検定（A≠B）：0.00323
- 事後分布
p(A<B) = 0.00221
- サンプルの平均値
ベータ分布
定式化
サンプルがベルヌーイ分布から生成されると仮定：
$$
\mathrm{Ber}(x|\theta)=\theta^x(1-\theta)^{1-x}, ~(x=0,1)
$$
- $\theta$: 1が発生する確率（未知）
パラメータの事前分布：
$$
\mathrm{Beta}(\theta|\alpha, \beta)=\frac{\theta^{\alpha-1}(1-\theta)^{\beta-1}}{B(\alpha, \beta)}
$$
事後分布：
\begin{align*}
p(\theta)\prod_{r\in D} p(r | \theta)&=\mathrm{Beta}(\theta|\alpha, \beta)\prod_{r\in D}\mathrm{Ber}(r|\theta) \\&\propto \theta^{\alpha-1+n_{1}}(1-\theta)^{\beta-1+n_{0}}\\&\propto\mathrm{Beta}(\theta|\alpha+n_{1}, \beta+n_{0})
\end{align*}
- $n_0$: 0の個数
- $n_1$: 1の個数
勝率の計算：
$x_A \sim \mathrm{Beta}(\alpha_A, \beta_A), x_B \sim \mathrm{Beta}(\alpha_B, \beta_B)$のとき、$x_A - x_B$ の分布を求めたいが正規分布のようにきれいに求められなさそうなので、シミュレーションで求める
simulation_size=100000
x_a = np.random.beta(alpha_a, beta_a, size=simulation_size)
x_b = np.random.beta(alpha_b, beta_b, size=simulation_size)
print(f"p(A<B) = {np.mean(x_a<x_b)}")
実験
実験設定
| A群 | B群 | |
|---|---|---|
| サンプル数 | 200 | 190 |
| 真の平均値（未知） | 0.2 | 0.3 |
| 事前分布のalpha（決める） | 1 | 1 |
| 事前分布のbeta（決める） | 1 | 1 |
この実験を何回か実施したときに出てくるパターンを4つ抽出して記載します。
当然、等確率で登場するわけではなく、実行結果3,4は起きにくいです。
コード
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.stats.proportion import proportions_ztest
from scipy.stats import beta
num_a = 200
num_b = 190
true_mu_a = 0.2
true_mu_b = 0.3
data_a=np.random.binomial(n=1, p=true_mu_a, size=num_a)
data_b=np.random.binomial(n=1, p=true_mu_b, size=num_b)
mean_a = data_a.mean()
mean_b = data_b.mean()
print(f"{mean_a=}")
print(f"{mean_b=}")
sum_a = data_a.sum()
sum_b = data_b.sum()
# ABテスト
z_score, p_value = proportions_ztest(
[sum_a, sum_b], [num_a, num_b],
alternative="two-sided"
)
print(f"{p_value=}(two-sided)")
z_score, p_value = proportions_ztest(
[sum_a, sum_b], [num_a, num_b],
alternative="smaller"
)
print(f"{p_value=}(smaller)")
# ベイジアンAB
x = np.linspace(0, 0.5, 1000)
y_a = beta.pdf(x, 1+sum_a, 1+num_a-sum_a)
y_b = beta.pdf(x, 1+sum_b, 1+num_b-sum_b)
plt.plot(x, y_a, label="A")
plt.plot(x, y_b, "--", label="B")
plt.legend()
plt.savefig("beta.png")
# サンプリングによる勝率計算
simulation_size=100000
x_a = np.random.beta(1+sum_a, 1+num_a-sum_a, size=simulation_size)
x_b = np.random.beta(1+sum_b, 1+num_b-sum_b, size=simulation_size)
print(f"p(A<B) = {np.mean(x_a<x_b)}")
-
実行結果1（有意差出て正しく採択されるケース）
- サンプルの平均値
- A: 0.18000
- B: 0.31579
- p値
- 片側検定（A<B）：0.00093
- 両側検定（A≠B）：0.00186
- 事後分布
p(A<B) = 0.99898
- サンプルの平均値
-
実行結果2（勝ち傾向はあるが有意差出ないケース）
- サンプルの平均値
- A: 0.20500
- B: 0.24211
- p値
- 片側検定（A<B）：0.18949
- 両側検定（A≠B）：0.37897
- 事後分布
p(A<B) = 0.81010
- サンプルの平均値
-
実行結果3（負け傾向になったとき）
- サンプルの平均値
- A: 0.24500
- B: 0.22632
- p値
- 片側検定（A<B）：0.66800
- 両側検定（A≠B）：0.66399
- 事後分布
p(A<B) = 0.33603
- サンプルの平均値
-
実行結果4（大きく負けたとき）
- サンプルの平均値
- A: 0.29500
- B: 0.20526
- p値
- 片側検定（A<B）：0.97943
- 両側検定（A≠B）：0.04114
- 事後分布
p(A<B) = 0.02041
- サンプルの平均値
まとめと感想
| 実行結果 | p値 | p(A<B) |
|---|---|---|
| 正規分布（実行結果1） | 0.00044 | 0.99959 |
| 正規分布（実行結果2） | 0.10507 | 0.89916 |
| 正規分布（実行結果3） | 0.69309 | 0.31513 |
| 正規分布（実行結果4） | 0.99839 | 0.00221 |
| ベータ分布（実行結果1） | 0.00093 | 0.99898 |
| ベータ分布（実行結果2） | 0.18949 | 0.81010 |
| ベータ分布（実行結果3） | 0.66800 | 0.33603 |
| ベータ分布（実行結果4） | 0.97943 | 0.02041 |
- 1-p値 と p(A<B) の値が近い
- 主義の違いでやっていることは正当化されるかもしれないが、結局有意水準を緩めているのと近いのかもしれない
当然、有意水準と勝率は意味合いが異なりますが、結果として似た判断を下すことになる可能性に注意したいと思いました。