終生顧客購買行動予測模型

• 
• 
• 
• 

はじめに

こんにちは。GMOインターネット 嘉山です。

新卒1年目で参画したデータ分析系プロジェクトにて、「Buy Till You Die（BTYD）モデル」を調査する機会があり、今回はそこで学んだ知見を紹介させていただきます。

BTYDモデルとは、手元の顧客購買履歴データから「将来を含めた顧客の合計購買金額」を推定するための統計モデリング手法の一つです。このモデルは「ベイズ統計学」をベースとした構造を有し、ただデータが持つ特性を確率分布で表現するだけでなく、モデル自体が内包する「不確実性」も捉えることができます。

余談ですが、この “Buy Till You Die” というモデル名、冷静に考えるとだいぶ攻めてますよね。。。直訳すると、”死ぬまで買え”。普通にホラーです。

ということで私もこのモデル名に負けないよう、人間の可読性を軽視した呪文のような数式をただひたすら書き起こし、名実ともに狂気じみた記事をお届けできるよう頑張りたいと思います。

次章より早速意味不明な数式が連発し、じわじわと読者に不快感と苦痛を与えていく構成となっております。早々に離脱したくなる方が続出する未来が予想されますが、新卒2年目の若気の至りに最後までお付き合いいただけますと幸いです。

Buy Till You Die（BTYD）モデルとは

概要

BTYDモデルとは、顧客の過去の購買行動データを基に、将来の購買行動を予測する統計モデルです。購買行動は Customer Lifetime Value（CLV）によって定量化されます。ここでCLVとは、将来の任意の期間における購買回数と平均購買金額の積から算出される合計購買金額を指します。

また、BTYDモデルは階層ベイズモデルから構成されており、顧客個人差を考慮したモデリングが可能です。ここで階層ベイズモデルとは、出力分布のパラメータに「事前分布」を設定するベイズモデルを指します。

モデル構造

BTYDモデルは、次の2つの統計モデルを組み合わせた構造になっています。

1. Beta Geometric（BG） / Negative Binomial Distribution（NBD）モデル

・BGモデル： 顧客が購買を離脱する確率（離脱率）をモデル化します。

　　例：「この顧客はすでに離脱しているのか？」

・NBDモデル： 顧客の購買回数をモデル化します。

　　例：「顧客Aは今後3か月間で何回購入するか？」

2. Gamma-Gammaモデル

顧客が1回の購買で支払う金額（平均購買金額）をモデル化します。

　　例：「顧客Bは1回の購入で平均いくら支払うか？」

BG/NBDモデル

NBDモデル：購買回数の予測

モデル構造

出力分布：ポアソン分布

概要

平均発生回数が $λ$ である事象が、ある期間内に $x$ 回発生する確率をモデリングする確率分布です。

　　例：1か月に平均2回購入する顧客が、3か月間で5回購入する確率

NBDモデルでは、観測期間 $(0,T]$ において、平均購買回数 $λ$ に従う顧客の購買回数が $x$ である確率を表現します。

確率質量関数（PMF）

$$
f_{\text{poisson}}(x;\lambda,T) = \frac{(\lambda T)^x e^{-\lambda T}}{x!}, \quad (x = 0, 1, 2, \ldots)
$$

変数

$x$：購買回数

$T$：観測期間

パラメータ

$\lambda$：平均購買回数

　　・分布の平均と分散を制御します。

　　・値が大きいほど、購買回数 $x$ の平均値と分布のばらつきも拡大します。

　　・値が小さいほど、購買回数 $x$ の平均値と分布のばらつきも縮小します。

事前分布：ガンマ分布

概要

$\alpha$ 時間あたりに $r$ 回発生する確率をモデリングする確率分布です。

ポアソン分布の平均発生回数 $λ$ に対する事前情報を提供する確率分布としても利用可能です。

NBDモデルでは、顧客間の平均購買回数 λ の個人差を表現します。

確率密度関数（PDF）

$$f_{\text {gamma}}(\lambda; r,\alpha) = \frac{\alpha^r}{\Gamma(r)} \lambda^{r – 1} e^{-\alpha \lambda}, \quad (\lambda > 0)$$

パラメータ

$r$：形状パラメータ

　　・分布の形状を制御します。

　　・値が大きいほど、分布は正規分布のように左右対称の形に近づくため、平均購買回数 λの平均値と分散が大きくなり、購買回数の個人差が大きくなります。

　　・値が小さいほど、分布は0付近に集中するため、平均購買回数 λの平均値と分散が小さくなり、購買回数の個人差が小さくなります。

$\alpha$：レートパラメータ

　　・分布の減衰速度を制御します。

　　・値が大きいほど分布は急速に減衰し分布の集中度が高くなるため、平均購買回数 λの平均値と分散が小さくなり、購買回数の個人差が小さくなります。

　　・値が小さいほどなだらかな減衰になり分布の裾野が長くなるため、平均購買回数 λの平均値と分散が大きくなり、購買回数の個人差が大きくなります。

NBDモデル

平均発生回数 $λ$ に対し周辺化することで、ガンマ分布を事前分布として有するポアソン分布混合モデルの確率密度関数を導出できます。

周辺化した数式は、負の二項分布（Negative Binomial Distribution：NBD）の形状に一致します。

