はじめに

仕事で、ベイジアンABテストを利用するのが良さそうだなと思い、調べてみましたのでまとめてみます。

ベイジアンABテスト

通常のABテストの懸念点

頻度論に基づく仮説検定では、帰無仮説における母数 $\theta_0$ に対して得られた標本 $\theta_1$ のがどの程度起こりうるかを算出する。

また、A/Bテストを開始する前に有意水準 $\alpha$ 、検出力 $1-\beta$ 、効果量 $d=\frac{\mu_0-\mu_1}{\sigma}$ を予め定め、サンプルサイズを決定しておかなけらばならない。

それぞれの指標の決定は、少なからず統計の知識が必要であり、万人には受け入れられにくい。

そして、設定したサンプルサイズが集まるまで、検定を行うまで待たなければならない。

ビジネスの現場において通常のABテストを適切に運用しようとすると、上記のようなことが制約になることが少なくない。

ベイズ統計を用いた検定

ベイズ統計では $\theta$ 自体を確率変数とし、データxを得た上での $\theta$ の事後分布 $p(\theta | x)$ を次のベイズの定理を用いて推定する。

$\displaystyle{ p(\theta |x) = \frac{p(x|\theta)p(\theta)}{p(x)} }$

ここで、 $p(\theta | x)$ は $\theta$ を条件付けた $x$ の確率分布（尤度）、 $p(\theta)$ は事前分布、p(x)は周辺尤度である。

ベイジアンABテストでは、A群とB群それぞれのイベント発生率の事後分布を算出する。

そして、これら二つの分布を比較することで、A群とB群に差があるかどうかを決定する。

イベント発生率の推定

ここであるイベントが発生するかどうかを考えます。

ある群におけるイベントの発生は、ベルヌーイ分布で表現することができる。

$\displaystyle{ x \sim Bernoulli(\theta)=\theta^{x_i}(1-\theta)^{1-x_i} }$

この時の尤度は次にようになります。

$\displaystyle{ f(x | \theta)=\prod_{i=1}^N\theta^{x_i}(1-\theta)^{1-x_i} }$

事前分布としてベータ分布を用いる。

$\displaystyle{ p(\theta)=f(\theta | \alpha,\beta) = \frac{1}{B(\alpha,\beta)}\theta^{\alpha-1}(1-\theta)^{\beta-1} }$

ベータ分布分布はベルヌーイ分布の共役事前分布である。

共役事前分布とは、事前分布を元に事後分布を計算した際に、それらが一致する分布のこと。

この共役事前分布を用いることで、ベイズ更新が行えるようになる。

なお、ベータ分布のパタメータ $\alpha, \beta$ をぞれぞれ１とすると、一様分布 $Uniform(0,1)$ に一致する。

これらより事後分布を計算すると次のようにある。

$\displaystyle{ \begin{align*} p(\theta\mid x) &\propto p(x\mid\theta)p(\theta)\\ &\propto \left(\prod_{i=1}^N\theta^{x_i}(1-\theta)^{1-x_i}\right)\theta^{\alpha -1}(1-\theta)^{\beta-1}\\ &=\theta^{\sum^N_{i=1}x_i+\alpha-1}(1-\theta)^{\sum^N_{i=1}(1-x_i)+\beta-1}\\ &= \theta^{\alpha'-1}(1-\theta)^{\beta'-1} \end{align*} }$

事後分布は $Beta(\alpha',\beta')$ となる。なお、 $\alpha', \beta'$ は次のようになる。

$\displaystyle{ \alpha'=\sum_{i=1}^Nx_i + \alpha \\ \beta'=\sum_{i=1}^N(1-x_i) + \beta }$

事前分布の設定

事前分布を適切に設定することでサンプルサイズを少なくすることができる。

事前分布の設定の考え方について、こちらでは次にように書いていました。

We do not recommend using a prior distribution so strong that it overwhelms any data that is observed. It’s generally good practice to choose priors that are a bit weaker than what the historical data suggest.

