ill-identified diary

所属組織の見解などとは一切関係なく小難しい話しかしません

[stan][R] RFM分析と階層ベイズ法 (解決編)

概要

  1. 前回の『[python] [stan] 潜在変数と階層ベイズ法と RFM 分析 [未完成] - ill-identified diary』の完成版. 忙しくて1年近く放置していた……

  2. パラメータを推定し顧客ごとの生涯顧客価値 (CLV) の計算まで実行できた.

  3. stan は 2.14.0 を利用. 前回のは 2.9 で, 2.10 以降は構文が大きく変わっているので注意.

  4. 前回の「プログラム」以外のセクションを読んでからこちらを読むことをおすすめする.

  5. 文章量は4ページ (画像とプログラム除く)

反省点

実は, こちらですでに RF 分析についての stan の一部正解コードが書かれている.

abrahamcow.hatenablog.com

前回の一番の問題点は, 原理上離散的なパラメータを扱えないハミルトニアンモンテカルロ (以下, HMC) 法で\zeta を無理やり離散パラメータとして扱おうとしたことだろう. これは stan の仕様ではなく, 尤度関数の勾配ベクトルの計算が必要となる HMC 法 (とその改良版である NUTS アルゴリズム) の原理的な問題である. 対数尤度関数をパラメータで微分可能でなければ勾配ベクトルは定義できないので, 離散値をとるパラメータは HMC 法および NUTS では扱えない. 言語仕様の隙をつく裏ワザ的な方法であっても不可能なはずなのだ. 上記の記事や Gelman, A. et al. (2013) で指摘されているように, 今回のような, 2値の潜在変数などの離散パラメータを持つモデルを HMC で推定する場合, 離散パラメータを周辺化して打ち消すか, モデルをうまく離散パラメータと連続パラメータの部分とに分離して, 連続部分を HMC 法で, 離散部分を Gibbs サンプラーなど離散値を扱える他のアルゴリズムで計算する, などの方法がある. ただし, 後者は現状の stan では不可能である. というわけで stan でやるには指摘されているように各顧客が以後も生存するかどうかのパラメータ \(\zeta\) は周辺化で消去する必要がある.

しかし, 今回の問題はこれで終わりではない. 上記の記事では各パラメータが顧客に対し個別に設定されていない. つまり顧客の平均的な傾向として, パラメータを推定していることになる. 元の論文では, 顧客ひとりひとりの生涯顧客価値 (CLV) を計算するのが目的であり, そのためにパラメータを顧客個人レベルごとに用意しているから, stan で再現するにあたってもこれを守りたい.

修正版

モデルの解説

一応, 簡単モデルに解説する. モデルを p(R,F,M) というふうに表す. \(R\), \(F\), \(M\) はそれぞれ RFM に対応する確率分布と, その分布の形状を決めるパラメータと考えると, \(R\), \(F\) の部分が離散パラメータ \(\zeta\) に依存するので, \begin{align}
p(R,F,M)= & p(R,F|\zeta)p(\zeta)p(M)\\
\ln p(R,F,M)= & \ln\left[p(R,F|\zeta)p(\zeta)\right]+\ln p(M)\end{align}
と表現できるが, \zeta はゼロか1の2通りなので, その周辺化は \(p(\cdots)p(\zeta=1)+p(\cdots)p(\zeta=0)\) で簡単に表せる. よって対数尤度は

\begin{align}
\ln p(R.F,M)= & \ln p(R,F|\zeta)p(\zeta)+\ln p(M)\\
= & \ln\left(p(R,F|\zeta=1)p(\zeta=1)+p(\cdots|\zeta=0)p(\zeta=0)\right)+\ln p(M)\end{align}

