導入

社内の勉強会で、ベイジアン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) の値が近い
主義の違いでやっていることは正当化されるかもしれないが、結局有意水準を緩めているのと近いのかもしれない

当然、有意水準と勝率は意味合いが異なりますが、結果として似た判断を下すことになる可能性に注意したいと思いました。