確かに、テストを行う前に得られている事前情報を表現するのもであるため、確信度が高いものを設定するのは適切ではなさそう。

RでベイジアンA/Bテスト

Rで手軽にABテストを行うには、bayesABpackageが便利です。二つのあるパラメータを指定したベルヌーイ分布から得られたデータを取得します。

A_binom <- rbinom(2000, 1, .62)
B_binom <- rbinom(4000, 1, .60)

bayesTest関数を用いてベイジアンA/Bテストが実行できます。 distribution引数で検定のパターン(分布)を、prior引数には事前分布のパラメータを指定します。

AB1 <- bayesTest(A_binom, 
                 B_binom, 
                 priors = c('alpha' = 1, 'beta' = 1), 
                 n_samples = 1e5, 
                 distribution = 'bernoulli')

Bの分布のパラメータよりAのパラメータの方が大きい確率を確認する。

Aの方の比率が高いと言える可能性はありますが、確信度としては高くありません。

> summary(AB1)[2]$probability$Probability
[1] 0.68993

結果を可視化してみます。以下のようにすると3種類の結果が出力されます。

plot(AB1)

一つ目は、事前分布の形状です。 f:id:saltcooky:20200722012207p:plain

二つ目は、それぞれの群の推定された事後分布です。 f:id:saltcooky:20200722012220p:plain

三つ目は、(A-B)/Bのデータの分布(相対効果値のヒストグラム)です。 f:id:saltcooky:20200722012231p:plain

AとBに差があるとは主張できないと考えられます。また、差があることを検出するにはサンプルサイズが足りないかもしれません。

事前分布を確信度が高いものにしてみました。

AB1 <- bayesTest(A_binom,
                 B_binom,
                 priors = c('alpha' = 90, 'beta' = 80), 
                 n_samples = 1e5, 
                 distribution = 'bernoulli')

Aの方の比率が高いと言える可能性は低くなった。

AとBに差があることを主張するにはサンプルサイズが足りない可能性がある方向により傾いたことになる。

事前分布を適切に設定することにより、よりよい意思決定に繋がりそうです。(言ってることが怪しい...)

> summary(AB1)[2]$probability$Probability
[1] 0.59808

サンプルサイズと確率

サンプルサイズと推定されるパラメータ $\theta$ の変化を確認します。

set.seed(1234)
N <- 10000
A_all_binom <- rbinom(N, 1, .61)
B_all_binom <- rbinom(N, 1, .60)

data.frame(samplesize = seq(100, N, 100)) %>% 
  mutate(A = map(samplesize, ~mean(A_all_binom[1:(.)])) %>% unlist(),
         B = map(samplesize, ~mean(B_all_binom[1:(.)])) %>% unlist()) %>% 
  pivot_longer(col = -samplesize, names_to = "data", values_to = "rate" ) %>% 
  ggplot(aes(x = samplesize, y = rate, group = data, col = data))+
    #geom_point()+
    geom_line()+
    theme(legend.position=c(0.8,0.8))

f:id:saltcooky:20200724190602p:plain

次に、サンプルサイズと推定されたA群の $\theta$ の方が高い確率の変化も確認しました。

test_list <- 
  lapply(1:100, 
        A = A_all_binom, 
        B = B_all_binom, 
        FUN = loop_bayesTest) %>% 
  unlist()

data.frame(p = test_list, samplesize = 1:N) %>% 
  ggplot(aes(x=samplesize, y=p))+
    geom_line()+
    geom_point()

f:id:saltcooky:20200724190614p:plain

5000サンプルぐらい集まれば確信を持ってA群の方が良いとして、テストを停止することができそうです。

テストの停止基準

上記の状態では、明確なサンプルサイズに関する設定がない場合ため、いつテストを停止するかの決定が行いにくい。

調べてみると次にような記事がありました。

参考1：towardsdatascience.com

参考2：medium.com

ここでは、テストの停止を判断するための損失関数を紹介していました。

