f:id:vasilyjp:20180927090637j:plain

データサイエンティストの中村です。

webで発生するトランザクション（購買など）の中には、確率分布を仮定することで抽象化できる物があります。今回は、トランザクションが発生する現象をモデリングする手法のひとつであるBG/NBDモデルと、この手法にもとづいて将来発生するトランザクションの回数を予測するためのライブラリであるlifetimesを紹介します。

トランザクションのモデリングについて

1987年にSchmittlein等によってPareto/NBDというモデルが提案されました。これは顧客の継続的に発生する購買行動に確率分布を当てはめ抽象化する手法で、結果として将来発生する購買を予測することに成功しました。顧客が離脱したか否かの判断や顧客生涯価値の見積もりが可能になるという点で、Pareto/NBDモデルは顧客分析における非常に強力なツールのひとつです。

Pareto/NBDをベースとして改良されたモデルはいくつか存在し、Fader等によって提案されたBG/NBDモデルもそのひとつです。BG/NBDモデルはPareto/NBDに用いられたデータ表現にアレンジを加え、精度を損なわずシンプルに改良したモデルです。

本来これらのモデルは購買行動を分析対象として設計されたものですが、後述する仮定に嵌まるものであれば、購買行動以外の現象にも適用できます。以降はwebのデータを想定して、用語をユーザーとトランザクションに統一します。

BG/NBDモデル

データ

観測期間中に発生した、あるユーザーに由来するトランザクションの回数を $x$ とします（最初のトランザクションを0回目とカウントする）。 $j(≤x)$ 回目のトランザクションが発生した時刻を、0回目のトランザクションを起点として $t_j$ で表します。すなわち、観測期間中の最後のトランザクションは $t_x$ となります。また、0回目のトランザクションから観測終了までの時間を $T$ とおきます。

データに対する仮定

ユーザー個人に由来するトランザクションに対して、以下の仮定を敷きます。

i. ユーザーは生存/離脱を示す2値の状態を持ち、生存のときのみトランザクションが発生する。状態は観測されない。

最初のトランザクションが観測された時点で生存となり、この状態である限りはトランザクションがランダムに発生します。ユーザーがサービスを離脱した時、以降トランザクションは発生せず、このユーザーは再び生存状態になることはありません。生存/離脱を正確に見積もることはターゲティングの観点から非常に重要です。

ii. ユーザーが生存中は、観測期間中に発生するトランザクションの回数は単位時間を $t$ としてパラメータ $\lambda t$ のPoisson分布に従う。

同時に、Poisson分布と指数分布の関係から、トランザクションの時間間隔は指数分布に従うことが言えます。

$\displaystyle P(t_j|t_{j-1};\lambda) = \lambda e^{-\lambda (t_j - t_{j-1})}, \ \ \ \ t \gt t_{j-1} \geq 0$

$\lambda$ はモデルのパラメータで、トランザクションの発生率のような量を表します。

iii. パラメータ $\lambda$ はGamma分布に従う。

$\displaystyle P(\lambda | r, \alpha) = \frac{\alpha^{r} \lambda^{r-1} e^{-\lambda \alpha}}{\Gamma(r)}$

iv. トランザクションの発生後にユーザーが離脱する確率は一定で、 $p$ とする。

すなわち、 $j$ 回目のトランザクション後にユーザーが離脱する確率は幾何分布で表現できます。

$\displaystyle P(\text{inactive immediately afterafter $j$th transacton}) = p(1-p)^{j-1}$

正直、確率変数に $p$ を使いたくないですが、表記は元の論文に合わせます。

v. パラメータ $p$ はBeta分布に従う。

$\displaystyle P(p|a,b) = \frac{p^{a-1}(1-p)^{b-1}}{B(a,b)}$

そして、ユーザー全体に対する仮定として、ユーザーの特性を示すパラメータ $(λ,p)$ は独立とします。

$i$ をユーザーのインデックス、 $N$ をユーザー数として、モデルの入力は $\left\{ x_i, t_{x_i}, T_i \right\}_{i=1}^{N}$ 、ユーザー個人の行動を制御するパラメータは $\left\{ \lambda_i, p_i \right\}_{i=1}^{N}$ 、モデルのパラメータは $(r, \alpha, a, b)$ となります。

