ライフタイムバリューを予測する (BG/NBDモデル)

今回は、顧客のライフタイムバリューを予測する方法の1つとして、BG/NBDモデルを紹介します。

モチベーション

顧客のライフタイムバリューを予測できると、その顧客に対してどれだけ投資して良いかがわかります。ここで言う投資は、クーポンやポイントなどのインセンティブ、DMや電話などの営業努力のことです。

また、離脱 (churn) の予測にも転用できます。過去に売上に大きく貢献してきた人が、未来のライフタイムバリューはそれと比べると極めて小さい、または、ゼロとなることが予測されるケースです。
新規顧客の獲得は既存顧客の維持よりも遥かに大変なことはよく知られてます。また、過去には大きく売上に貢献していたので、引き止められれば優良顧客まで戻る見込みは高いはずです。

BG/NBDモデル

このライフタイムバリューを予測するモデルの1つに、BG/NBD (Beta-Geometric / Negative Binomial Distribution) が2005年に提案されています (PDF) 。

repository.upenn.edu

先行していたPareto/NBD (提案論文) の簡易版として提案されたモデルです。

顧客ごとの過去のアクション (ここではトランザクションという) を使う
- 期間Tにおけるトランザクション数
  - Tは最初のトランザクションの時刻 $t_0$ を起点
- 最後にトランザクションした時間 $t_x$
予測したいことは3つ
- 長さtの時間のトランザクション数の (全体の) 期待値 $E[X(t)$ ]
- 長さtの時間でトランザクション数xとなる (全体の) 確率 $P[X(t)=x$ ]
- (T, T+t] の時間における、ある顧客 (X=x, t, T) のトランザクション数の期待値 $E(\Upsilon(t) | X=x, t, T)$

ある顧客Aをイメージするとこんな感じです。期間T中にx回トランザクションが発生し、最後のトランザクションは時間 $t_x$ でした。では、Tの後の時間tの間に何回トランザクションが発生するか、を予測したいのです。

前提

BG/NBDモデルは5つの前提を置いています。

顧客はアクティブ (active) と離脱 (inactive) の2つの状態を持つ
- アクティブな顧客のみがトランザクションを行う

前提1. アクティブな顧客はポアソン分布に従ってトランザクションする

$f(t_j | t_{j-1}, \lambda) = \lambda e^{t_j - t_{j-1}}, t_{j} > t_{j-1} \geq 0$

前提2. λはガンマ分布に従う

$f(\lambda | r, \alpha) = \frac{\alpha^{r} \lambda^{r-1} e^{-\lambda\alpha}}{\Gamma(r)}, \lambda>0$

前提3. トランザクション後に確率pで離脱する

Pareto/NBDと異なる前提を置いています。Pareto/NBDでは、離脱は任意のタイミングで起きる前提に立っています。

P(inactive immediately after j-th transaction) = $p(1-p)^{j-1}, j=1, 2, ...$

前提4. pはベータ分布に従う

B(a, b)はベータ関数です。

$f(p|a, b) = \frac{p^{a-1} (1-p)^{b-1}}{B(a, b)}, 0 \leq p \leq 1$

$B(a, b) = \frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}$

前提5. λとpは顧客間で独立

そのため前提1.〜4.で挙げたパラメータがわかれば、他の顧客を考慮すること無く、λとpを計算できます。

パラメータの推定

上の前提をまとめると、求める必要があるパラメータは $r, \alpha, a, b$ の4つとなります。尤度関数を定義してパラメータを求めます。

尤度関数Lは以下の式となります。 $\delta_{x>0}$ は x>0のときに1、x=0で0となります。したがって、まだ1回もトランザクションしていない場合は、第1項だけで計算されます。

$L(r \alpha, a, b | X=x, t_x, T) = \frac{B(a, b + x)}{B(a, b)} \frac{\Gamma(r + x)a^r}{\Gamma(r)(\alpha + T)^{r + x}} + \delta_{x>0} \frac{B(a+1, b+x-1)}{B(a, b)} \frac{\Gamma(r + x)\alpha^r}{\Gamma(r)(a + t_x)^{r+x}}$

最大化する対数尤度関数LLは以下の式です。Nは顧客数です。

$LL(r, \alpha, a, b) = \sum_{i=1}^{N} ln[L(r, \alpha, a, b | X_i = x_i, t_{x_i}, T_i)$ ]

予測

上で求まったパラメータを使って、3つの値を予測します。

トランザクション数がxとなる確率 (全体)

