🦄

ビジネス判断のためのベイジアンABテスト—KPIが連続変数のとき

2024/03/25に公開

Python

分析環境

本記事では、下記の環境で分析を実行していきます。

	バージョン
macOS	14.1
xcode-select C++コンパイラ	version 2403
Python	3.11.4
CmdStan	2.34.1
cmdstanpy	1.2.1

CmdStanやcmdstanpyの導入方法は、cmdstanpyの公式ドキュメントのInstallationページを参照してください。また、CmdStanを利用する場合、Stan言語をコンパイルするためのC++コンパイラが必要になる点にご留意ください。

問題設定

あなたはネット通販事業を行う企業でデータ分析業務を担当しています。売上を伸ばすための施策で、クーポンを配布する群と配布しないコントロール群で100人ずつに無作為に割り当ててABテストを実施しています。顧客1人あたりのクーポンのコストは850円です。取り扱っている商品・サービスの利益率は0.5と仮定します。このクーポンの効果について下記のようなビジネス上の疑問に回答する要望があります。どのような分析を行うべきでしょうか？

顧客1人あたり平均売上金額に対するクーポンの効果が0円以上である確率はどの程度でしょうか？
顧客1人あたり平均売上金額に対するクーポンの効果が、コストの850円以上である確率はどの程度でしょうか？
ROIはどの程度でしょうか？

データ

事前に今後利用するライブラリを読み込みます。

import numpy as np
import pandas as pd
import plotly.express as px
import cmdstanpy

下記のコードで架空のデータを生成します。

generate_data.py

# クーポン群とコントロール群の顧客数
N_coupon = 100
N_control = 100

# クーポン群とコントロール群の顧客1人あたりの平均購入金額
mu_coupon = 5500
mu_control = 4500

# クーポン群とコントロール群の顧客の購入金額の標準偏差
sigma_coupon = 750
sigma_control = 550

# 乱数の種を指定
np.random.seed(123)

# クーポン群とコントロール群の顧客の購入金額
# ここではそれぞれ異なるパラメータを持つ正規分布から生成されていることを仮定
y_coupon = np.random.normal(mu_coupon, sigma_coupon, N_coupon).round()
y_control = np.random.normal(mu_control, sigma_control, N_control).round()

ここで生成したデータをデータフレームに格納して、ヒストグラムを描画して分布の確認を行います。

visualize_data.py

df = pd.DataFrame({"Coupon": y_coupon, "Control": y_control}).melt(
    var_name="group", value_name="y"
)

# ヒストグラムを描画
px.histogram(
    df,
    x="y",
    color="group",
    marginal="box",
    opacity=0.8,
    width=600,
    height=400,
    barmode="overlay",
    labels={"y": "購入金額", "group": "グループ"},
)

上記コードで描画したヒストグラムは以下のようなものになります。分布の重なりはありつつも、設定したとおりクーポン群の方が購入金額が平均的に大きい傾向が読み取れます。

データのヒストグラム
クーポン群とコントロール群の顧客1人あたりの購入金額の分布

ベイジアンABテストによる分析

ベイジアンABテストとは

ベイジアンABテストは、ベイズ統計学に基づいた群間比較を行うための分析方法です。ビジネス場面ではCVRの分析に頻繁に利用されています。ベイズ統計学では、下記のベイズの定理を理論的な基礎にしています。

\begin{equation} P(A|B) = \frac{P(B|A)P(A)} {P(B)} \end{equation}

このベイズの定理の $A$ を推測したいパラメータ（ $\theta$ ）に、 $B$ を観測されたデータ（ $x$ ）に置き換えて考えると、

\begin{equation} P(\theta|x) = \frac{P(x|\theta)P(\theta)} {P(x)} \end{equation}

のようになります。
この式の意味を下記の表にまとめました。

要素	意味
$P(\theta\|x)$	事後確率（分布）。今回のデータが観測されたときに、パラメータが特定の値である確率（分布）。
$P(x\|\theta)$	尤度（ゆうど）。パラメータが特定の値であるときの、今回のデータの観測されやすさ。
$P(\theta)$	事前確率（分布）。今回のデータが観測される前の状態において、パラメータが特定の値である確率。
$P(x)$	正規化定数。事後確率（分布）が確率分布となるようにするための定数。

このように事後分布を求めることで、データが観測された後のパラメータとその不確実性について議論することが可能になります。

Stan言語について

本記事では、ベイジアンABテストの実行にあたってStan言語を利用します。Stan言語は、事後分布からの乱数を生成するための確率的プログラミング言語です。

Stan言語では、主に下記の4つのブロックに推定したモデルに関する記述を行います。