尤度関数の導出

本記事では導出過程を一部省略します。全体が気になる方はFader05をご覧ください。

観測期間中に $x$ 回のトランザクションが観測され、 $(t_x,T$ ]にトランザクションが発生していない時の尤度は、以下のようになります。

$x=0$ の時、 $(0,T$ ]にトランザクションが発生しないので、

$\displaystyle \mathcal{L}(x=0, T|\lambda) = 1-\int_{0}^{T} P(t|0;\lambda) dt = e^{-\lambda T}$

$x>0$ の時、 $t_x$ で離脱した場合と $T$ まで生存しながら $(t_x,T$ ]にトランザクションが発生しない場合を考慮して、

$\displaystyle \mathcal{L}(x, t_1,\ldots,t_x,T|\lambda,p) = p (1-p)^{x-1} \lambda^{x} e^{-\lambda t_x} + (1-p)^{x} \lambda^{x} e^{-\lambda T}$

上記2式をクロネッカーのデルタを使って書き直します。

$\displaystyle \mathcal{L}(x, t_x,T|\lambda,p) = (1-p)^{x} \lambda^{x} e^{-\lambda T} + \delta_{x>0} p (1-p)^{x-1} \lambda^{x} e^{-\lambda t_x}$

次に、パラメータ $\lambda, p$ の周辺化を試みます。事前分布を用いて

$\displaystyle \mathcal{L}(x, t_x, T|r, \alpha, a, b) \\ = \int \mathcal{L}(x, t_x,T|\lambda, p) P(\lambda|r, \alpha) P(p|a, b) d\lambda dp \\ =\frac{B(a, b+x)}{B(a, b)} \frac{\Gamma(r+x) \alpha^{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)(\alpha+t_x)^{r+x}}$

となります。式の見た目は屈強ですが計算は意外と簡単です。 $\lambda$ に関する積分と $p$ に関する積分は別々に計算できますし、確率密度関数と正規化項の関係を利用して楽ができます。

よって目的関数は以下のようになります。

$\displaystyle \mathop{\rm arg\,max}\limits_{r,\alpha,a,b} \sum_{i=1}^{N} \log \mathcal{L}(x_i, t_{x_i}, T_i|r, \alpha, a, b)$

推論時に重要な指標

Fader05ではトランザクションの回数に関する確率である $P(x|r,\alpha,a,b)$ やその期待値 $E[x|r,\alpha,a,b$ ]の導出が述べられています。式を書き写すだけになるので本記事では扱いません。

$\lambda,p$ の条件付き事後分布

それぞれ混合分布の形になると思います。私の手計算につき間違えている可能性もあるので、おかしな部分を発見した場合はご指摘ください。

$x=0$ のとき、

$\displaystyle P(\lambda|p,x,t_x,T;r,\alpha) = Gamma(r, T+\alpha) \\ P(p|\lambda,x,t_x,T;a,b) = Beta(a,b)$

$x>0$ のとき、

$\displaystyle P(\lambda|p,x,t_x,T;r,\alpha) = \frac{Z1_\lambda}{Z1_\lambda+Z2_\lambda} Gamma(x+r, T+\alpha) + \frac{Z2_\lambda}{Z1_\lambda+Z2_\lambda} Gamma(x+r, t_x+\alpha) \\ P(p|\lambda,x,t_x,T;a,b) = \frac{Z1_p}{Z1_p+Z2_p} Beta(a,x+b) + \frac{Z2_p}{Z1_p+Z2_p} Beta(a+1,x+b-1) \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Z1_{\lambda} = (1-p)^{x} (t_x+\alpha)^{x+r} \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Z2_{\lambda} = p (1-p)^{x-1} (T+\alpha)^{x+r} \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Z1_p = B(a,x+b) \lambda^{x}e^{-\lambda T} \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Z2_p = B(a+1,x+b-1) \lambda^{x}e^{-\lambda t_x}$

lifetimesについて

lifetimesはトランザクションデータからユーザーの生存/離脱を推定し、未来に発生するトランザクションやユーザーの生涯価値を予測するためのライブラリです。生データをほぼそのまま食わせることができるため非常に簡単に使えます。

github: https://github.com/CamDavidsonPilon/lifetimes
docs: http://lifetimes.readthedocs.io/en/latest/