$$\begin{align*}f_{\text {nbd}}(x;r,\alpha)&=\int_{0}^{\infty}f_{\text {poisson}}(x;\lambda,T)f_{\text {gamma}}(\lambda; r,\alpha)\,d\lambda\\&=\frac{\Gamma(r+x)}{\Gamma(r)\Gamma(x+1)}\Big(\frac{\alpha}{\alpha+T}\Big)^r\Big(\frac{T}{\alpha+T}\Big)^x\end{align*}$$

BGモデル：離脱率の予測

モデル構造

出力分布：幾何分布

概要

発生率 $\rho$ に従う事象が、初めて発生するまでの試行回数 $x$ をモデリングする確率分布です。

BGモデルでは、離脱率 $\rho$ に従う顧客がちょうど $x$ 回目の購買直後に離脱する確率を表現します。

確率質量関数（PMF）

$$f_{\text {geometric}}(x;\rho) = (1 – \rho)^{x – 1} \rho, \quad (x = 1, 2, 3, \ldots)$$

パラメータ

$\rho$：離脱率

　　・値が大きいほど分布は 0 付近に集中し、急激に減衰します。

　　・離脱率が高いほど、多く試行回数を重ねる前に少ない試行回数で離脱が発生します。

　　・値が小さいほど分布は全範囲に分散し、緩やかに減衰します。

　　・離脱率が低いほど、多く試行回数を重ねた後に離脱が発生する可能性が高くなります。

事前分布：ベータ分布

概要

幾何分布の事象の発生率 $\rho$ に対する事前情報を提供する確率分布です。

BGモデルでは、顧客間の離脱率 $\rho$ の個人差を表現します。

確率密度関数（PDF）

$$f_{\text {beta}}(\rho; a, b) = \frac{\Gamma(a + b)}{\Gamma(a) \Gamma(b)} \rho^{a – 1} (1 – \rho)^{b – 1}, \quad (0 < \rho < 1)$$

パラメータ

$a$：形状パラメータ

　　・値が大きいほど分布は1付近に集中するため、離脱率 $\rho$ が高い顧客が多くなります。

　　・値が小さいほど分布は0付近に集中するため、離脱率 $\rho$ が低い顧客が多くなります。

$b$：形状パラメータ

　　・値が大きいほど分布は0付近に集中するため、離脱率 $\rho$ が低い顧客が多くなります。

　　・値が小さいほど分布は1付近に集中するため、離脱率 $\rho$ が高い顧客が多くなります。

補足

　　・$a>b$：離脱率 $\rho$ が高い顧客が多くなります。

　　・$a=b$：左右対称で離脱率 $\rho=0.5$ 付近に最も確率が集中します。

　　・$a<b$：離脱率 $\rho$ が低い顧客が多くなります。

BGモデル

離脱率 $\rho$ に対し周辺化することで、ベータ分布を事前分布に有する幾何分布混合モデルの確率密度関数を導出することができます。

考慮すべき事象は下記の2パターンが考えられます。

1. $x$回目の購買後にも生存している確率：顧客データの最終購買日以降もアクティブ状態である確率

$$\begin{align*}f_{\text {bg}_{\text{active}}}(x;a,b)&=\int_{0}^{1}f_{\text{geometric}}(x;\rho)f_{\text{beta}}(\rho; a, b)\,d\rho\\&=\int_{0}^{1}(1-\rho)^x\frac{\Gamma(a + b)}{\Gamma(a) \Gamma(b)} \rho^{a – 1} (1 – \rho)^{b – 1}\,d\rho\\&=\frac{B(a,b+x)}{B(a, b)}\end{align*}$$

2. $x$回目の購買後に離脱する確率：顧客データの最終購買日以降は非アクティブ状態になる確率

$$\begin{align*}f_{\text {bg}_{\text {inactive}}}(x;a,b)&=\int_{0}^{1}f_{\text{geometric}}(x;\rho)f_{\text{beta}}(\rho; a, b)\,d\rho\\&=\int_{0}^{1}\rho(1-\rho)^{x-1}\frac{\Gamma(a + b)}{\Gamma(a) \Gamma(b)} \rho^{a – 1} (1 – \rho)^{b – 1}\,d\rho\\&=\frac{B(a+1,b+x-1)}{B(a, b)}\end{align*}$$

パラメータ推定

尤度関数

ベイズモデルのパラメータ推定は、モデルの尤もらしさを表す尤度関数を最大化するパラメータを推定します。尤度関数は、パラメータを与えた時の観測データの同時確率から算出され、あるパラメータを設定した時の確率分布が、データの特徴をよく表している場合に尤度関数は大きい値を取り、あまりデータの特徴と一致しない場合に尤度は小さい値を取ります。

BG/NBDモデルの場合、考慮できる事象として「最終購買日後も生存している場合」と「最終購買と共に離脱する場合」の２パターンが考えられ、それらの和の確率により個々の顧客の尤度が算出されます。それらを全顧客に対し計算し、総積を取ることで全体の尤度を計算します。