f:id:saltcooky:20200727000355p:plain

ここで、 $\alpha , \beta$ はA群とB群のメトリックである。

$x$ は選択する群を意味する。

$\alpha$ が $\beta$ より小さい時にAを選択すると、損失は $\beta -\alpha$ 、 $\alpha$ が $\beta$ より大きい時にAを選択すると損失は0となる。

$\alpha , \beta$ の真の値を観測できないため、事後密度 $f(\alpha , \beta )$ の予想損失を計算する。

f:id:saltcooky:20200727000747p:plain

これは精度が悪い選択肢を間違って選択した場合の損失を表現している。

AとBのどちらかの群の予測損失が閾値 $\epsilon$ を下回ると、実験を停止する。

この閾値は精度が悪い選択肢を間違って選択した場合の損失の許容度を表現している。

また、サンプルサイズと検出精度のトレードオフになってる。

損失のシミュレーション

この損失関数を利用したベイジアンABテストの停止がどのような挙動を示すのかを確認しました。

基本的にはこちらのコードをそのまま動かしただけです。

github.com

シミュレーションの設定は次にようにしました。

A群の真値：0.015
B群の真値：0.018
最大サンプルサイズ：20000
試行回数：500
停止基準： $\epsilon=0.0005$

sim1 = run.sims(a.cvr = 0.015, 
                b.cvr = 0.018, 
                max.obs = 10000,
                obs.increments = 100,
                n.simulations = 500)

sim1 = getloss(sim1)

はじめに、ある試行におけるサンプルサイズに対するA群とB群の損失の変化を確認しました。

threshold = 0.001

sim1 %>% 
  filter(n.sim==3) %>% 
  pivot_longer(cols = c("loss.a","loss.b"), names_to="group", values_to ="loss") %>% 
  ggplot(aes(x = checkin, y = loss, group=group, col=group))+
    geom_line()+
    xlab("サンプルサイズ")+
    geom_hline(yintercept = threshold, linetype = 2, aes(colour = "Threshold")) +
    theme_gray(base_family = "HiraKakuPro-W3")+
    theme(legend.position=c(0.8,0.85))

f:id:saltcooky:20200729020032p:plain

Bの損失が1000付近で $\epsilon=0.0005$ を切っていました。

次に、サンプルサイズが増えていったときにテストが停止する割合を確認してみました。

outcomes %>%
  ggplot(aes(x = checkin * 2)) + # Multiply by two because both variants had this many obs
    stat_ecdf() +
    ggtitle("テストの停止に関する経験分布関数") +
    labs(x = "サンプルサイズ", y = "テストを停止した割合")+
 　 theme_gray(base_family = "HiraKakuPro-W3")

f:id:saltcooky:20200729020127p:plain

サンプルサイズが5000程度で全ての試行は停止しています。

小さい効果量の割に少ないサンプルサイズでテストを停止することができているようです。

500回の試行のうちどれだけB群を採択してしたのかを確認しました。

outcomes <-
  sim1 %>% 
  mutate(threshold.met = (loss.a < threshold | loss.b < threshold)) %>%
  filter(threshold.met | checkin == max(checkin)) %>%
  group_by(n.sim) %>%
  filter(checkin == min(checkin)) %>%
  mutate(selection = case_when(
    loss.a < loss.b ~ "A",
    loss.b < loss.a ~ "B"
  )) 

outcomes %>%
  group_by(selection) %>%
  count(name = "count") %>%
  ungroup %>% 
  mutate(freq = count / sum(count)) %>%
  ggplot(aes(x = selection, fill = selection, y = freq, label = round(freq, 2))) +
    geom_bar(stat = "identity") +
    geom_text(vjust = -0.75, size = 4) +
    ggtitle("各群の採択割合") +
    labs(x = "", y = "採択割合") +
    coord_cartesian(ylim = c(0, 1)) +
    theme_gray(base_family = "HiraKakuPro-W3") +
    theme(legend.position = "none")