$P(X(t)=x | r, \alpha, a, b) = \frac{B(a, b+x)}{B(a, b)} \frac{\Gamma(r+x)}{\Gamma(r)x!} (\frac{\alpha}{\alpha+t})^r (\frac{1}{\alpha+t})^x + \delta_{x>0} \frac{B(a+1, b+x-1)}{B(a, b)} (1-(\frac{\alpha}{\alpha+t})^r (\sum_{j=0}^{x-1} \frac{\Gamma(r+j)}{\Gamma(r)j!}(\frac{t}{\alpha+t})^j))$

トランザクション数の期待値 (全体)

$E(X(t) | r, \alpha, a, b) = \frac{a+b-1}{a-1}(1-(\frac{\alpha}{\alpha+t})^r {}_2F_1(r, b; a+b-1; \frac{t}{\alpha+t}))$

${}_2F_1()$ はガウス超幾何関数 (Gaussian hypergeometric function) です。このあたりは勉強不足でよくわからんですが、以下で計算できるようです。

$_2F_1(r, b; a+b-1; \frac{t}{\alpha+t}) = \sum_{n=0}^{\infty} \frac{(r)_n (b)_n}{(a+b-1)_n} \frac{(\frac{t}{\alpha+t})^n}{n!}$

$(x)_n$ はポッホハマー記号で $(x)_0 := 1$ , $(x)_n := \prod_{k=0}^{n-1}(x+k)$ で定義されます。

ある顧客のトランザクション数の期待値

$E(\Upsilon(t) | X=x, t_x, T, r, \alpha, a, b) = \frac{a+b+x-1}{a-1} \frac{1-(\frac{\alpha + T}{\alpha + T + t})^{r+x}{}_2F_1(r + x, b + x; a+b+x-1; \frac{t}{\alpha + T + t})}{1 + \delta_{x>0}\frac{a}{b+x-1}(\frac{\alpha + T}{\alpha + t_x})^{r+x}}$

複雑な式に見えますが、パラメータさえ定まってしまえば、簡単に計算できます。以下のパラメータ (原文Figure 1から抜粋) を使ってPythonで計算してみます。

パラメータ	値
r	0.243
α	4.414
a	0.793
b	2.426

import math
from scipy.special import hyp2f1  # 超幾何級数2F1の実装

# パラメータ
r, alpha, a, b = 0.243, 4.414, 0.793, 2.426

# トランザクション数の期待値
def expect_transaction_count(T, x, t_x, t):
    term1 = (a + b + x - 1)/(a - 1)
    
    numerator = 1 - ((alpha + T) / (alpha + T + t))**(r + x) * hyp2f1(r + x, b + x, a+b+x-1, t / (alpha + T + t))
    denominator = 1 + (a / (b + x - 1))*((alpha + T) / (alpha + t_x))**(r + x)
    term2 = numerator / denominator
    
    return term1 * term2

3パターンの顧客でトランザクション数の期待値を比較してみます。いずれもT=t=30、つまり過去の30期間 (例えば日数や週数) のxとt_xを使って、未来の30期間のトランザクション数を予測します。

顧客1と顧客2を比較すると、過去のトランザクション数xが多いほうが、未来のトランザクション数も増えると予想されることがわかります
顧客1と顧客3を比較すると、最後のトランザクションt_xが遠いほうが、未来のトランザクション数も減ると予想されることもわかります

# 顧客1
T, x, t_x, t = 30, 4, 25, 30
expect_transaction_count(T, x, t_x, t)
# 2.3526752183407695

# 顧客2
T, x, t_x, t = 30, 8, 25, 30
expect_transaction_count(T, x, t_x, t)
# 4.388865813295875

# 顧客3
T, x, t_x, t = 30, 4, 15, 30
expect_transaction_count(T, x, t_x, t)
# 1.1367813390968273

所感

BG/NBDモデルについての、私の所感をまとめます。

期間T内のトランザクション数と最終トランザクション日時しか見ていないため、周期性のある顧客については表現できていない
- 周期性のある顧客は、tの設定によっては、予想されるトランザクション数と実際のトランザクション数が大きく異なることが考えられる
〜のトランザクション発生日時は捨ててしまって良いのだろうか？
- 例えば過去30日間で見たとき、同じx=10だとしても、1〜9日目と30日目にトランザクションした人と、平均的にトランザクションした人は同じと扱って良いかは疑問が残る (後者の方が継続して利用しそう)
あくまでトランザクション数であるため、ビジネスサイドが求めるライフタイムバリューと揃えるためには何らかの計算が必要
- その顧客の過去の平均金額と掛け算する、など

まとめ

顧客のライフタイムバリューを表現するBG/NBDモデルを紹介しました。