$$\begin{align*}L(r,\alpha,a,b)&=\prod_{i=1}^{N}\Big\{{f_{\text{bg}_{\text{active}}}(x_i;a,b)f_{\text{nbd}}(x_i;r,\alpha)+f_{\text{bg}_{\text{inactive}}}(x_i;a,b)f_{\text{nbd}}(x_i;r,\alpha)}\Big\}\\&=\prod_{i=1}^{N}\Big\{\underbrace{\frac{B(a,b+x_i)}{B(a, b)}\frac{\Gamma(r+x_i)\alpha^r}{\Gamma(r)(\alpha+T_i)^{r+x_i}}}_{\text アクティブ状態の項}+\underbrace{\delta_{x_i>0}\frac{B(a+1,b+x_i-1)}{B(a, b)}\frac{\Gamma(r+x_i)\alpha^r}{\Gamma(r)(\alpha+t_{x_{i}})^{r+x_i}}\Big\}}_{\text 非アクティブ状態の項}\end{align*}$$

変数

$t_{x}$：生存期間（初回購買日から最終購買日）

$N$：合計顧客数

対数尤度関数

一般的に、演算の安定性を担保するために尤度関数を対数化した対数尤度関数に対し最適化を行います。下記の対数尤度関数に対し、MCMCサンプリング、及び、MAP推定を行いパラメータの推定を行います。MAP推定は、対数尤度関数を最大化するパラメータを推定するため「点推定」に当たり、MCMCサンプリングは、対数尤度関数を指標にパラメータをサンプリングするため、「事後分布の推定」に当たります。

$$\begin{align*}LL(r,\alpha,a,b)&=\sum_{i=1}^{N}\Bigg[\log\Gamma(r+x_i)-\log\Gamma(r)+\log\Gamma(a+b)\\&+\log\Gamma(b+x_i)-\log\Gamma(b)-\log\Gamma(a+b+x_i)\\&+r\log\alpha-(r+x_i)\log(\alpha+t_{x_{i}})\\&+\log\Big\{\big(\frac{\alpha+t_{x_{i}}}{\alpha+T_{i}}\big)^{r+x_i}+\delta_{x_i>0}\frac{a}{b+x_i-1}\Big\} \Bigg]\end{align*}$$

予測

予測ステップでは、推定されたパラメータ $\{ \hat{r}, \hat{\alpha}, \hat{a}, \hat{b} \}$ を用いて購買回数 $X(t)$ の期待値 $\mathbb{E}(X(t))$を計算します。

顧客全体における購買回数の期待値

予測期間 $t$ とした時の、期間$(T,T+t]$における顧客全体の購買回数の期待値 $\mathbb{E}(X(t)|\hat{r},\hat{\alpha},\hat{a},\hat{b})$ は下記で定義されます。

$$\begin{align*}\mathbb{E}(X(t)|\hat{r},\hat{\alpha},\hat{a},\hat{b}) &= \int_{0}^{\infty}\int_{0}^{1}\mathbb{E}(X(t)|\lambda,\rho)f_{\text{gamma}}(\lambda;\hat{r},\hat{\alpha})f_{\text{beta}}(\rho;\hat{a},\hat{b})\,d\rho\,d\lambda \\&= \frac{\hat{a}+\hat{b}-1}{\hat{a}-1} \Bigg\{1-\Big(\frac{\hat{\alpha}}{\hat{\alpha}+t}\Big)^{\hat{r}}{}_2F_{1}\Big(\hat{r},\hat{b};\hat{a}+\hat{b}-1;\frac{t}{\hat{\alpha}+t}\Big)\Bigg\}\end{align*}$$

変数

$t$：予測期間

推定後のパラメータ：$\{ \hat{r}, \hat{\alpha}, \hat{a}, \hat{b} \}$

導出概要

パラメータ $\lambda$ 、$\rho$ が与えられた時、購買回数の期待値 $\mathbb{E}(X(t)|\lambda,\rho)$ は下記で定義されます。

$$\begin{align*}\mathbb{E}(X(t)|\lambda,\rho)&=\lambda t\cdot P(\tau>t)+\int_{0}^{t}\lambda \tau \cdot g(\tau;\lambda,\rho)\,d\tau\\&=\lambda t\cdot e^{-\rho\lambda t}+\lambda^2\rho\int_{0}^{t}\tau e^{-\lambda \rho t}\,dt\\&=\frac{1}{\rho}\Big(1-e^{-\lambda \rho t}\Big)\end{align*}$$

ここで、$j$ 回目の購買を行い、なおかつ離脱しない確率 $P(\tau>t)$ は下記で定義されます。

$$\begin{align*}P(\tau>t)&=\sum_{j=0}^{\infty}(1-\rho)^j\cdot\frac{(\lambda t)^j e^{-\lambda t}}{j!}\\&=e^{-\rho\lambda t}\end{align*}$$

また、ちょうど時刻 $\tau$ で離脱する確率 $g(\tau;\lambda,\rho)$を下記で定義されます。

$$g(\tau;\lambda,\rho)=\lambda p e^{-\lambda \rho t}$$

ただし $\lambda$ と $\rho$ は観測できないため、事前分布 $f_{\text {gamma}}({\lambda}; \hat{r},\hat{\alpha})$と$f_{\text {beta}}(\rho; \hat{a}, \hat{b})$を与え、$\lambda$ と $\rho$ に対して周辺化することで期待値 $\mathbb{E}(X(t)|\hat{r},\hat{\alpha},\hat{a},\hat{b})$ を算出します。