ブロック	記述する内容
data{}	観測データの変数宣言を行う。たとえば、CV数や、売上金額データ、観測データ数など。
parameters{}	観測データから推定するパラメータの変数宣言を行う。たとえば、CVRや、売上の母平均など。
model{}	観測データと推定するパラメータの関係を記述する（尤度関数）。また、事前分布の記述も行う。たとえば、施策A群のCV数（data{}で宣言）は、施策A群のCVRをパラメータ(parameters{}で宣言)にもつ二項分布から生成されるなど。
generated quantities{}	推定するパラメータを利用して、計算したい量に関する変数宣言と定義を行う。たとえば、CVRの差や比、推定した施策Aの売上の母平均が施策Bの売上の母平均よりも大きい確率など。

Stan言語およびベイズ統計学の詳細は、たとえば下記の資料を参照してください。

分析の実行

今回は顧客の購入金額が正規分布に従っていることを仮定して分析を行いたいと思います。頻度論的仮説検定におけるウェルチのt検定と同等の両群の分散が異なるというモデルを仮定します。

本記事ではモデルのわかりやすさを考慮して、購入金額を生成する正規分布のパラメータの事前分布に範囲が十分に広い一様分布を利用します。より汎用的な事前分布の設定である効果量の事前分布に標準コーシー分布を置くモデルについては下記の書籍の第八章「ガウス分布の平均の比較」などを参照してください。

また、本記事ではビジネス場面のベイジアンABテストで主流のパラメータ推定に基づいた分析を行いたいと思います。パラメータ推定に基づいた分析とモデル比較に基づいた分析の比較は下記の資料に記載があります。

モデル比較に基づいた分析を行いたい場合は、たとえば下記のような資料を参照してください。

本記事で利用するモデル（尤度と事前分布）を数式で表現すると下記になります。

\begin{aligned} \mu_{coupon} &\sim Uniform(0, 100000) \\ \mu_{control} &\sim Uniform(0, 100000) \\ \sigma_{coupon} &\sim Uniform(0, 100000) \\ \sigma_{control} &\sim Uniform(0, 100000) \\ y_{coupon} &\sim Normal(\mu_{coupon}, \sigma_{coupon}) \\ y_{control} &\sim Normal(\mu_{control}, \sigma_{control}) \end{aligned}

また、問題設定の項に記載したビジネス上の疑問に回答するために、下記のようにパラメータを利用した生成量を計算します。

\begin{aligned} \delta &= \mu_{coupon} - \mu_{control} \\ \delta_{cost} &= \mu_{coupon} - \mu_{control} - Cost_{coupon} \\ \pi_{coupon} &= \mu_{coupon} \cdot r - Cost_{coupon} \\ roi &= \pi_{coupon} / Cost_{coupon} \\ \end{aligned}

ここで、 $\pi_{coupon}$ はクーポン群の利益金額で、 $Cost_{coupon}$ はクーポンの費用を表しています。

ここで計算した $\delta$ および $\delta_{cost}$ を利用して

$\delta$ が0以上の確率
$\delta_{cost}$ が0以上の確率

についても生成量で計算します。

ここまでに整理したモデルをStan言語で表現すると下記になります。

coupon_ab_testing.stan

data {
  int N_coupon; // クーポン群の数
  int N_control; // コントロール群の数
  vector[N_coupon] y_coupon; // クーポン群のデータ
  vector[N_control] y_control; // コントロール群のデータ
  real<lower=0, upper=1> profit_rate; // 利益率
  real cost_coupon; // クーポンのコスト
}

parameters {
  real mu_coupon; // クーポン群の平均
  real mu_control; // コントロール群の平均
  real<lower=0> sigma_coupon; // クーポン群の標準偏差。0以上
  real<lower=0> sigma_control; // コントロール群の標準偏差。0以上
}

model {
  // 事前分布
  mu_coupon ~ uniform(0, 100000); // 0以上100000以下の一様分布
  mu_control ~ uniform(0, 100000); // 0以上100000以下の一様分布
  sigma_coupon ~ uniform(0, 100000); // 0以上100000以下の一様分布
  sigma_control ~ uniform(0, 100000); // 0以上100000以下の一様分布

  // 尤度
  y_coupon ~ normal(mu_coupon, sigma_coupon);
  y_control ~ normal(mu_control, sigma_control);
}