となる. それぞれの対数確率は前回, 式で表したものをそのまま使える. まず, \zeta は, \begin{align}
p(\zeta=1)= & \dfrac{1}{1+\frac{\mu}{\lambda+\mu}\left\{ \exp\left[(\lambda+\mu)R\right]-1\right\} }\end{align}
であり, これをパラメータとするベルヌーイ分布を尤度に使えばよい. 続いて, \begin{align}
\ln p(R,F|\zeta=1)= & \ln\frac{\lambda^{x}t^{x-1}}{\Gamma(x)}\mathrm{e}^{-(\lambda+\mu)T}\\
= & \left[x\ln\lambda+(x-1)\ln t-\ln\Gamma(x)\right]-(\lambda+\mu)T,\\
\ln p(R,F|\zeta=0)= & \dfrac{\lambda^{x}t^{x-1}}{\Gamma(x)}\mu\mathrm{e}^{-(\lambda+\mu)\tau}\\
= & \left[x\ln\lambda+(x-1)\ln t-\ln\Gamma(x)\right]+\ln\mu-(\lambda+\mu)\tau,\\
\ln p(M)= & \ln\mathcal{N}(\ln\eta_{i},\,\omega^{2})\end{align}
となる. ただし,
\begin{align}
\left[x\ln\lambda+(x-1)\ln t-\ln\Gamma(x)\right]\end{align}
の部分は \zeta がいずれの値でも同じなので, stan のコード上では周辺化とは独立して記述できる. また, \(x\) はリピート購買回数を表すので, 1度しか買い物をしない顧客はゼロになる. \(\Gamma(0)\) は定義できないが, この場合は \(t=0\) つまり初めから生存していないという扱いになるので, 生存期間を表すガンマ分布に由来するこの部分の尤度を加える必要がない. よって, 周辺化消去は,
\begin{align}
\ln p(R,F|\zeta=1)= & -(\lambda+\mu)T,\\
\ln p(R,F|\zeta=0)= & \ln\mu-(\lambda+\mu)\tau\end{align}
p(\zeta) とで行う.

Stan コーディング

次に, stan のコードを書く際の技術的な問題を解説する.

ζ の扱い(39, 42~50行目)

まず確率 p(\zeta_{i}=1) を上記の式にしたがって generated parameters 項で,パラメータ \(\xi_{i}\)xi として計算する (42行目~) *1. xi は確率分布にしたがう乱数ではなく, 確率の大きさそのものを表している.

p(\zeta_{i}=1) の式を再確認すると, この確率 xi は, , つまり直近まで買い物をしたことがある, という状態ではかならず 1 となる. 整合のとれる定義としては問題ないのだが, このような関数では HMC 法に必要な勾配ベクトルが定義できなくなる. そこで, アップデートで追加された if 構文を活用して, \(R=0\) の場合のみ定数で1を代入するようにした.

τ の扱い (28, 38, 40, 66行目)

次に, \tau だが, これは \([t,T]\) の範囲に限定する必要がある. しかし, stan のパラメータの値の範囲制約は少し面倒で, パラメータの宣言部で real <lower=0, upper=1> tau ; のように宣言し, かつ対応する分布も切断分布などで範囲外の値を取らないようにする必要がある. stan の仕様上, パラメータが範囲外の値となった場合, 対応する対数尤度が \(-\infty\) となるからである*2. ところで, 今 \(\tau_{i}\) は顧客の数だけ必要だが, stan 2.14 時点では要素ごとに異なる制約を与えられない. 全て別の変数として羅列するのはあまりにも面倒なので, ややトリッキーだが以下のように対処した.

  1. (28行目) 各顧客に対応した \tau_{i} を 0 から 1 の範囲に統一した tau_rand として, 一様分布 \(U(0,1)\) から生成させる.

  2. (38, 40行目)generated parameters ブロックで \begin{align}
\tau_{i}= & t_{i}+(T_{i}-t_{i})\times\mathrm{tau\_rand}_{i}\end{align} という変換式で, 個々で異なる範囲 [t_{i},T_{i}] にある一様乱数\tau_{i} に変換する.

  3. (66行目) model ブロックで, ~ exponential() ではなく target += exponential_lpdf(tau ... ) を用いて \tau_{i} の尤度を与える*3.