$$\begin{align*}\mathbb{E}(X(t)|\hat{r},\hat{\alpha},\hat{a},\hat{b}) &= \int_{0}^{\infty}\int_{0}^{1}\mathbb{E}(X(t)|\lambda,\rho)f_{\text{gamma}}(\lambda;\hat{r},\hat{\alpha})f_{\text{beta}}(\rho;\hat{a},\hat{b})\,dp\,d\lambda \\&= \frac{\hat{a}+\hat{b}-1}{\hat{a}-1} \Bigg\{1-\Big(\frac{\hat{\alpha}}{\hat{\alpha}+t}\Big)^{\hat{r}}{}_2F_{1}\Big(\hat{r},\hat{b};\hat{a}+\hat{b}-1;\frac{t}{\hat{\alpha}+t}\Big)\Bigg\}\end{align*}$$

特定の顧客における購買回数の期待値

観測期間$(0,T_i]$において購買回数が$x_i$、生存期間が$t_{x_{i}}$である特定の顧客のデータが与えられた時、予測期間 $t$ とした時の期間$(T_i,T_i+t]$におけるその顧客の購買回数の期待値 $\mathbb{E}(X_i(t) \mid x_i, t_{x_i}, T_i; \hat{r}, \hat{\alpha}, \hat{a}, \hat{b})$ は下記で定義されます。

$$\begin{align*}\mathbb{E}(X_i(t) \mid x_i, t_{x_i}, T_i; \hat{r}, \hat{\alpha}, \hat{a}, \hat{b}) &= \int_{0}^{\infty} \int_{0}^{1} \mathbb{E}(X_i(t)|x_i,t_{x_{i}},T_i;\lambda,\rho)f(\lambda,\rho|x_i,t_{x_{i}},T_i;\hat{r},\hat{\alpha},\hat{a},\hat{b})\,d\rho\,d\lambda \\&= \frac{\frac{\hat{a} + \hat{b} + x_i – 1}{\hat{a} – 1}\Big\{1 – \Big(\frac{\hat{\alpha} + T_i}{\hat{\alpha} + T_i + t}\Big)^{\hat{r} + x_i} {}_2F_{1}\Big(\hat{r} + x_i, \hat{b} + x_i; \hat{a} + \hat{b} + x_i – 1; \frac{t}{\hat{\alpha} + T_i + t}\Big)\Big\}}{1 + \delta_{x_i > 0}\frac{\hat{a}}{\hat{b} + x_i – 1} \Big(\frac{\hat{\alpha} + T_i}{\hat{\alpha} + t_{x_i}}\Big)^{\hat{r} + x_i}}\end{align*}$$

導出概要

顧客データ $\{x_i,t_{x_{i}},T_i\}$ が与えられた時の購買回数の期待値 $\mathbb{E}(X_i(t)|x_i,t_{x_{i}},T_i;\lambda,\rho)$ は下記で定義されます。

$$\mathbb{E}(X_i(t)|x_i,t_{x_{i}},T_i;\lambda,\rho)=\mathbb{E}(X_i(t)|\lambda,\rho) \cdot P(\text {active at }T_i|x_i,t_{x_{i}},T_i,\lambda,\rho)$$

ここで、パラメータ $\lambda$ 、$\rho$ が与えられた時、購買回数の期待値 $\mathbb{E}(X_i(t)|\lambda,\rho)$ は下記で定義されます。

$$\mathbb{E}(X_i(t)|\lambda,\rho)=\frac{1}{\rho}\Big(1-e^{-\lambda \rho t}\Big)$$

また、顧客データ $\{x_i,t_{x_{i}},T_i\}$ とパラメータ $\lambda$ 、$\rho$ が与えられた時の $x$ 回目の購買後の生存率は下記で定義されます。

$$P(\text {active at }T_i|x_i,t_{x_{i}},T_i,\lambda,\rho)=\frac{(1-\rho)^x \lambda^x e^{-\lambda T}}{L(\lambda,\rho|x_i,t_{x_{i}},T_i)}$$

ここで、尤度 $L(\lambda,p|x_i,t_{x_{i}},T_i)$ は下記で定義されます。

$$L(\lambda,\rho|x_i,t_{x_{i}},T_i)=(1-\rho)^x\lambda^x e^{-\lambda T_i}+\delta_{x>0}\rho(1-\rho)^{x-1}\lambda^x e^{-\lambda t_{x_{i}}}$$

ただし $\lambda$ と $\rho$ は観測できないため、事後分布 $f(\lambda,\rho|x_i,t_{x_{i}},T_i;\hat{r},\hat{\alpha},\hat{a},\hat{b})$ から $\lambda$ と $\rho$ を得ます。ベイズの定理より事後分布 $f(\lambda,\rho|x_i,t_{x_{i}},T_i;\hat{r},\hat{\alpha},\hat{a},\hat{b})$ は下記で定義されます。

