こんにちは。 データチームの後藤です。
A/Bテストはサービス改善のための施策の効果測定に欠かせないツールですが、最近のVASILYでは、運用するサービスが増えてきたことに伴いA/Bテストの内容も多様化してきました。今回はそのA/Bテストにベイズ推論を用いた具体的な例を紹介します。
問題設定
あるサービスのコンバージョン率を上げるため、コンバージョンの前提となる行動Xを増やすための修正を実施しました。ここで、「コンバージョン」は商品の購入などの成約を、「コンバージョン率」は利用数に対して成約に結びついた割合を、「行動X」は買い物カゴに商品を入れるなどコンバージョンの前提となる行動のことを指すことにします。修正がコンバージョン率の上昇に寄与したのかをデータから判断する必要があります。
修正後のページ(パターンA)を表示したグループと、修正前のページ(パターンB)を表示したグループの行動ログを集計し、以下の3つの仮説を検証します。各パターンでは、表示したページ以外の条件は同じであると仮定します。
仮説
- パターンAはパターンBより行動Xを取りやすい
- パターンAはパターンBよりコンバージョン率が高い
- パターンAはパターンBよりコンバージョンに至るまでの時間が短い
これらの仮説を検証する際、通常は平均値の差の検定や母比率の差の検定などの仮説検定が用いられます。今回は、ベイズ推論のフレームワークを利用しビジネスサイドが状況を判断しやすい結論を出していきます。
ベイズ推論を使うメリット
ベイズ推論を利用した場合、従来の仮説検定と比べて以下のメリットが考えられます。
- 前提知識や常識を事前分布として反映できる
- 推定値の分布に正規性を仮定していないので、より柔軟な区間推定ができる
- データの生成過程をモデリングすることにより、得られた結果からストーリーを語ることができる
PyMC3とは
PyMC3はPythonでベイズ推論を実行できるフレームワークです。 http://docs.pymc.io/index.html
公式ドキュメントのExampleをみるとわかるように、様々な機能が実装されています。今回のようなA/Bテストに用いるだけでなく、GLMやMixture Modelの推論など様々な用途に利用できます。今回はPyMC3の生成モデルの定義とMCMCによるサンプリングの機能を利用します。
仮説1: パターンAはパターンBより行動Xを取りやすい
現実的には、ユーザーはサービスの利用開始から様々な過程を通じて行動Xに至る(あるいは、至らない)と考えられます。問題を単純にするために、それらを一旦無視して、利用開始から24時間以内に行動Xを1回でも取ったかどうかを焦点とします。
行動Xを取ったユーザーを1、行動Xを取らなかったユーザーを0と数値化することにより、これらのサンプルが確率のベルヌーイ分布から生成されると考えます。またに関しての情報はないので、事前分布として一様分布を採用します。
データの集計結果は以下のようになりました。
- グループA - サンプル数: 1602 - 行動Xを取った数: 740 - 平均: 0.462 - グループB - サンプル数: 2121 - 行動Xを取った数: 933 - 平均: 0.440
グループAのほうが行動Xを取った割合が高いので、修正の効果が出ているかもしれません。
PyMC3のフレームワークを使えば、今回のモデルを以下のように記述できます。行動Xを取る確率が増えたかどうかも知りたいので、というとの差を新たな生成量として定義します。
モデル
import pymc3 as pm with pm.Model() as model: # 事前分布として一様分布を採用 p_A = pm.Uniform('p_A', 0, 1.0) p_B= pm.Uniform('p_B', 0, 1.0) # 各行動はベルヌーイ分布から生成される obs_A = pm.Bernoulli("obs_A", p_A, observed=sample_A) obs_B = pm.Bernoulli("obs_B", p_B, observed=sample_B) # 新たな生成量として行動を起こす確率の差を定義 delta_prob = pm.Deterministic('delta_prob', p_B - p_A)
推論
with model: start = pm.find_MAP(fmin=optimize.fmin_powell) step = pm.Slice() trace = pm.sample(20000, step=step, start=start) _ = pm.traceplot(trace) _ = pm.plot_posterior(trace)
、とその差の標本は均衡分布に収束していることが確認できます。各変数の分布も左右対称のベル型をとっていることもわかります。
pm.summary(trace)
mean | sd | mc_error | hpd_2.5 | hpd_97.5 | n_eff | Rhat | |
---|---|---|---|---|---|---|---|
p_A | 0.461861 | 0.012475 | 0.000059 | 0.437458 | 0.485777 | 39737.0 | 0.999980 |
p_B | 0.439937 | 0.010784 | 0.000055 | 0.418653 | 0.460909 | 37254.0 | 1.000117 |
delta_prob | 0.021924 | 0.016421 | 0.000081 | -0.009180 | 0.055029 | 38586.0 | 1.000014 |
解釈
実際にとなるサンプルの数を数えて、Aのほうが行動Xを取りやすいと言える確率を算出します。サンプリングした、を用いて、が成り立つ場合に1を、成り立たない場合に0をとる変数を定義し平均を取ります。
(trace['p_A'] - trace['p_B'] > 0).mean() 0.90927500000000006
- グループAのほうがグループBに比べて行動Xを取りやすいと言える確率は約91%で、高い確率で修正の効果があったと言えそうです。
- との平均値は約2%ほど差があるので、修正の結果行動Xをとるユーザーを2%ほど増やせた可能性が高いです。
仮説2: パターンAはパターンBよりコンバージョン率が高い
仮説1ではコンバージョンの前提となる行動Xをとるユーザーの割合について扱いましたが、コンバージョン率についても同様の問題設定で対応することができます。
モデル
with pm.Model() as model: # 事前分布として一様分布を採用 p_A = pm.Uniform('p_A', 0, 1.0) p_B= pm.Uniform('p_B', 0, 1.0) # 各行動はベルヌーイ分布から生成される obs_A = pm.Bernoulli("obs_A", p_A, observed=sample_A) obs_B = pm.Bernoulli("obs_B", p_B, observed=sample_B) # 新たな生成量として行動を起こす確率の差を定義 delta_prob = pm.Deterministic('delta_prob', p_A-p_B)
推論
with model: start = pm.find_MAP(fmin=optimize.fmin_powell) step = pm.Slice() trace = pm.sample(20000, step=step, start=start) _ = pm.traceplot(trace) _ = pm.plot_posterior(trace)
pm.summary(trace)
mean | sd | mc_error | hpd_2.5 | hpd_97.5 | n_eff | Rhat | |
---|---|---|---|---|---|---|---|
p_A | 0.114902 | 0.004827 | 0.000025 | 0.105367 | 0.124348 | 39673.0 | 0.999975 |
p_B | 0.106540 | 0.007789 | 0.000036 | 0.091163 | 0.121687 | 38234.0 | 0.999989 |
delta_prob | 0.008362 | 0.009185 | 0.000043 | -0.009795 | 0.026113 | 38091.0 | 0.999987 |
(trace['p_A'] - trace['p_B'] > 0).mean() 0.82094999999999996
解釈
パターンAがパターンBより優れている確率の推定値は約82%となりました。パターンAが優れていると完全に言い切るのは難しい結果になりました。一方で、パターンAはパターンBと比べて極端に劣る可能性は低いと考えられます。仮説1での主張と合わせることで、今後はパターンAを採用するという判断ができるかもしれません。
別解
サンプル数を、コンバージョン数をとすると、Xは二項分布に従います。事前分布に二項分布の共役事前分布であるベータ分布を採用すると、事後分布が解析的に求まるため、MCMCを使うことなくサンプルすることができます。事前分布に、のベータ分布(一様分布)を用いると、事後分布は以下の式に従います。
from scipy.stats import beta visitors_to_A = 4297 conversion_from_A = 397 visitors_to_B = 1584 conversion_from_B = 156 alpha_prior = 1 beta_prior = 1 posterior_A = beta(alpha_prior + conversion_from_A, beta_prior + visitors_to_A - conversion_from_A) posterior_B = beta(alpha_prior + conversion_from_B, beta_prior + visitors_to_B - conversion_from_B) samples = 20000 samples_posterior_A = posterior_A.rvs(samples) samples_posterior_B = posterior_B.rvs(samples) prob = (samples_posterior_A < samples_posterior_B).mean()
仮説3: パターンAはパターンBよりコンバージョンに至るまでの時間が短い
長期間コンバージョンに至らないユーザーをコンバージョンするまで追い続けるのは現実的でないです。そのため、データを用意する際に何らかの工夫が必要になります。ここでは、M日以内にコンバージョンしない場合はデータセットから除外するという処理を行いました。この制限を設けたことにより、データの分布はM日目より先が打ち切られた若干歪な分布を持つことになります。
今回はM=9としてデータセットを作成しました。頻度の低い現象が発生するまでにかかる時間間隔の分布を取り扱っているため、指数分布に従っていると考えられます。まずは分布を確認します。
コンバージョンに至るまでの時間の分布(グループA)
コンバージョンに至るまでの時間の分布(グループB)
sample_A_size: 146 sample_B_size: 391 mean_A: 1.87 mean_B: 1.78
まずはこのデータが指数分布から生成されたと仮定して、ベイズ推論を行います。指数分布のパラメータの事前分布は一様分布を採用しています。
モデル
with pm.Model() as model: # Priors for unknown model parameters lambda_A = pm.Uniform('lambda_A', 0, 5.0) lambda_B= pm.Uniform('lambda_B', 0, 5.0) obs_A = pm.Exponential("obs_A", lambda_A, observed=sample_A) obs_B = pm.Exponential("obs_B", lambda_B, observed=sample_B) ave_A = pm.Deterministic('ave_A', 1/lambda_A) ave_B = pm.Deterministic('ave_B', 1/lambda_B) delta = pm.Deterministic('delta', ave_A - ave_B)
推論
with model: start = pm.find_MAP(fmin=optimize.fmin_powell) step = pm.Slice() trace = pm.sample(20000, step=step, start=start) _ = pm.traceplot(trace)
結果
mean | sd | mc_error | hpd_2.5 | hpd_97.5 | n_eff | Rhat | |
---|---|---|---|---|---|---|---|
lambda_A | 0.536024 | 0.043165 | 0.000222 | 0.450473 | 0.619218 | 40000.0 | 0.999990 |
lambda_B | 0.583912 | 0.041367 | 0.000202 | 0.504196 | 0.666190 | 38614.0 | 0.999978 |
ave_A | 1.877783 | 0.152392 | 0.000807 | 1.589442 | 2.181460 | 39591.0 | 0.999995 |
ave_B | 1.721235 | 0.122699 | 0.000598 | 1.483570 | 1.961930 | 38377.0 | 0.999986 |
delta | 0.156548 | 0.196036 | 0.001019 | -0.213939 | 0.555821 | 39552.0 | 1.000006 |
ave_Aとave_Bのmeanと実際のデータの平均が一致していており一見良い推定結果が得られたように見えますが、モデルにデータの打ち切りの過程を含められていないため、推定されるパラメータにバイアスがかかっていると考えられます。
そこで推定する分布関数そのものをデータに合わせて定義し直してみます。具体的には指数分布の9日以降の部分を切ったTruncatedExponentialを定義します。
TruncatedExponential class
from pymc3.distributions.dist_math import bound from pymc3.distributions import draw_values, generate_samples import theano.tensor as tt import scipy.stats.distributions class TruncatedExponential(pm.Continuous): def __init__(self, lambda_, *args, **kwargs): super().__init__(*args, **kwargs) self.mode = tt.minimum(tt.floor(lambda_).astype('float32'), 1) self.lambda_ = lambda_ = tt.as_tensor_variable(lambda_) def te_cdf(self, lambda_, size=None): lambda_ = np.asarray(lambda_) dist = scipy.stats.distributions.expon() lower_cdf = dist.cdf(0) upper_cdf = 1 nrm = upper_cdf - lower_cdf sample = np.random.random(size=size) * nrm + lower_cdf return dist.ppf(sample) def random(self, point=None, size=None): lambda_ = draw_values([self.lambda_], point=point) return generate_samples(self.te_cdf, lambda_, dist_shape=self.shape, size=size) def logp(self, value): lambda_ = self.lambda_ p = pm.math.log(lambda_) - lambda_ * value - pm.math.log(1 - pm.math.exp(-9*lambda_)) log_prob = bound( p, lambda_ >= 0, value >= 0) # Return zero when mu and value are both zero return tt.switch(1 * tt.eq(lambda_, 0) * tt.eq(value, 0), 0, log_prob)
推論
with pm.Model() as model: # Priors for unknown model parameters lambda_A = pm.Uniform('lambda_A', 0.0, 5.0) lambda_B= pm.Uniform('lambda_B', 0.0, 5.0) TE('obs_A', lambda_=lambda_A, observed=sample_A) TE('obs_B', lambda_=lambda_B, observed=sample_B) ave_A = pm.Deterministic('ave_A', 1/lambda_A) ave_B = pm.Deterministic('ave_B', 1/lambda_B) delta = pm.Deterministic('delta', ava_A - ave_B) start = pm.find_MAP(fmin=optimize.fmin_powell) step = pm.Slice() trace = pm.sample(20000, step=step, start=start)
結果
mean | sd | mc_error | hpd_2.5 | hpd_97.5 | n_eff | Rhat | |
---|---|---|---|---|---|---|---|
lambda_A | 0.314618 | 0.060387 | 0.000309 | 0.194953 | 0.431916 | 37317.0 | 1.000036 |
lambda_B | 0.364973 | 0.038134 | 0.000196 | 0.288430 | 0.438147 | 40000.0 | 1.000004 |
ave_A | 3.310105 | 0.724105 | 0.003919 | 2.142754 | 4.715704 | 30976.0 | 1.000038 |
ave_B | 2.770713 | 0.298075 | 0.001547 | 2.207913 | 3.351493 | 37732.0 | 1.000012 |
delta | 0.539392 | 0.782097 | 0.004153 | -0.799306 | 2.125230 | 31161.0 | 1.000069 |
指数分布を仮定した場合に比べて、推定された平均値が高くなっていることがわかります。データ取得の事情を反映したTruncatedExponentialを使って指数分布のパラメータを推定したことによる効果だと考えられます。
との分布は、右側に伸びた左右非対称なものとなっています。9日以降のデータが無いために平均値の高い側が定まりにくいという事情を反映していると考えられます。
コンバージョンまでにかかる平均時間のグループ間の差は0を中心に分布していることがわかります。この結果から、コンバージョンに至るまでの時間短縮への寄与があるとは言えないと判断できます。
まとめ
PyMC3を使って、ベイズ推論によるA/Bテストを行いました。今回のデータセットを分析した結果、修正により行動Xは増加しており、コンバージョン率は上がる可能性が高いが、コンバージョンに至るまでの時間の短縮には繋がっているとは言えない、と判断できました。
データの生成過程をより単純なものに置き換えたり、分布関数を定義し直すなどの工夫を取り入れ、各命題に直感的な解釈を与えることができたかと思います。
最後に
弊社では、画像データや行動ログなど多様なデータを扱い、ビジネスサイドの意思決定を支える分析を行います。 多様なデータの分布を捉え、収益へ貢献できるデータサイエンティストを募集しています。