📈

PyMC Marketingで学ぶBG/NBD入門

に公開
Python
marketing
bayesian
pymc
clv
tech

はじめに

PyMC MarketingのBG/NBDの実装をチュートリアルに沿って試してみました。
自分の解釈を加えている部分もあるため、間違いなどありましたらコメントいただけますと幸いです。

BG/NBDモデルとは

BG/NBD（Beta-Geometric/Negative Binomial Distribution）モデルは、非契約型ビジネスにおける顧客の購買行動を予測するための確率的モデル。

契約型ビジネスと非契約型ビジネス

特徴 例
契約型 契約, 更新, 解約のプロセスがあり、顧客がアクティブかどうかが明確 保険、動画配信サービスの有料会員
非契約型 取引ごとに完結しているビジネス。顧客がアクティブかどうかがわからない 小売、飲食

非契約型ビジネスでは、ふつう顧客が「離脱したのか、たまたま買っていないだけか」を観測から直接判定できません。BG/NBDは、この不確実性を明示的にモデル化し、将来の購入回数とアクティブ確率（生存確率） を推定することで、意思決定に必要な量を提供することができます。
これらの推定値は、CLV（Customer Lifetime Value）の算出やリテンション/再アクティベーション施策の設計に直結するため、マーケティング分析で役に立ちます。

BG/NBDモデルでは以下の5つの仮定を置いて、顧客の購買行動を確率的モデルとして説明します。

  1. 顧客が"アクティブ"な場合、購入から次の購入までの期間は、購買パラメータ \lambda の指数分布に従う。
f(t_j | t_{j-1}; \lambda) = \lambda e^{-\lambda(t_j - t_{j-1})}, \quad t_j\ge t_{j-1} \ge 0
  1. 顧客ごとの \lambda の異質性はガンマ分布に従う。
f(\lambda \mid r, \alpha) = \frac{\alpha^{\,r}}{\Gamma(r)}\, \lambda^{r-1} e^{-\alpha \lambda}, \quad \lambda > 0
  1. 毎回の購入の直後に、顧客は確率 p で"離脱"する
  2. 顧客ごとの p の異質性はベータ分布に従う。
f(p|a, b) = \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}p^{a-1}(1-p)^{b-1}, \quad 0\le p \le 1
  1. \lambda と p は顧客ごとに独立である。

チュートリアル

必要なライブラリのimport

import arviz as az
import matplotlib.pyplot as plt
import pandas as pd
import xarray as xr
from fastprogress.fastprogress import progress_bar

from pymc_marketing import clv

サンプルデータを読み込む

data_path = "https://raw.githubusercontent.com/pymc-labs/pymc-marketing/main/data/clv_quickstart.csv"

df = pd.read_csv(data_path)

df.head()
index frequency recency T monetary_value
0 2 30.43 38.86 22.35
1 1 1.71 38.86 11.77
2 0 0.00 38.86 0.00
3 0 0.00 38.86 0.00
4 0 0.00 38.86 0.00

このデータセットはCDNOWという音楽通販サイトの実購買履歴（1997年公開）に基づいており、CLV分析の代表的な教材として利用されます。
一レコードが一人の顧客に対応しており、顧客の購買行動がfrequency, recency, T, monetary_valueという値に要約されています。
それぞれの意味は...

  • frequency: 購買回数 -1 (-1は初回の購買を除外するためで、リピート回数と解釈できる)
  • recency: 最初に購買した日から最後に購買した日までの期間
  • T: 最初に購買した日から参照時点(ex.モデル作成日)までの期間
  • monetary_value: 購買金額の合計

RFM分析の文脈で「Recency」は、「分析時点から直近の購入日までの期間」として定義されることもあり、定義が異なる点に注意が必要です。
同様の値は、T - recency で表現されます。

BG/NBDモデルでは、パラメータの推定のためにfrequency, recency, T の3つを利用します。

# pymc-marketingの実装に合わせるため顧客ID列を"customer_id"という名前で作成します

data = (
    df.reset_index()
    .rename(columns={"index": "customer_id"})
    .drop(columns="monetary_value")
)

パラメータ推定

# モデルのインスタンス化
model = clv.BetaGeoModel(data=data)

# build_model を行うことでモデルの設定を確認できる(省略化)
model.build_model()
model