$$\begin{align*}f(\lambda,\rho|x_i,t_{x_{i}},T_i;\hat{r},\hat{\alpha},\hat{a},\hat{b})=\frac{L(\lambda,\rho|x_i,t_{x_{i}},T_i)f_{\text {gamma}}(\lambda|x_i,t_{x_{i}},T_i;\hat{r},\hat{\alpha})f_{\text {beta}}(\rho|x_i,t_{x_{i}},T_i;\hat{a},\hat{b})}{\int_{0}^{\infty}\int_{0}^{1}L(\lambda,\rho|x_i,t_{x_{i}},T_i)f_{\text {gamma}}(\lambda|x_i,t_{x_{i}},T_i;\hat{r},\hat{\alpha})f_{\text {beta}}(\rho|x_i,t_{x_{i}},T_i;\hat{a},\hat{b})}\end{align*}$$

購買回数の期待値 $\mathbb{E}(X_i(t)|x_i,t_{x_{i}},T_i;\lambda,\rho)$ と 事後分布 $f(\lambda,\rho|x_i,t_{x_{i}},T_i;\hat{r},\hat{\alpha},\hat{a},\hat{b})$ に対し、 $\lambda$ と $\rho$ で周辺化することで、顧客データ $\{x_i,t_{x_{i}},T_i\}$ 、及び推定パラメータ $\{ \hat{r}, \hat{\alpha}, \hat{a}, \hat{b} \}$ が与えられた時の購買回数の期待値 $\mathbb{E}(X_i(t) \mid x_i, t_{x_i}, T_i; \hat{r}, \hat{\alpha}, \hat{a}, \hat{b})$ は下記で定義されます。

$$\begin{align*}\mathbb{E}(X_i(t) \mid x_i, t_{x_i}, T_i; \hat{r}, \hat{\alpha}, \hat{a}, \hat{b}) &= \int_{0}^{\infty} \int_{0}^{1} \mathbb{E}(X_i(t)|x_i,t_{x_{i}},T_i;\lambda,\rho)f(\lambda,\rho|x_i,t_{x_{i}},T_i;\hat{r},\hat{\alpha},\hat{a},\hat{b})\,d\rho\,d\lambda \\&= \frac{\frac{\hat{a} + \hat{b} + x_i – 1}{\hat{a} – 1}\Big\{1 – \Big(\frac{\hat{\alpha} + T_i}{\hat{\alpha} + T_i + t}\Big)^{\hat{r} + x_i} {}_2F_{1}\Big(\hat{r} + x_i, \hat{b} + x_i; \hat{a} + \hat{b} + x_i – 1; \frac{t}{\hat{\alpha} + T_i + t}\Big)\Big\}}{1 + \delta_{x_i > 0}\frac{\hat{a}}{\hat{b} + x_i – 1} \Big(\frac{\hat{\alpha} + T_i}{\hat{\alpha} + t_{x_i}}\Big)^{\hat{r} + x_i}}\end{align*}$$

計算例

変数の設定

　　・パラメータ

　　　　・NBDモデル

　　　　　　・形状パラメータ：$\hat{r}=0.8$

　　　　　　・レートパラメータ：$\hat{\alpha}=4$

　　　　・BGモデル

　　　　　　・形状パラメータ：$\hat{a}=1.2$

　　　　　　・形状パラメータ：$\hat{b}=2.0$

顧客データ

　　・購買回数：$x_i=4$（回）

　　・生存期間：$t_{x_{i}}=60$（日）

　　・観測期間：$T_i=90$（日）

　　・予測期間：$t=365$（日）

予測

　　・分子