この方法ならば, generated parameters での変換によって \tau_{i} は必ず \([t_{i},T_{i}]\) の範囲におさまる. ので, あとはこの値に対応した対数確率密度を与えればよい. さらに, target += ... で切断分布を与える場合は, さらに手動で \(\ln p(\cdot)=\ln f(x)-\ln(F(b)-F(a))\) の形になるよう調整する必要がある. しかし, \([t_{i},T_{i}]\) の範囲とするのは \(\zeta_{i}=0\) の場合だけで, それ以外は範囲制約は不要になる(ただし \tau>t_{i} にはなるはず). \(\zeta_{i}\) は周辺化消去しただけで値は確定していないため, 場合分けができない. そこで切断をしないことにした. 確率分布としては不適切だが, 尤度を評価するだけなら問題ないはず*4.

あとは, 再び if 文を使い, x_{i}>0 の時だけ,
\begin{align}
\left[x_{i}\ln\lambda_{i}+(x_{i}-1)\ln t_{i}-\ln\Gamma(x_{i})\right]\end{align}
の対数尤度を代入し, \zeta_{i} 絡みの部分を周辺化する.

なお, 事前分布は 松浦 (2016) を参考にして変更している.

最後に, generated quantities ブロックで CLV を計算する. また, 多変量正規分布を事前分布としたため, \lambda,\mu,\eta を1つのベクトルとしてさらに対数変換する必要があったので, 参照しやすいようここで lambda, mu, eta という名称で変換している.

以下がそのプログラム. データセットの加工は前回と全く同じなので省略した.

RFに対応する部分だけが解析的に事後分布を求められないため, 阿部 (2011) では RF 部分のみをメトロポリス-ヘイスティングス法で計算し, MAP推定値を得たうえで残りのMや線形回帰モデルの部分は共役性を利用してギブスサンプラーで計算しているようだが, このプログラムではRFMすべての要素を同時にHMC (NUTS) で推定している.


gist00cbfe03c6469a659d6f8c3d676da5a9

以下がRのプログラム. 上の stan ソースコードrfm_hierarchical1702.stan という名前で保存して呼び出している.


gist6690e6e75e9edd4ba4d1b0fa30007684

Rのプログラムにもいくつか変更を加えている まず, stan プログラムをコンパイルする際に, c++ のライブラリを使っている. stan_model() 関数の boost_lib=eigen_lib= はそのヘッダファイルの場所である. いくつかの処理が多少早くなるが, 導入が面倒な人はこれらの引数を削除して使うべし*5. また, これも最近知ったのだが, stan()sampling() では pars=include= を使うことで出力結果に残す変数を指定できるらしい. またしてもマニュアルはちゃんと読みましょうという話だった.

以上を実行したところ, 今度はうまく収束した. \boldsymbol{\theta}_{0},\boldsymbol{\Sigma}_{0} の収束は traceplot と Rhat で確認したが, 個別パラメータは全部見るのが面倒なので Rhat が 1.1未満であることをのみを確認している.プログラムは一応ggmcmcを使って全てのパラメータの traceplot とヒストグラムを描画しているが, これが意外と時間とメモリを要するので注意。


\theta_0
f:id:ill-identified:20170226230420p:plain
f:id:ill-identified:20170226230448p:plain

\Gamma_0
f:id:ill-identified:20170226230519p:plain
f:id:ill-identified:20170226230532p:plain
なお, 共分散行列\Gamma_0 はプログラム内部ではコレスキー分解を用いて処理しているため, 下側三角成分のみ意味を持ち, それ以外の成分は内部的にすべて1を与えている.

\lambda_iの一部
f:id:ill-identified:20170226231038p:plain

参考文献

*1:数値計算の精度を良くするために stan の特殊な関数をいくつか使用している. 詳しくは公式マニュアルを参照されたい.

*2:この仕様は stan の公式マニュアルに書いてあるため, ちゃんと読めば知っているはずである. ちゃんと読んでいなかったのでハマった.

*3:指数分布のパラメータはAbe (2009), 阿部 (2011) では\lambda+\mu としているが, 阿部 (2008) では単に \(\mu\) となっている. 前者は理論的な説明が不明瞭で, またパラメータの識別が難しくなるのではないかということから, 後者を採用した.

*4:いわゆる変則 (improper) 事前分布

*5:i7-4790K のマシンで1時間強かかるので, 性能の低いマシンでとりあえず試したいという場合は, 適当にデータを減らしてやってみるとよい.