BG/NBD
            alpha ~ Weibull(2, 10)
      phi_dropout ~ Uniform(0, 1)
    kappa_dropout ~ Pareto(1, 1)
                r ~ Weibull(2, 1)
                a ~ Deterministic(f(kappa_dropout, phi_dropout))
                b ~ Deterministic(f(kappa_dropout, phi_dropout))
recency_frequency ~ BetaGeoNBD(a, b, r, alpha, <constant>)

model.graphviz()

graph

上の図はBG/NBDモデルの構造を表しています。
最初の仮定でおいた α, r (ガンマ分布), a, b (ベータ分布)の他に、a, bのハイパーパラメータとしてphi_dropout, kappa_dropoutが追加されていることがわかります。
これらは、パフォーマンス向上のためにデフォルトで設定されているものだそうです
以下のようにすると、a, bに直接事前分布を設定することもできます。

# 直接事前分布を設定
model_config = {
    "a": Prior("Gamma", alpha=2, beta=1),  # aの事前分布
    "b": Prior("Gamma", alpha=2, beta=1),  # bの事前分布
}

モデルの当てはめ

NUTSサンプラーによってパラメータの事後分布をサンプリングします。

sample_kwargs = {
    "draws": 2_000,
    "chains": 4,
    "target_accept": 0.9,
    "random_seed": 42,
}

idata_mcmc = model.fit(**sample_kwargs)

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 2 jobs)
NUTS: [alpha, phi_dropout, kappa_dropout, r]


Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 10 seconds.

サンプリング結果の確認

model.fit_summary()
param mean sd hdi_3% hdi_97% r_hat
alpha 4.503 0.394 3.780 5.257 1.00
phi_dropout 0.249 0.020 0.211 0.286 1.00
kappa_dropout 3.188 0.942 1.721 4.935 1.00
r 0.245 0.013 0.220 0.270 1.00
a 0.782 0.195 0.462 1.141 1.00
b 2.407 0.756 1.169 3.733 1.00

MCMCの収束の基準として、r_hat が 1 に近いかどうかで判断します。
今回の結果は問題なさそうです。

トレースプロットからも収束しているかどうかをチェックします。

axes = az.plot_trace(
    data=model.idata,
    compact=True,
    kind="rank_bars",
    backend_kwargs={"figsize": (12, 9), "layout": "constrained"},
)
plt.gcf().suptitle("BG/NBD Model Trace", fontsize=18, fontweight="bold");

trace

左側のプロットで、チェインごとの事後分布がほぼ重なっているので、収束しているとして良さそうです

モデルの当てはまりチェック

以下のコードでは、
パラメータの事後分布を元に、学習データと同一条件のもとで、観測されたコホートの frequency 分布を事後予測で再生成し、実際の学習データのfrequencyと比較を行っています。

clv.plot_expected_purchases_ppc(model, ppc="posterior")

Sampling: [recency_frequency]

Output()

ppc

横軸がfrequency(リピート回数)、縦軸が人数です
モデルが予測するfrequencyの分布は、データとほぼ一致しており、顧客の購買行動をおおむね捉えられています。

モデルの活用

一部の顧客に対して、将来の購買回数やアクティブである確率を予測していきます

# 顧客をサンプリング
example_customer_ids = [1, 6, 10, 18, 45, 1412]

data_small = data.query("customer_id.isin(@example_customer_ids)")

data_small.head(6)
customer_id frequency recency T
1 1 1.71 38.86
6 1 5.00 38.86
10 5 24.43 38.86
18 3 28.29 38.71
45 12 34.43 38.57
1412 14 30.29 31.57

上記の6人の顧客についてみていきます。
データから読み取れることは、

  • 1, 6 はfrequencyが1, recencyが小さく、Tが大きい
    • 初期に1回だけリピートしたが、最後の購入から期間が経っている → 離脱している可能性が高い
  • 45, 1412 はfrequencyが大きく, recencyとTが大きく、両者の差が小さい
    • 観測期間の初期から頻繁に購入している → アクティブである可能性が高い

将来の購買回数予測

t時点後の購入回数を予測します

steps = 90

expected_num_purchases_steps = xr.concat(
    objs=[
        model.expected_purchases(
            data=data_small,
            future_t=t,
        )
        for t in progress_bar(range(steps))
    ],
    dim="t",
).transpose(..., "t")

