A/Bテストの評価をベイズ統計でやってみない？
遊びでA/Bテストの評価をベイズ統計でやってみたら、思いのほか面白かったので記事に残します。
用語の定義
コンバージョン
コンバージョンとは「Webサイト上で起きた最終的な成果」のことです。
具体的に何を意味するかはサイトの種類によっては様々です。
例えば、ECサイトでは商品の購入で、SNSでは会員登録などです。
コンバージョン率
コンバージョン率は「成果に繋がる最初の行動に対して実際に成果に繋がった割合」のことです。
ECサイトではある商品が購入された数をその商品が表示された数で割ったものがコンバージョン率となります。
（例えば、ある商品の表示数が100でその内で購入数が2ならばコンバージョン率は 2/100=0.02になります。）
A/Bテスト
A/Bテストとは複数の施策に対して良否をつけるためのテストです。
例えば、ウェブサイトのユーザビリティーを上げるためにデザインAとデザインBを同じ期間試してどちらが表示数が多いのか比べたり、
広告のパターンAとパターンBを出し分けてどちらのコンバージョン率が大きいか比べたりすることです。
今回は広告のA/Bテストを考えます。
広告の表示に対して購入が何件有ったのか（コンバージョン率）をA/Bテストします。
コンバージョンの確率分布
表示に対して購入が有ったか、無かったのかの二通りです。
それはBernoulli（ベルヌーイ）分布で表現できます。
購入が有った時x=1
購入が無かった時x=0
とすれば、がx=1になる時の確率なので
がコンバージョン率ということになります。
xは確率変数,はパラメータ（母数とも呼ばれる）と呼ばれています。
なぜベイズ統計を使うのか
割合の問題点
コンバージョン率は購入数÷表示数であることを説明しました。ここでパターンAのコンバージョン率を確かめる実験を二回行ったとします。
実験①: 表示数100, 購入数2 でした。コンバージョン率は2/100=0.02になります。
実験②: 表示数10000、購入数200でした。コンバージョン率は200/10000=0.02になります。
例えば実験② は実験①に比べて長い間表示してたので表示数も多かったとします。
どちらの方が信頼できますか？
実験①はデータが少なくてたまたまコンバージョン率が0.02になった様に思えますが、実験②は実験①よりは信頼できそうです。
しかし、両方とも同じコンバージョン率になります。
コンバージョン率の様な割合はこの信頼性を表現出来ていない事になります。
尤度と最尤法
ベイズ統計の前に最尤法という方法でパラメータの推定をしてみます。
尤度
確率分布はパラメータを固定した場合の確率変数の関数でしたが、同じ関数を確率変数を固定した場合のパラメータの関数と捉えた物を尤度と言います。
この場合の確率変数には観測データが入ります。観測データの表示数がNの場合に尤度は以下の様になります。
最尤法
尤度が最大値の時のパラメータが尤も（もっとも）適した値であるという方針のもとでパラメータを推定する方法です。
最大値の時の引数の値を求めるには微分が0のところを求めればいいのですが、計算を楽にするために尤度の対数を微分をします。（対数関数は単調増加関数なので引数が最大のときに対数関数の値も最大になるので問題ない。）
\begin{align*}
\frac{d\log L(\theta)}{d\theta}
&=\frac{d}{d\theta}\log\prod_{i=1}^N\theta^{x_i}(1-\theta)^{1-x_i}\\
&=\frac{d}{d\theta} \sum_{i=1}^N \left(x_i \log\theta +(1-x_i)\log(1-\theta)\right) \\
&=\sum_{i=1}^N \left(\frac{x_i}{\theta} - \frac{1-x_i}{1-\theta}\right) =0
\end{align*}
\begin{align*}
\sum_{i=1}^N \left(x_i(1-\theta) - (1-x_i)\theta \right) &=0\\
\sum_{i=1}^N(x_i -\theta) &=0\\
\sum_{i=1}^N x_i - N\theta &= 0\\
\theta = \frac{\sum_{i=1}^N x_i}{N}
\end{align*}
Nは表示数で、はコンバージョンが有ったときに1で無かったときは0なので、
の尤も適した値はコンバージョン率だとわかります。
ですので、最尤法を使っても割合の問題が出てきて信頼性を測るすべはありません。
ベイズ統計
ベイズ統計ではパラメータを１つの決まった値ではなくて確率分布として考えます。
こうする事で割合や最尤法では測れなかった信頼性が測れる事が後にわかります。
ベイズの定理
確率分布のベイズの定理は以下の様になります。
は事後分布と呼ばれます。これをパラメータ(コンバージョン率)の分布と捉えます。
は
に依存してないので今回はあまり気にする必要がありません。
は尤度でベルヌーイ分布と同じ形でしたね。
は事前分布と呼ばれています。
共役事前分布
事前分布の形を決定しないと事後分布の計算はできませんが形は謎です。
ですので、この分布は適当に決めてしまいましょう。
ただ、事前分布を適当に決めるにしても何らかの指針は欲しいものです。また適当に決めると計算できない場合がほとんどです。
そこで共役事前分布というものを考えます。この共役事前分布は事後分布が計算できるように選んだ事前分布デス。また、事前分布と事後分布は同じ形の確率分布になります。
どういうものがあるのかはwikiなどを見て欲しいのですが、尤度がベルヌーイ分布の場合は事前分布としてベータ分布を選べば良いことがわかっています。
ベータ分布
はベータ関数ですが、積分をして1とするための規格化係数なのであまり気にしなくていいです。
形は図1の様になります。
に依存して形が大きく変わっていることがわかります。
また、が共に１の時は一様分布になっているのもわかります。
事後分布の導出
共役事前分布のベータ分布を使って事後分布を計算して見ましょう。（事後分布もベータ分布になるはずですが。）
\begin{align*}
f(\theta\mid x) &\propto f(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*}
パラメータがのベータ分布になりました。
但し、
です。はコンバージョン数だけ増加して、
はコンバージョンされなかった数だけ増加しています。
事前分布のベータ分布のパラメータとして（一様分布になる）を選ぶことにしましょう。
事後分布のグラフ
必要なライブラリをpipで入れます。（pandasは描画には使いませんが後で使うので一緒にインストールします。）
pip install numpy scipy matplotlib pandas
事後分布を描画する関数を作ります。
import numpy as np import scipy.stats from matplotlib import pyplot as plt def show_beta(h, t): x = np.linspace(0, 1, 1000) plt.xlim(0, 1) plt.ylim(0, 5) plt.xlabel(r"$\theta$",fontsize=15) plt.ylabel(r"$f(\theta\mid x)$",fontsize=15) alpha = 1 beta = 1 alpha += h beta += t f = scipy.stats.beta(alpha, beta) m = f.mean() v = f.var() y = f.pdf(x) plt.title(r'$\alpha$=%d+1, $\beta$=%d+1 [mean=%.2f, var=%.7f]' % (h,t,m,v), fontsize=15) plt.plot(x, y) plt.show()
hにコンバージョン数、tに非コンバージョン数を入れると事後分布を描画する様にしました。
また、図のタイトルに平均(mean)と分散(var)も表示しています。
実験①：表示数100,購入数2の時を描画してみましょう。
show_beta(2,100-2)
出力
実験②：表示数10000,購入数200の時は
show_beta(200,10000-200)
出力
図2と図３を見比べると、図２の分布には幅があり不確定であることがわかります。一方で図３は幅がほとんどなく値がほぼ確定してることがわかります。この幅は分散に依存してるので、実際に分散（var）を見ると図３の方が全然小さいことがわかります。この様にパラメータ（コンバージョン率）を分布で表現すれば、その値の信頼性が分散という形で確認できることがわかりました。
これがベイズ統計を使う利点の一つです。
ベイジアンA/Bテストの実装
準備ができたので、ベイズ統計を使ってA/Bテストの結果を評価しましょう。
これはベイジアンA/Bテストと呼ばれています。
早速実装しまし。
コード
class BayesianAB: def __init__(self, params=None, size=2): """ params are a list of Beta distribution's initial params. e.g. [[alpha of A, beta of A], [alpha of B, beta of B]] """ if params is None: self.params = [] for _ in range(size): self.params.append([1,1]) self.size = len(self.params) self.data = [] for _ in range(size): self.data.append([0,0]) print("the number of comparison: ", self.size) self.sampling() def update(self, data, sampling=True): """ data are a list of pairs of impression and conversion. e.g. [[imp of A, conv of A], [imp of B, conv of B]] """ if self.size != len(data): print("No match of the size.") for p, current, new in zip(self.params, self.data, data): imp = new[0] conv = new[1] current[0] += imp current[1] += conv p[0] += conv p[1] += (imp - conv) if sampling: self.sampling() def mean_ver(self): """ return [(mean of A, variance of A), (mean of B, variance of B), ...] """ return [(posterior.mean(), posterior.var())for posterior in self.posterior_list] def sampling(self, n_samples=50000): print("num of samples: ", n_samples) self.posterior_list = [beta(*p).rvs(n_samples) for p in self.params] def show_beta(self, title="", save=False, labels=None): plt.figure(figsize=(10, 5)) plt.title("Posterior distribution "+ title) cmap = plt.get_cmap('jet') color_list= [] for i, posterior in enumerate(self.posterior_list): color =cmap(0.25*(i+1)) color_list.append(color) plt.hist(posterior, bins=100, histtype="stepfilled", normed=True, color=color, alpha=0.5) handles = [Rectangle((0,0),1,1,color=c, ec="k", alpha=0.5) for c in color_list] if labels is None: labels = [chr(65+i) for i in range(self.size)] # create A,B,... plt.legend(handles, labels) if save: plt.savefig("{}.png".format(title)) plt.show() def diff_prob(self, index_high, index_low): prob = (self.posterior_list[index_low] < self.posterior_list[index_high]).mean() if prob < 0.5: prob = 1 - prob return prob def show_metrics(self): print(self.data) def metrics(self, labels=None): if labels is None: labels = [chr(65+i) for i in range(self.size)] # create A,B,... return pd.DataFrame(self.data, index=labels, columns=["Impressions", "Conversions"])
使用例
初期化するときに事前分布のベータ分布のパラメータを設定できる様にしましたが、ほとんどの場合は一様分布で十分なのでそのまま初期化しましょう。
またsizeで比較する数を指定できます。
abtest = BayesianAB(size=2)
出力
the number of comparison: 2 num of samples: 50000
作成したクラスの使い方を説明するのですが、ここではダミーの行動データを作りましょう。
def dummy_genrator(conversion_rate, n_impression): return np.array(bernoulli(conversion_rate).rvs(n_impression)) a = dummy_genrator(0.015, 204) b = dummy_genrator(0.021, 200) imp_a = len(a) imp_b = len(b) conv_a = np.count_nonzero(a) conv_b = np.count_nonzero(b) print("A \n表示数: {}\n購入数: {}\n".format(imp_a, conv_a)) print("B \n表示数: {}\n購入数: {}\n".format(imp_b, conv_b))
出力
A 表示数: 204 購入数: 2 B 表示数: 200 購入数: 3
このデータを使って更新します。
abtest.update([[imp_a, conv_a], [imp_b, conv_b]])
これで更新されているので、平均と分散を確認すると、
abtest.mean_ver()
出力
[(0.014610260919694517, 7.003148406424335e-05), (0.019850764074203515, 9.568870563433339e-05)]
[(平均,分散),(平均,分散)]という形式で出力してるので
パターンAの分布の平均が0.014610260919694517, 分散が7.003148406424335e-05
パターンBの分布の平均が0.019850764074203515, 分散が9.568870563433339e-05
です。
図も描画して見ると
abtest.show_beta()
出力
updateの時の順番と同じインデックスを指定してやります。Aなら0,Bなら１です。
abtest.diff_prob(0,1)
出力
0.6682600000000001
完全に一致してる場合は0.5になります。ここでは0.668...なのでやはりあまり有意差がないということです。
更新したデータを見るにはshow_metricsを使います。
abtest.show_metrics()
出力
[[204, 2], [200, 3]]
[[Aの表示数,Aの購入数],[Bの表示数],[Bの表示数]]という形式で出力してくれます。
本番っぽい使い方
大体の使い方はわかったと思うので、本番っぽい使い方をしてみましょう。
今回はパターンは2パターンではなく３パターンの比較します。それぞれのパターンをA,B,Cと呼びます。
ある程度データが溜まった想定でダミーデータの生成します。
a = dummy_genrator(0.015, 1002) b = dummy_genrator(0.021, 1200) c = dummy_genrator(0.016, 1130) imp_a = len(a) imp_b = len(b) imp_c = len(c) conv_a = np.count_nonzero(a) conv_b = np.count_nonzero(b) conv_c = np.count_nonzero(c) print("A \n表示数: {}\n購入数: {}\n".format(imp_a, conv_a)) print("B \n表示数: {}\n購入数: {}\n".format(imp_b, conv_b)) print("C \n表示数: {}\n購入数: {}\n".format(imp_c, conv_c))
出力
A 表示数: 1002 購入数: 22 B 表示数: 1200 購入数: 22 C 表示数: 1130 購入数: 18
早速グラフをみてみましょう。
abtest = BayesianAB(size=3) #今回は3パターンなのでsizeは３ abtest.update([[imp_a, conv_a], [imp_b, conv_b],[imp_c, conv_c]]) abtest.show_beta()
出力
3つとも重なっている部分が大きいですね。
数値でも確認しましょう。
print(abtest.diff_prob(1,0)) # BとAの比較 print(abtest.diff_prob(1,2)) #BとCの比較 print(abtest.diff_prob(2,0)) #CとAの比較
出力
0.72712 0.66822 0.84508
有意差がわかりませんね。もっとデータが必要です。データが溜まるまで待ちましょう。
データが溜まった想定でダミーデータの生成します。
a = dummy_genrator(0.015, 25100) b = dummy_genrator(0.021, 22000) c = dummy_genrator(0.016, 23300) imp_a = len(a) imp_b = len(b) imp_c = len(c) conv_a = np.count_nonzero(a) conv_b = np.count_nonzero(b) conv_c = np.count_nonzero(c) print("A \n表示数: {}\n購入数: {}\n".format(imp_a, conv_a)) print("B \n表示数: {}\n購入数: {}\n".format(imp_b, conv_b)) print("C \n表示数: {}\n購入数: {}\n".format(imp_c, conv_c))
出力
A 表示数: 25100 購入数: 364 B 表示数: 22000 購入数: 469 C 表示数: 23300 購入数: 354
先ほど作ったBayesianABのインスタンスにデータの追加更新をしてやりましょう。
abtest.update([[imp_a, conv_a], [imp_b, conv_b],[imp_c, conv_c]])
abtest.show_beta()
print(abtest.mean_ver())
出力
[(0.014823136368344494, 5.586078465062482e-07), (0.0212059201285066, 8.839507467322899e-07), (0.01526854730312879, 6.223394152971891e-07)]
全ての分布が最初の時よりも幅が小さくなっていることがわかります。また、パターンBが明らかに離れているので有意差があることがわかります。
print(abtest.diff_prob(1,0)) # BとAの比較 print(abtest.diff_prob(1,2)) #BとCの比較 print(abtest.diff_prob(2,0)) #CとAの比較
出力
1.0 1.0 0.65888
A,CとBは100%違うことがわかります。そして、Bのコンバージョン率の平均0.021と一番高いのでこの広告を選ぶべきだとわかります。
しかし、一番効果がある広告を選ぶならこれでいいのですが、目的がA,B,Cのそれぞれの有意差を調べることならまだAとCの有意差ははっきりしていません。
ですので、もうちょっとデータを溜めましょう。
データが溜まった想定でダミーデータの生成します。
a = dummy_genrator(0.015, 251000) b = dummy_genrator(0.021, 220000) c = dummy_genrator(0.016, 233000) imp_a = len(a) imp_b = len(b) imp_c = len(c) conv_a = np.count_nonzero(a) conv_b = np.count_nonzero(b) conv_c = np.count_nonzero(c) print("A \n表示数: {}\n購入数: {}\n".format(imp_a, conv_a)) print("B \n表示数: {}\n購入数: {}\n".format(imp_b, conv_b)) print("C \n表示数: {}\n購入数: {}\n".format(imp_c, conv_c))
出力
A 表示数: 251000 購入数: 3808 B 表示数: 220000 購入数: 4544 C 表示数: 233000 購入数: 3811
追加更新しましょう
abtest.update([[imp_a, conv_a], [imp_b, conv_b],[imp_c, conv_c]])
abtest.show_beta()
print(abtest.mean_ver())
出力
[(0.015141034587728202, 5.332992221848848e-08), (0.0207078937798314, 8.243283715216886e-08), (0.016250559254538843, 6.199512187766625e-08)]
abtest.diff_prob(2,0)
出力
0.99964
Cの方が良いと言う事がわかりましたね。
このようにベイジアンA/Bテストでは予め有意差を測るためにデータ数がどの位必要なのかを知ってる必要はなくインクリメントに更新して評価することが出来きます。例えば、リアルタイムに更新をすれば有意差が出たと判断した時点ですぐにA/Bテストをやめて機会損失を最小限に抑える事ができます。
これもベイズ統計の強みです。
カイ二乗検定と比較
せっかくなのでカイ二乗検定と比較してみましょう。
from scipy.stats import chi2_contingency def chi2(imp_a, conv_a, imp_b, conv_b,title=""): print(title) metrics = [ [conv_a, imp_a - conv_a], [conv_b, imp_b - conv_b], ] print(metrics[0]) print(metrics[1]) x2, p, dof, expected = chi2_contingency(metrics) print("カイ二乗値: ", round(x2, 5)) print("p値: ", p) if p < 0.05: print("有意差あり。") else: print("有意差なし。") print("\n")
最初のデータのとき
imp_a_1 = 1002 conv_a_1 = 22 imp_b_1 = 1200 conv_b_1 = 22 imp_c_1 = 1130 conv_c_1 = 18 chi2(imp_a_1, conv_a_1, imp_b_1, conv_b_1,"###1nd A vs B ###") chi2(imp_a_1, conv_a_1, imp_c_1, conv_c_1,"###1nd A vs C ###") chi2(imp_c_1, conv_c_1, imp_b_1, conv_b_1,"###1nd B vs C ###")
出力
###1nd A vs B ### [22, 980] [22, 1178] カイ二乗値: 0.20435 p値: 0.6512354862810796 有意差なし。 ###1nd A vs C ### [22, 980] [18, 1112] カイ二乗値: 0.74604 p値: 0.38773258145798684 有意差なし。 ###1nd B vs C ### [18, 1112] [22, 1178] カイ二乗値: 0.08233 p値: 0.7741616454611604 有意差なし。
2回目にデータが溜まった時
imp_a_2 = 1002 + 25100 conv_a_2 = 22 + 364 imp_b_2 = 1200+ 22000 conv_b_2 = 22 + 469 imp_c_2 = 1130 + 23300 conv_c_2 = 18 + 354 chi2(imp_a_2, conv_a_2, imp_b_2,conv_b_2,"###2nd A vs B ###") chi2(imp_a_2, conv_a_2, imp_c_2,conv_c_2,"###2nd A vs C ###") chi2(imp_c_2, conv_c_2, imp_b_2,conv_b_2,"###2nd B vs C ###")
（ベイジアンA/Bテストでは追加更新してたので表示数と購入数は1回目のデータとの合計になる。）
出力
###2nd A vs B ### [386, 25716] [491, 22709] カイ二乗値: 28.2126 p値: 1.086949244640548e-07 有意差あり。 ###2nd A vs C ### [386, 25716] [372, 24058] カイ二乗値: 0.13625 p値: 0.7120340373703777 有意差なし。 ###2nd B vs C ### [372, 24058] [491, 22709] カイ二乗値: 23.24073 p値: 1.4293812107138311e-06 有意差あり。
A,BとCは有意差があって、AとCの有意差は分からない。ベイジアンA/Bテストと結論は同じ。
3回目のデータが溜まった時
imp_a_3 = 1002 + 25100 + 251000 conv_a_3 = 22 + 364 + 3803 imp_b_3 = 1200 + 22000 + 220000 conv_b_3 = 22 + 469 + 4544 imp_c_3 = 1130 + 23300 + 233000 conv_c_3 = 18 + 354 + 3811 chi2(imp_a_3, conv_a_3, imp_b_3,conv_b_3,"###3rd A vs B ###") chi2(imp_a_3, conv_a_3, imp_c_3,conv_c_3,"###3rd A vs C ###") chi2(imp_c_3, conv_c_3, imp_b_3,conv_b_3,"###3rd B vs C ###")
出力
###3rd A vs B ### [4189, 272913] [5035, 238165] カイ二乗値: 231.76377 p値: 2.4586705996516234e-52 有意差あり。 ###3rd A vs C ### [4189, 272913] [4183, 253247] カイ二乗値: 11.01697 p値: 0.0009028171800658043 有意差あり。 ###3rd B vs C ### [4183, 253247] [5035, 238165] カイ二乗値: 137.02089 p値: 1.1932289659422658e-31 有意差あり。
AとCに有意差がある事がわかりましたね。
有意差を調べると言う意味では結果は同じになりましたね。