公式のdocumentationにて使い方が丁寧に紹介されています。詳細はそちらを参照いただくとして、本記事ではかいつまんで紹介します。

ちなみに、本記事では購買データの代わりにIQONのコーデ投稿のログを用いています。

データの加工

ほとんど必要ありません。生ログからユーザーIDとdateのみ抽出してpandas.DataFrameにロードします。

In [4]: transactions.head()
Out[4]: 
   transaction_id  user_id        date
0          100001  2533078  2017-08-18
1          100002  2489799  2017-08-18
2          100003  2595561  2017-08-18
3          100004  2391692  2017-08-18
4          100005  2515961  2017-08-18

トランザクションの集計はlifetimes.utils.summary_data_from_transaction_dataやlifetimes.utils.calibration_and_holdout_dataにDataFrameを食わせるだけです。

In [5]: summary_calib_holdout = calibration_and_holdout_data(transactions, 'user_id', 'date',
   ...:                                                    calibration_period_end='2017-06-30',
   ...:                                                    observation_period_end='2017-12-31',
   ...:                                                    freq='D')
   ...: summary_calib_holdout.head()
   ...:                                                    
Out[5]: 
         frequency_cal  recency_cal  T_cal  frequency_holdout  \
user_id                                                         
7                  0.0          0.0   93.0                0.0   
77                 0.0          0.0  156.0                1.0   
260                0.0          0.0   93.0                0.0   
1468              11.0        132.0  140.0                0.0   
1938              40.0        134.0  175.0               12.0   

         duration_holdout  
user_id                    
7                     184  
77                    184  
260                   184  
1468                  184  
1938                  184

学習

scikit-learnと同じようにfitを実行します。

bgf = BetaGeoFitter(penalizer_coef=0.0)
bgf.fit(summary_calib_holdout.frequency_cal, summary_calib_holdout.recency_cal, summary_calib_holdout.T_cal)

サポートしているアルゴリズムはgithubをご確認ください。

可視化

可視化の為の関数は充実しています。例えば、キャリブレーションとホールドアウトそれぞれの期間で発生したトランザクションのトラジェクトリーを描画するには、以下の関数に学習済みモデルを渡すだけです。

plot_calibration_purchases_vs_holdout_purchases(bgf, summary_cal_holdout, n=20)

f:id:vasilyjp:20180319212921p:plain

オレンジがモデルの予測で青が実測です。モデルの予測と実際のトランザクション数が概ね一致している事が確認できます。

サンプルコード

というほどのものでもないですが、実験に使ったJupyter Notebookを載せておきます。後半では、周辺化した $\lambda,p$ をMCMCで復元してみました。モデルパラメータ $r,\alpha,a,b$ はユーザー全体の特徴を表現する値でしたが、 $\lambda,p$ を復元することでユーザー個人を分析することができます。

MCMCのスクリプトはNotebookと同じgistに公開しました。

モデルの制約

Pareto/NBDやBG/NBDモデルでは、

定常的に発生するトランザクションをモデル化したもので、周期性のある購買行動や広告などの影響は説明できない
一度離脱したユーザーは二度と復帰しない

というかなり強い制約があります。したがってこの手法を適用する前はデータの性質をよく観察しておく必要があります。

モデルの制約を緩める方向に発展させた手法は存在し、例えばRemoro13はユーザーの状態が変化することを許容し、HMMで状態遷移をモデル化しています。

まとめ

VASILYでは、最新の研究にアンテナを張りながら、同時にユーザーの課題解決を積極的に行うメンバーを募集しています。興味のある方はこちらからご応募ください。

https://www.wantedly.com/projects/109351www.wantedly.com

参考

David C. Schmittlein, Donald G. Morrison, Richard Colombo, Counting Your Customers: Who-Are They and What Will They Do Next?, 1987.
Peter S. Fader, Bruce G. S. Hardie, Ka Lok Lee, “Counting Your Customers” the Easy Way: An Alternative to the Pareto/NBD Model, 2005.
Jaime Romero, Ralf van der Lans, Berend Wierenga, A Partially Hidden Markov Model of Customer Dynamics for CLV Measurement, 2013.

ZOZO TECH BLOG

TECH BLOG

将来発生するトランザクション数を予測する方法