fig, axes = plt.subplots(
    nrows=len(example_customer_ids)//2,
    ncols=2,
    figsize=(12, 8),
    sharex=True,
    sharey=True,
    layout="constrained",
)

for i, customer_id in enumerate(example_customer_ids):
    ax = axes[i//2][i%2]
    customer_expected_num_purchases_steps = expected_num_purchases_steps.sel(
        customer_id=customer_id
    )
    az.plot_hdi(
        range(steps),
        customer_expected_num_purchases_steps,
        hdi_prob=0.94,
        color="C0",
        fill_kwargs={"alpha": 0.3, "label": "$94 \\%$ HDI"},
        ax=ax,
    )
    az.plot_hdi(
        range(steps),
        customer_expected_num_purchases_steps,
        hdi_prob=0.5,
        color="C0",
        fill_kwargs={"alpha": 0.5, "label": "$50 \\%$ HDI"},
        ax=ax,
    )
    ax.plot(
        range(steps),
        customer_expected_num_purchases_steps.mean(dim=("chain", "draw")),
        color="C0",
        label="posterior mean",
    )
    ax.legend(loc="upper left")
    ax.set(title=f"Customer {customer_id}", xlabel="t", ylabel="purchases")

axes[-1][0].set(xlabel="steps")
axes[-1][1].set(xlabel="steps")
plt.gcf().suptitle("Expected Number of Purchases", fontsize=18, fontweight="bold")

expected-purchases

図は横軸: t時点後の、縦軸: 購入回数の期待値を表しています。
頻繁に購入していた顧客45, 1412は、将来にも購入回数が増えていくことが予測されます
一方で、一見客だった1, 6 は将来の購入回数もほぼゼロです

顧客のアクティブである確率

次にt時点後に顧客がアクティブ(離脱していない)かどうかの確率を出してみます

steps = 90

future_alive_all = []

for t in progress_bar(range(steps)):
    future_data = data_small.copy()
    future_data["T"] = future_data["T"] + t
    future_alive = model.expected_probability_alive(data=future_data)
    future_alive_all.append(future_alive)

expected_probability_alive_steps = xr.concat(
    objs=future_alive_all,
    dim="t",
).transpose(..., "t")

fig, axes = plt.subplots(
    nrows=len(example_customer_ids)//2,
    ncols=2,
    figsize=(12, 8),
    sharex=True,
    sharey=True,
    layout="constrained",
)

# axes = axes.flatten()

for i, customer_id in enumerate(example_customer_ids):
    ax = axes[i//2][i%2]
    customer_expected_probability_alive_steps = expected_probability_alive_steps.sel(
        customer_id=customer_id
    )
    az.plot_hdi(
        range(steps),
        customer_expected_probability_alive_steps,
        hdi_prob=0.94,
        color="C1",
        fill_kwargs={"alpha": 0.3, "label": "$94 \\%$ HDI"},
        ax=ax,
    )
    az.plot_hdi(
        range(steps),
        customer_expected_probability_alive_steps,
        hdi_prob=0.5,
        color="C1",
        fill_kwargs={"alpha": 0.5, "label": "$50 \\%$ HDI"},
        ax=ax,
    )
    ax.plot(
        range(steps),
        customer_expected_probability_alive_steps.mean(dim=("chain", "draw")),
        color="C1",
        label="posterior mean",
    )
    ax.legend(loc="upper right")
    ax.set(title=f"Customer {customer_id}", ylabel="probability alive", ylim=(0, 1))

axes[-1][0].set(xlabel="steps")
axes[-1][1].set(xlabel="steps")
plt.gcf().suptitle(
    "Expected Probability Alive over Time", fontsize=18, fontweight="bold"
);

prob-alive

図は横軸: t時点後の、縦軸: アクティブ確率を表しています。
この図は t時点後まで、購入が起きなかったときのアクティブ確率を予測している点に注意する必要があります。

そのため、元々頻繁に購入していた45や1412の購入が途絶えると、急激にアクティブ確率が減少するという予測になります。

おわりに

最後まで読んでいただきありがとうございました。
今回はドキュメントのノートブックに沿ってBG/NBDモデルを実装してみました。
統計モデルやMCMCについて、まだまだ理解できていないことが多いですが、少しずつ勉強していきたいと思います。

参考

Discussion