$$\begin{align*}\text {numerator}&=\frac{1.2+2.0+4-1}{1.2-1}\Bigg\{1-\Big(\frac{4+90}{4+90+365}\Big)^{0.8+4}\\&\cdot{}_2F_1\Big(0.8+4,2.0+4;1.2+2.0+4-1;\frac{365}{4+90+365}\Big)\\&=30.94\end{align*}$$

　　・分母

$$\text {denominator}=1+\frac{1.2}{2.0+4-1}\Bigg(\frac{4+90}{4+60}\Bigg)^{0.8+4}=2.52$$

　　・購買回数の期待値

$$\mathbb{E}(X_i(t) \mid x_i, t_{x_i}, T_i; \hat{r}, \hat{\alpha}, \hat{a}, \hat{b}) =\frac{\text {numerator}}{\text {denominator}}=12.26（回/年）$$

Gamma-Gammaモデル

Gamma-Gammaモデル：平均購買金額の予測

※ 前提として、BTYDモデルでは顧客の「購買回数」と「平均購買金額」は独立であると仮定します。

モデル構造

出力分布：ガンマ分布

概要

NBDモデルの事前分布と同一の確率分布で、右に長い裾野を持ち、常に正の連続値を取る確率分布です。

Gamma-Gammaモデルでは、ある顧客の平均購買金額が $\overline{z}$ になる確率を表現します。

確率密度関数（PDF）

$$g_{\text {gamma}_{1}}(\overline{z};p,\nu;x)=\frac{(\nu x)^{px}\overline{z}^{px-1}e^{-\nu x\overline{z}}}{\Gamma(px)}$$

変数

$\overline{z}$ ：平均購買金額

$p$：形状パラメータ

　　・値が大きいほど、分布は平均値付近に集中し、ばらつきが小さくなります。

　　・値が小さいほど、分布は広がり、ばらつきが大きくなります。

$\nu$：レートパラメータ

　　・値が大きいほど、分布は急速に減衰します。

　　・値が小さいほど、分布の減衰は緩やかになり、分散が大きくなります。

事前分布：ガンマ分布

概要

Gamma-Gammaモデルでは、顧客の平均購買金額 $\overline{z}$ のばらつきの個人差を表現します。

確率密度関数（PDF）

$$g_{\text {gamma}_{2}}(\nu;q,v)=\frac{v^{q}\nu^{q-1}e^{-\nu v}}{\Gamma(q)}$$

変数

$q$：形状パラメータ

　　・値が大きいほど、出力分布のレートパラメータ $\nu$ の分布は尖鋭化し、個人差は小さくなります。

　　・値が小さいほど、出力分布のレートパラメータ $\nu$ の分布の裾が広がり個人差が大きくなります。

$v$：レートパラメータ

　　・値が大きいほど、レートパラメータ $\nu$ の分布は0付近で急速に減衰するため、平均購買金額レベルが高くなり、個人差のばらつきも小さくなります。

　　・値が小さいほど、レートパラメータ $\nu$ は緩やかに減衰するため、平均購買金額が安くなり、個人差のばらつきも大きくなります。

Gamma-Gammaモデル

レートパラメータ $\nu$ に対し周辺化することで、ガンマ分布を事前分布として有するガンマ分布の確率密度関数を導出できます。

$$\begin{align*}g_{\text {gg}}(\overline{z};p,q,v;x)&=\int_{0}^{\infty}g_{\text {gamma}_{1}}(\overline{z};p,\nu)g_{\text {gamma}_{2}}(\nu;q,v)\,d\nu\\&=\frac{\Gamma(px+q)}{\Gamma(px)\Gamma(q)}\frac{\overline{z}^{px-1}x^{px}v^q}{(x\overline{z}+v)^{px+q}}\end{align*}$$

パラメータ推定

尤度関数

$$L(p,q,v)=\prod_{i=1}^{N}g_{\text {gg}}(\overline{z_i};p,q,v;x_i)$$

対数尤度関数

$$\begin{align*}LL(p,q,v)&=\sum_{i=1}^{N}\text{log}g_{\text {gg}}(\overline{z_i};p,q,v;x_i)\\&=\sum_{i=1}^{N}\Big\{\log\Gamma(px_i+q)-\log\Gamma(px_i)-\log\Gamma(q)\\&+q\log v+px_i\log x_i+(px_i-1)\log \overline{z_i}-(px_i+q)\log(x_i\overline{z_i}+v)\Big\}\end{align*}$$

予測

予測ステップでは、推定されたパラメータ $\{ \hat{p}, \hat{q}, \hat{v}\}$ を用いて平均購買金額 $\bar{Z}$ の期待値 $\mathbb{E}(\bar{Z})$を計算します。

顧客全体における購買金額の期待値

顧客全体の購買金額の期待値 $\mathbb{E}(\overline{Z}|\hat{p},\hat{q},\hat{v})$ は下記で定義される。

$$\mathbb{E}(\overline{Z}|\hat{p},\hat{q},\hat{v})=\frac{\hat{p}\hat{v}}{\hat{q}-1}$$

導出概要

$\bar{Z}\sim \text{Gamma}(p,\nu)$ に従い、ガンマ分布の期待値の定義より平均購買金額 $\bar{Z}$ の期待値は下記で定義されます。

$$\mathbb{E}(\bar{Z})=\mathbb{E}\Big(\frac{\hat{p}}{\nu}\Big)$$

ここで平均購買金額 $\bar{Z}$ に対する潜在変数 $\zeta$ を導入します。

$$\zeta=\frac{\hat{p}}{\nu}より、\nu=\frac{\hat{p}}{\zeta}、\Big|\frac{d\nu}{d\zeta}\Big|=\Big|-\frac{\hat{p}}{\zeta^2}\Big|=\frac{\hat{p}}{\zeta^2}$$

また、レートパラメータ $\nu$ は下記で定義されます。

$$g_{\text {gamma}_{2}}(\nu;\hat{q},\hat{v})=\frac{\hat{v}^{\hat{q}}\nu^{\hat{q}-1}e^{-\nu \hat{v}}}{\Gamma(\hat{q})}$$

変数変換の法則より、潜在変数 $\zeta$ の分布 $f_{\text{inverse-gamma}}(\zeta;\hat{p},\hat{q},\hat{v})$ は下記の逆ガンマ関数で定義されます。

$$\begin{align*}f_{\text{inverse-gamma}}(\zeta;\hat{p},\hat{q},\hat{v}) &= g_{\text{gamma}_{2}}(\nu;\hat{q},\hat{v}) \cdot \left|\frac{d\nu}{d\zeta}\right| \\&= \frac{\hat{v}^{\hat{q}}\nu^{\hat{q}-1}e^{-\nu \hat{v}}}{\Gamma(\hat{q})} \cdot \frac{\hat{p}}{\zeta^2} \\&= \frac{(\hat{p}\hat{v})^{\hat{q}}\zeta^{\hat{q}-1}e^{-\frac{\hat{p}\hat{v}}{\zeta}}}{\Gamma(\hat{q})}\end{align*}$$

逆ガンマ関数の期待値の定義より、潜在変数 $\zeta$ の期待値 $\mathbb{E}(\zeta|\hat{p},\hat{q},\hat{v})$は下記で定義できます。

$$\mathbb{E}(\zeta|\hat{p},\hat{q},\hat{v})=\frac{\hat{p}\hat{v}}{\hat{q}-1}$$

特定の顧客における購買金額の期待値

観測期間$(0,T]$において購買回数が $x_i$ 、平均購買金額が $\overline{z_i}$ である特定の顧客のデータが与えられた時、その顧客の平均購買金額の期待値 $\mathbb{E}(\overline{Z_i}|\overline{z_i},x_i;\hat{p},\hat{q},\hat{v})$ は下記で定義されます。

$$\mathbb{E}(\overline{Z_i}|\overline{z_i},x_i;\hat{p},\hat{q},\hat{v})=\Big(\frac{\hat{q}-1}{\hat{p}x_i+\hat{q}-1}\Big)\Big(\frac{\hat{p}\hat{v}}{\hat{q}-1}\Big)\Big(\frac{\hat{p}x_i}{\hat{p}x_i+\hat{q}-1}\Big)\bar{z}_i$$

導出概要

顧客データ $\{x_i,\bar{z_i}\}$ と推定パラメータ $\{ \hat{p}, \hat{q}, \hat{v}\}$ 与えられた時、レートパラメータ $\nu$ の事後確率は下記で定義されます。

$$\begin{align*}f(\nu|\bar{z_i},x_i;\hat{p},\hat{q},\hat{v}) &= \frac{g_{\text {gamma}_{1}}(\overline{z_i};\hat{p},\nu;x_i)\cdot g_{\text {gamma}_{2}}(\nu;\hat{q},\hat{v})}{\int_{0}^{\infty}g_{\text {gamma}_{1}}(\overline{z_i};\hat{p},\nu;x_i)\cdot g_{\text {gamma}_{2}}(\nu;\hat{q},\hat{v})d\nu} \\&= \Bigg\{\frac{\nu^{\hat{p}x_i}\bar{z_i}^{\hat{p}x_i-1}e^{-\nu \bar{z_i}x_i}}{\Gamma(\hat{p}x_i)}\frac{\hat{v}^{\hat{q}}\nu^{\hat{q}-1}e^{-\nu \hat{v}}}{\Gamma(\hat{q})}\Bigg\}\Bigg/ \Bigg\{\frac{\Gamma(\hat{p}x_i+\hat{q})}{\Gamma(\hat{p}x_i)\Gamma(\hat{q})}\frac{\bar{z_i}^{\hat{p}x_i-1}x_i^{\hat{p}x_i}\hat{v}^{\hat{q}}}{(\hat{v}+x_i\bar{z_i})^{\hat{p}x_i+\hat{q}}}\Bigg\} \\&= \frac{(\hat{v}+x_i\bar{z_i})^{\hat{p}x_i+\hat{q}}\nu^{\hat{p}x_i+\hat{q}-1}e^{-\nu(\hat{v}+x_i\bar{z_i})}}{\Gamma(\hat{p}x_i+\hat{q})}\end{align*}$$

また、ここで潜在変数 $\zeta$ を下記のように定義します。

$$\zeta=\frac{1}{\nu}より、\nu=\frac{1}{\zeta}、\Big|\frac{d\nu}{d\zeta}\Big|=\Big|-\frac{1}{\zeta^2}\Big|=\frac{1}{\zeta^2}$$

変数変換の法則より、潜在変数 $\zeta$ の分布 $f_{\text{inverse-gamma}}(\zeta|\bar{z_i},x_i;\hat{p},\hat{q},\hat{v})$ は下記の逆ガンマ関数で定義されます。

$$\begin{align*}f_{\text{inverse-gamma}}(\zeta|\bar{z_i},x_i;\hat{p},\hat{q},\hat{v}) &= f(\nu|\bar{z_i},x_i;\hat{p},\hat{q},\hat{v}) \cdot \Bigg|\frac{d\nu}{d\zeta}\Bigg| \\&= \frac{(\hat{v}+x_i\bar{z_i})^{\hat{p}x_i+\hat{q}}\zeta^{-(\hat{p}x_i+\hat{q})-1}e^{-\frac{\hat{v}+x_i\bar{z_i}}{\zeta}}}{\Gamma(\hat{p}x_i+\hat{q})}.\end{align*}$$

逆ガンマ関数の期待値の定義より、潜在変数 $\zeta$ の期待値 $\mathbb{E}(\zeta|\overline{z_i},x_i;\hat{p},\hat{q},\hat{v})$は下記で定義できます。

$$\begin{align*}\mathbb{E}(\zeta|\overline{z_i},x_i;\hat{p},\hat{q},\hat{v})&=\frac{\hat{p}(\hat{v}+x_i\overline{z_i})}{\hat{p}x_i+\hat{q}-1}\\&=\Big(\frac{\hat{q}-1}{\hat{p}x_i+\hat{q}-1}\Big)\Big(\frac{\hat{p}\hat{v}}{\hat{q}-1}\Big)\Big(\frac{\hat{p}x_i}{\hat{p}x_i+\hat{q}-1}\Big)\bar{z}_i\end{align*}$$

計算例

変数の設定

　　・パラメータ

　　　　・ガンマ分布（出力分布）

　　　　　　・形状パラメータ：$\hat{p}=6$

　　　　・ガンマ分布（事前分布）

　　　　　　・形状パラメータ：$\hat{q}=5$

　　　　　　・形状パラメータ：$\hat{v}=2,000$

顧客データ

　　・購買回数：$x_i=4$（回）

　　・平均購買金額：$\overline{z_i}=5,000$（円）

予測

$$\mathbb{E}(\overline{Z_i}|\overline{z_i},x_i;\hat{p},\hat{q},\hat{v})=\frac{6\cdot(2,000+4\cdot5,000)}{6\cdot4+5-1}=4,714.29（円）$$

Customer Lifetime Value（CLV）

予測

顧客の将来の購買回数と購買金額を掛け合わせ、割引率を考慮してCLVを計算します。

将来の割引率を考慮した期間$(T,T+t]$におけるCLVの期待値

$$\text{CLV}_i = \sum_{t=1}^T \frac{\mathbb{E}(X_i(t) \mid x_i, t_{x_i}, T_i; \hat{r}, \hat{\alpha}, \hat{a}, \hat{b}) \cdot \mathbb{E}(\overline{Z_i}|\overline{z_i},x_i;\hat{p},\hat{q},\hat{v})}{(1 + d)^t}$$

変数

$\mathbb{E}(X_i(t) \mid x_i, t_{x_i}, T_i; \hat{r}, \hat{\alpha}, \hat{a}, \hat{b})$：将来の購買回数の期待値（BG/NBDモデル）

$\mathbb{E}(\overline{Z_i}|\overline{z_i},x_i;\hat{p},\hat{q},\hat{v})$：平均購買金額の期待値（Gamma-Gammaモデル）

$d$：割引率

割引率とは？

　・時間価値の減少

　　・「お金の現在価値（Present Value）」という概念に基づき、将来得られる収益は現在の収益よりも価値が低いと考えます。例えば、10年後の100万円は、現在の100万円と等価ではないため、割引係数を用いて調整します。

・リスクの考慮

　　・顧客が将来も継続して取引を行う保証はないため、将来の収益には不確実性が伴うことを考慮します。

計算例

変数の設定

　　・パラメータ

　　　　・年間割引率：$d=5$（%）

　　　　・日次割引率：$d_{\text {daily}}=(1+0.05)^{1/365}-1=0.0001337$（%）

　　　　・予測期間：$T=365$（日）

顧客データ

　　・年間購買回数の期待値：$\mathbb{E}(X_i(t)|x_i,t_{x_{i}},T_i;r,\alpha,a,b)=12.26$（回）

　　・購買金額の期待値：$\mathbb{E}(\overline{Z_i}|\overline{z_i},x_i;p,q,v)=4,714.29$（円）

予測

　　・日次購買回数

$$\lambda_i=\frac{\mathbb{E}(X_i(t)|x_i,t_{x_{i}},T_i;r,\alpha,a,b)}{365}=0.0336（回/日）$$

　　・日次期待購買金額

$$m_i=\lambda_i\cdot\mathbb{E}(\overline{Z_i}|\overline{z_i},x_i;p,q,v)=158.43（円/日）$$

　　・割引率総和

$$D=\sum_{t=1}^{365}\frac{1}{(1+d_{\text {daily}})^{t}}=\frac{1-(1+d_{\text {daily}})^{-365}}{d_{\text {daily}}}=364.8$$

　　・CLV

$$\text{CLV}_i=m_i\cdot D=57,912（円/年）$$

参考文献

・PyMC-Marketing Github Repository：https://github.com/pymc-labs/pymc-marketing

・RFM指標と顧客特性に基づくLTV予測モデル：https://www.jstage.jst.go.jp/article/pjsai/JSAI2019/0/JSAI2019_2Q1J202/_pdf/-char/ja

・“Counting Your Customers” the Easy Way:An Alternative to the Pareto/NBD Model：https://www.brucehardie.com/papers/bgnbd_2004-04-20.pdf

・Computing P(alive) Using the BG/NBD Model：https://www.brucehardie.com/notes/021/palive_for_BGNBD.pdf

・Overcoming the BG/NBD Model’s #NUM! Error Problem：https://brucehardie.com/notes/027/bgnbd_num_error.pdf

・A Step-by-Step Derivation of the BG/NBD Model：https://www.brucehardie.com/notes/039/bgnbd_derivation__2019-11-06.pdf

・The Gamma-Gamma Model of Monetary Value：https://www.brucehardie.com/notes/025/gamma_gamma.pdf

・RFM and CLV: Using Iso-value Curves for Customer Base Analysis：https://www.brucehardie.com/papers/rfm_clv_2005-02-16.pdf

ブログの著者欄

採用情報

関連記事

• 

  ローカル5G デモ / パブリックな4Gとローカル5Gの比較

  技術情報

• 

  SREが明かす！システム監視における動的閾値設定の適応例

  技術情報

• 

  「クラウドって色々あるけど、結局どう違うの？」AWS・GCP・Azure・GPUクラウドを解説！

  技術情報

• 

  ChatGPTを利用してTeamsでチャットボットを開発してみる

  技術情報

• 

  新卒技術研修 GMOテクノロジーブートキャンプ2021を実施しました【前編】

  技術情報

• 

  VM管理をもっと手軽に！WebVirtCloud利用ガイド

  技術情報