generated quantities {
  real delta; // クーポン群とコントロール群の平均の差
  real delta_cost; // クーポン群とコントロール群の平均の差からコストを引いた値
  real<lower=0, upper=1> p_delta_over_zero; // クーポン群とコントロール群の平均の差が0以上の確率
  real<lower=0, upper=1> p_delta_over_cost; // クーポン群とコントロール群の平均の差がコスト以上の確率
  real<lower=0> profit_coupon; // クーポン群の利益
  real<lower=0> roi; // ROI

  delta = mu_coupon - mu_control; // クーポン群とコントロール群の平均の差
  delta_cost = delta - cost_coupon; // クーポン群とコントロール群の平均の差からコストを引いた値
  p_delta_over_zero = step(delta); // クーポン群とコントロール群の平均の差が0以上の確率
  p_delta_over_cost = step(delta_cost); // クーポン群とコントロール群の平均の差がコスト以上の確率
  profit_coupon = (mu_coupon * profit_rate) - cost_coupon; // クーポン群の利益
  roi = profit_coupon / cost_coupon; // ROI
}

このモデルをコンパイルします。

compile_model.py

# モデルのコンパイル
model = cmdstanpy.CmdStanModel(stan_file="./coupon_ab_testing.stan")

モデルをコンパイルしたあとは、下記のコードでMCMC法による事後分布からのサンプリングを行います。

mcmc_sampling.py

# データの設定
data = {
    "N_coupon": N_coupon,
    "N_control": N_control,
    "y_coupon": y_coupon,
    "y_control": y_control,
    "profit_rate": 0.5,
    "cost_coupon": 850,
}

# サンプリング
fit = model.sample(data=data, seed=123, chains=4, iter_warmup=2000, iter_sampling=2000)

事後分布からのサンプリングが終了したら、事後分布の収束を確認します。

check_convergence.py

# 収束の確認
print(fit.diagnose())

ここでは問題なく収束しているという記述が出力されていると思います。

収束の診断

また、トレースプロットを確認しても、問題なく事後分布からのサンプリングが行えていることが確認できます。

traceplot.py

# 事後分布からのサンプルを抽出
samples = fit.draws_pd(
    vars=[
        "chain__",
        "iter__",
        "mu_coupon",
        "mu_control",
        "sigma_coupon",
        "sigma_control",
    ]
)

samples = samples.melt(id_vars=["chain__", "iter__"], var_name="parameter")

px.line(
    samples,
    x="iter__",
    y="value",
    color="chain__",
    facet_col="parameter",
    facet_col_wrap=2,
    width=600,
    height=400,
    labels={"iter__": "反復回数", "value": "パラメータ"},
)

事後分布のトレースプロット
推定したパラメータごとのトレースプロット：うまく推定できていることが読み取れる

MCMCの収束については、下記の資料を参照してください。

事後分布の収束を確認できたので、結果の確認を行いたいと思います。

# 結果のサマリーを表示
fit.summary().round(3)

結果のサマリー

結果の表から、 $\mu_{coupon}, \mu_{control}, \sigma_{coupon}, \sigma_{control}$ それぞれの設定通りに上手く推定できていることが読み取れます。

ここで、問題設定に記載したビジネス上の疑問に回答していきたいと思います。

顧客1人あたり平均売上金額に対するクーポンの効果が0円以上である確率はどの程度でしょうか？

表の p_delta_over_zero の行から100%の確率であることが読み取れます。強い確信を持って、クーポンを配布された顧客の平均購入金額は、配布されていない顧客の平均購入金額よりも大きいことが言えます。

顧客1人あたり平均売上金額に対するクーポンの効果が、コストの850円以上である確率はどの程度でしょうか？

表の p_delta_over_cost の行から、およそ96.4%の確率でコストを上回る効果があるということができそうです。また、表の delta_cost はおよそ181.744円程度であることがわかります。この確率や金額が大きいのか小さいのかについては、具体的なビジネス状況に依存する問題だと思います。

ROIはどの程度でしょうか？

表の roi の行から、ROIはおよそ2.247倍であることがわかります。こちらの値が大きいのか小さいのかについても、具体的なビジネス状況に依存する問題だと思います。

結論

ベイジアンABテストを活用してデータの分析を行うことで、天気予報の降水確率と同じような粒度の情報になっていて、ビジネス判断に活用しやすい材料を用意することが可能です。
また、たとえばROIがある特定の基準値よりも大きい確率など自由に生成量を定義することで、さまざまなビジネス上の判断に活用できる量を計算することができます。

要素	意味
$P(\theta\|x)$	事後確率（分布）。今回のデータが観測されたときに、パラメータが特定の値である確率（分布）。
$P(x\|\theta)$	尤度（ゆうど）。パラメータが特定の値であるときの、今回のデータの観測されやすさ。
$P(\theta)$	事前確率（分布）。今回のデータが観測される前の状態において、パラメータが特定の値である確率。
$P(x)$	正規化定数。事後確率（分布）が確率分布となるようにするための定数。