f:id:vasilyjp:20180927112657j:plain

こんにちは。データチームの後藤です。

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と数値化することにより、これらのサンプルが確率 $p$ のベルヌーイ分布から生成されると考えます。また $p$ に関しての情報はないので、事前分布として一様分布を採用します。

データの集計結果は以下のようになりました。

- グループA
    - サンプル数: 1602
    - 行動Xを取った数: 740
    - 平均: 0.462
- グループB
    - サンプル数: 2121
    - 行動Xを取った数: 933
    - 平均: 0.440

グループAのほうが行動Xを取った割合が高いので、修正の効果が出ているかもしれません。

PyMC3のフレームワークを使えば、今回のモデルを以下のように記述できます。行動Xを取る確率が増えたかどうかも知りたいので、 $\delta_{prob}$ という $p_B$ と $p_A$ の差を新たな生成量として定義します。

モデル

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)

$p_A$ 、 $p_B$ とその差 $\delta_{prob}$ の標本は均衡分布に収束していることが確認できます。各変数の分布も左右対称のベル型をとっていることもわかります。 f:id:vasilyjp:20180219092705p:plain

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

解釈

実際に $p_A > p_B$ となるサンプルの数を数えて、Aのほうが行動Xを取りやすいと言える確率を算出します。サンプリングした $p_A$ 、 $p_B$ を用いて、 $p_A > p_B$ が成り立つ場合に1を、成り立たない場合に0をとる変数 $x_i$ を定義し平均を取ります。

$E[p_A > p_B ] \approx \sum_{i=0}^N x_i /N$

(trace['p_A'] - trace['p_B'] > 0).mean()

0.90927500000000006

グループAのほうがグループBに比べて行動Xを取りやすいと言える確率は約91％で、高い確率で修正の効果があったと言えそうです。
$p_A$ と $p_B$ の平均値は約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)

f:id:vasilyjp:20180219082016p:plain

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を採用するという判断ができるかもしれません。

別解

サンプル数を $N$ 、コンバージョン数を $X$ とすると、Xは二項分布に従います。事前分布に二項分布の共役事前分布であるベータ分布を採用すると、事後分布が解析的に求まるため、MCMCを使うことなくサンプルすることができます。事前分布に $\alpha = 1$ 、 $\beta=1$ のベータ分布（一様分布）を用いると、事後分布は以下の式に従います。

$Beta(\alpha + X, \beta + N - X)$

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）

f:id:vasilyjp:20180219085620p:plain

コンバージョンに至るまでの時間の分布（グループB）

f:id:vasilyjp:20180219085626p:plain

sample_A_size: 146
sample_B_size: 391
mean_A: 1.87
mean_B: 1.78

まずはこのデータが指数分布から生成されたと仮定して、ベイズ推論を行います。指数分布のパラメータ $\lambda$ の事前分布は一様分布を採用しています。

モデル

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)

f:id:vasilyjp:20180219090648p:plain

結果

	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)

f:id:vasilyjp:20180219093236p:plain

結果

	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を使って指数分布のパラメータ $\lambda$ を推定したことによる効果だと考えられます。

$ave_A$ と $ave_B$ の分布は、右側に伸びた左右非対称なものとなっています。9日以降のデータが無いために平均値の高い側が定まりにくいという事情を反映していると考えられます。
コンバージョンまでにかかる平均時間のグループ間の差は0を中心に分布していることがわかります。この結果から、コンバージョンに至るまでの時間短縮への寄与があるとは言えないと判断できます。