Ateam Commerce Techのアドベントカレンダー 15日目は @o93 が担当します！
はじめに
1年程前，書籍を参考にベイジアンABテストを実装してみたのですが，簡単にJupyter Notebook作っただけでそのまま社内のgitlabに放置していました...今回アドベントカレンダーを機会に記事にすることにしました．ベイズ統計に関してまだまだ理解が追い付いていない身ですが，事象の曖昧さを定量的に表現するアプローチとして，仕事を超えて自分の考え方に大きく影響を受けていたりします，，，
難しい数式などを省いて，できるだけ解りやすくエンジニアライクに説明してみたつもりです．興味ある方は是非ご覧ください！
使用している技術
- Python - プログラミング言語
- JupyterLab - プログラムの実行環境
- NumPy - ABテストのシミュレートと計算
- Matplotlib - グラフの描画
- Figma - 図描いただけ..
- ウェブ最適化ではじめる機械学習 - 参考書籍
記載していること
- 通常のABテストとベイジアンABテストの違い
- ベイジアンABテストの実装方法
- ベイジアンABテストの評価方法
記載していないこと
- JupterLabの環境構築，PythonやNumPyのプログラミング手法などには触れません
- PyMC3などのライブラリは用いません
- 通常のABテストについては詳しく説明しません
- ベイズ統計学に関する詳しい説明はしません，数式も用いません
注意点
難しい数式や言葉をなるべく用いずに記事を書いたことで，厳密な表現などで誤解を生んでしまっているかもしれません...明らかな間違いや我慢出来ない点はコメントにてご指摘頂けると幸いです．
通常のABテストとの違い
まず，一般的なABテストの流れを示します．
- 利用者がサイトを訪れる
- デザインA，Bをランダムで選んで表示
- 利用者がターゲットをクリックしたかを記録
- データが溜まったら集計して，AとBそれぞれのクリック率を算出
- Aのクリック率が
0.4，Bが0.6なので，デザインBを採用しよう！！
図を説明するとこのようになるかと思いますが，最後の意思決定デザインBを採用は正しい判断でしょうか．クリックのところを確認すると，5 + 5で合計10アクセスの結果から判断しています．また，それぞれの利用者の性別や年齢はバラバラですが，10人だと利用者の特徴が偏ってしまう可能性があります．
この状況で，サイト運営者の多くは「10アクセスじゃ判断できない．もう少しアクセスログが溜まるのを待とう..」という判断を下しているかもしれません．では100アクセスで5クリックの差があればデザインBだと言えるでしょうか，，，，これらの判断は運営者の過去の経験に基づく感覚的な判断だと言えます．
以下は通常のABテストでのクリック率の確からしさを表そうと頑張ってみた図です．
横軸はクリック率，縦軸はクリック率がその値になる確率です．アクセス数が幾つの時点であれ，クリック率はAとBそれぞれの2つの値だけが算出されて確定してしまいます．この問題に対してベイズ的アプローチで挑んだのが，ベイジアンABテストです．通常のABテストとの違いは以下です．
- A，Bそれぞれのクリック率について，一定の範囲内に収まる確率が算出できる
- AよりBのほうが優れている確率が算出できる
- ↑のどちらもアクセスが増えると本来の値に近くなっていく様子が見える
ABテストの目的であるどちらが優れているかの値だけでなく，その時点での値の曖昧さを，常に定量的に表し意思決定に使うことができます！では，次項でベイジアンABテストを実装して，上図がどのようになるか試してみます．
ベイジアンABテストの実装
なるべく文章を使って説明しているので，Pythonが良く分からなかったらコード読み飛ばしても何とかなる気がしています...Google CoLabにコードを貼り付ければそのまま動作するかと思います．
Import
必要なライブラリをインポートします．
import numpy as np
from matplotlib import pyplot as plt
データの準備
ABテストに必要なデータを準備します．
# クリック率の連続値（結果はクリック率の確率分布が得られるので，0〜1までのクリック率の配列を定義）
click_rates = np.linspace(0, 1, 1001)
print('クリック率:', click_rates)
# 都度更新される確率分布（まだアクセスログが無い状態の結果を事前分布と呼ぶらしい）
a = np.array([1 / len(click_rates) for _ in click_rates])
print('確率分布:', a)
クリック率: [0. 0.001 0.002 ... 0.998 0.999 1. ]
確率分布: [0.000999 0.000999 0.000999 ... 0.000999 0.000999 0.000999]
２つの配列を作っています．
- クリック率
click_rates：先程のグラフの横軸 - 確率分布
a: 先程のグラフの縦軸
aは，それぞれのクリック率がその値になる確率を保持しています．以下のように全てを足し合わせると必ず1.0になります．
np.sum(a)
1.0000000000000004
（誤差は気にせず先に進めます..）
１〜６の目があるサイコロのそれぞれの確率が1 / 6で，足し合わせると1になるのと同じです．
グラフ表示してみます．
plt.plot(click_rates, a, color='#FF4444')
plt.show()
まだアクセスログのデータを反映していないテスト前なので，一様分布で初期化されたままになっています．
ここで注意したいのは，このグラフに表示されている赤い線は，クリック率そのものではないということです．クリック率はX軸であり，Y軸は，「ABテストの結果としてAのクリック率が０〜１の間の何処かの値になるが，まだテスト前だから，1つの値に初期化しとこう」という状態のクリック率が取り得る値の確率分布です．最初私はこの部分が理解できなかったため，ゆっくり時間をかけて考えました．
Bも同様に配列を作ります．
b = np.array([1 / len(click_rates) for _ in click_rates])
これでデータの準備は完了です！
アクセス発生による確率分布の更新
利用者によるアクセスが発生したときに確率分布を更新するための関数を準備します．
# 尖度関数（結果を得るのに使用する関数）
likelihood = lambda r: click_rates if r else (1 - click_rates)
# 事後分布を求める関数（ABテストのアクセスログを入力して評価結果を更新する関数）
def posterior(r, prior):
lp = likelihood(r) * prior
return lp / lp.sum()
複雑な関数が出てきましたね...この辺りの説明は省略します．詳細は書籍を参照ください！
デザインBでは1回目のアクセスからいきなりクリックが発生しているので，Bから試してみましょう．
# Bへのアクセスとクリックが発生
b = posterior(1, b)
# AとBの確率分布を描画
plt.plot(click_rates, a, label='a', color='#FF4444')
plt.plot(click_rates, b, label='b', color='#4444FF')
plt.legend()
plt.show()
まだアクセスログを1件しか入れていない中，いきなり利用者がクリックしたため，クリック率100%になる確率が高いかも..と読み取れるような確率分布になりました．Aの確率分布はまだ更新していないのでそのままです．1番上の図に記載した通りのアクセスログを使用して，AとBを更新してみましょう．
# Aのアクセスとクリック
a = posterior(0, a)
a = posterior(0, a)
a = posterior(1, a)
a = posterior(0, a)
a = posterior(1, a)
# Bのアクセスとクリック（最初のアクセスはさっき更新したので，残りの４回）
b = posterior(0, b)
b = posterior(1, b)
b = posterior(1, b)
b = posterior(0, b)
# AとBの確率分布を描画
plt.plot(click_rates, a, label='a', color='#FF4444')
plt.plot(click_rates, b, label='b', color='#4444FF')
plt.legend()
plt.show()
AとBでグラフに違いが出てきましたが，確率分布の重なり具合を見ると，Bが優れているとは判断できないように思えます，，，このようにアクセスが増える度に，その時点での判断材料となる確からしさを更新していくことができます．
テスト結果の検証
このままアクセスログが増えていき，以下のような結果になったとします．
| デザイン | アクセス | クリック |
|---|---|---|
| A | 88 | 20 |
| B | 83 | 35 |
これまでは1アクセスずつ確率分布を更新していましたが，一気に確率分布を求める方法で，結果を検証してみます．
# ベータ分布
def betaf(alpha, beta):
numerator = click_rates ** (alpha - 1) * (1 - click_rates) ** (beta - 1)
return numerator / numerator.sum()
# 事後分布
def posterior2(a, N, prior):
return betaf(a + 1, N - a + 1)
# Aの確率分布
a = 1 / len(click_rates)
# アクセスログを入力
a = posterior2(20, 88, a)
# Bの確率分布
b = 1 / len(click_rates)
# アクセスログを入力
b = posterior2(35, 83, b)
# AとBの確率分布を描画
plt.plot(click_rates, a, label='a', color='#FF4444')
plt.plot(click_rates, b, label='b', color='#4444FF')
plt.legend()
plt.show()
クリック数とアクセス数を入れるだけで確率分布が表示できるようになりました．この結果をHDIという区間を表す統計量を使用して評価します．詳しい説明は省略しますが，確率分布の一定範囲を合計（積分）することで，クリック率がその範囲内に収まる確率を求めています．
# HDIを算出
def hmv(xs, ps, alpha=0.95):
xps = sorted(zip(xs, ps), key=lambda xp: xp[1], reverse=True)
xps = np.array(xps)
xs = xps[:, 0]
ps = xps[:, 1]
return np.sort(xs[np.cumsum(ps) <= alpha])
# グラフを描画
def plot_hm(name, rates, v, hm, c, ay=0):
plt.plot(rates, v, c=c, label=name)
plt.annotate(
'', xy=(hm.min(), ay), xytext=(hm.max(), ay),
arrowprops=dict(color=c, shrinkA=0, shrinkB=0, arrowstyle='<->', linewidth=1.5))
plt.annotate('%.3f' % hm.min(), xy=(hm.min(), ay), ha='right', va='bottom', c=c)
plt.annotate('%.3f' % hm.max(), xy=(hm.max(), ay), ha='left', va='bottom', c=c)
region = (hm.min() < rates) & (rates < hm.max())
plt.fill_between(rates[region], v[region], 0, alpha=0.2, color=c)
# AとBそれぞれのHDIを算出
a_hm = hmv(click_rates, a, alpha=0.9)
b_hm = hmv(click_rates, b, alpha=0.9)
# AとBのグラフを描画
plot_hm('a', click_rates, a, a_hm, 'tab:red', 0.0016)
plot_hm('b', click_rates, b, b_hm, 'tab:blue', 0.0022)
plt.tight_layout()
plt.legend()
plt.show()
算出したHDIから以下のことが分かります．
- Aのクリック率は90％の確率で，0.161〜0.305の範囲内に収まる
- Bのクリック率は90％の確率で，0.336〜0.510の範囲内に収まる
しきい値を90％としたときの定量的な検証結果から，デザインAよりデザインBのクリック率が高い，と判断することができます！通常のABテストのときに描いたグラフと比べると，アクセスによる曖昧さが定量的に見える化できているはずです．
AよりBのクリック率が高い確率
最後に，AよりBのクリック率が高い確率を定量的に算出してみます．まず，ABそれぞれベータ分布に従ったサンプルを作成します．
# ベータ分布に従うサンプルを１０００００ずつ作成
a = np.random.beta(20, 88, size=100000)
b = np.random.beta(35, 83, size=100000)
# AとBのヒストグラムを描画
plt.hist(a, bins=60, alpha=0.5, color='tab:red', label='a')
plt.hist(b, bins=60, alpha=0.5, color='tab:blue', label='b')
plt.legend()
plt.show()
今度は確率分布ではなく，上のアクセス数のパラメータに従うベータ分布のサンプルデータを作成しています．偏りが無いように，十分な数をランダムに確保します．前回のグラフと微妙にグラフの形が異なっている理由は，前回はAとBで合計の数が異なっていますが，今回はAとBの差分から確率変数を作るためにABで同じ数のサンプルを確保しているからです．
次にB - Aを行うことで差分の確率変数を作ります．
# B - A（AよりBのほうが大きい確率）
delta = b - a
plt.hist(delta, bins=60, alpha=0.5, color='tab:blue')
plt.legend()
plt.show()
確率変数deltaが正の値を取る割合を以下のようにして求めた結果が，上記グラフの面積＝確率に相当します．
print((delta > 0).mean())
0.97621
これがデザインAよりデザインBのほうがクリック率が高い確率になります．97.62％と定量的に表現できています．
また，グラフの0より小さい部分の面積と，0より大きい部分の面積を比較すると，0より大きい部分が圧倒的に大きいことから，視覚的にも分かりやすく表現できているかと思います．
追記_2022/1/8
おわりに
今回はベイジアンABテストという手法を，複雑な数式を用いることなくプログラムだけで説明してみました．ベイジアンABテストを使用した分析を行うことで，確からしさ・曖昧さを見える化したり，定量的な値を求めることで，より速く，適切な意思決定を行うことができます．Google Optimizeなどの内部でも使用されている技術です．
個人的には，「事実」として確定値を点で算出する通常のアプローチから，1つの事実から「解釈」として確率分布を線で推論するベイズ的アプローチに思考の幅が広がる感覚を覚えました．まだ理解できていないことばかりですが，因果推論や量子力学など，未だ中二病の私にはたまらない妄想の題材の1つになっており，，今後も妄想と勉強を続けていこうと思います！
Ateam Commerce Techのアドベントカレンダー 16日目は @NamedPython が担当します！
ここまでお読み頂き，ありがとうございました！