1. 前回の振り返り

前回の記事では、顧客生涯価値(Life Time Value、以下LTV) の推定に広く利用されてきたBTYDモデルを紹介するとともに、その限界を指摘しました。その限界とは「顧客の一時的な離脱」をとらえることができないことだったので、それを捉えることのできるモデルを提案しました。以下この提案モデルを、BTYDTI モデルとでも呼ぶことにします(Buy Till You Die + Temporarily Inactive の省略です)。設定を再掲すると以下のようになります。

BTYDTI-1. 顧客にはデータでは直接観察されない「生存状態」、「死亡状態」に加え「飽き状態」の2つが存在する。
BTYDTI-2 全顧客は最初の観測時点において「生存状態」か「飽き状態」のいずれかあり、それぞれの確率を $\pi, 1-\pi$ と表現する。。
BTYDTI-3 生存状態であれば、毎日 $p$ の確率で購買を行うが、それ以外の状態では購買は発生しない。
BTYDTI-4 生存状態にある顧客は毎日 $s_1$ の確率で飽き状態、 $s_2$ の確率で死亡状態に移行する。飽き状態の顧客は確率 $r$ の確率で生存状態に復帰するが、死亡状態の顧客は生存状態にも飽き状態にも二度と復帰することはない。

　さて、この BTYDTI モデルは実は隠れマルコフモデルという枠組みで捉えることができます。したがって、今回は隠れマルコフモデルとは何か？そして、それはどのように推定が行われるのか？についてお話します。

2. 隠れマルコフモデルとは

　今、時間を $t$ の添え字で表すこととし、 $t = 1,2,...,T$ であるとします。我々が観測できる値(BTYDTIモデルにおける購買の有無)を $\boldsymbol{x} = (x_1,x_2,...,x_T)$ と表すことにし、それぞれの実現値を $1,2,...,M$ と離散的に表せるとします。例えばBTYDTIモデルでは購買が発生しなかったら1、発生したら2と表現するわけです。

さらに、観測できない状態（BTYDTIモデルでいうところの生存状態・飽き状態・死亡状態）を $\boldsymbol{q} = (q_1,q_2,...,q_T)$ で表すことにします。さらに、それぞれの状態を $1,2,...K$ と表すことにします。例えばBTYDTIモデルでは、生存状態を1、飽き状態を2、死亡状態を3とするわけです。

ここで、隠れマルコフモデルとは、これらの $(\boldsymbol{x}, \boldsymbol{q})$ を確率変数とみなしたモデルで、以下の2つの性質が成立することを仮定するモデルとなります。

仮定1-A. 「ある日の状態」は、「一時点前の状態」にのみ依存して決定される。すなわち、$t = 2, 3, ..., T$ について以下が成立する。

\begin{align}
P(q_t|q_{t-1},x_{t-1},q_{t-2},x_{t-2},...,q_1, x_1) = P(q_t|q_{t-1})
\end{align}

例えば、今日生存状態かどうかは「昨日の生存状態か、飽き状態か、死亡状態かのみに依存する」ということになります(BTYDTI-4 と整合的であることを確認してください)。一昨日この人がどんな状態だったか、であるとか、昨日この人は購買を行ったか、という情報は、「昨日の状態が何であったか」以上の情報を一切もたらさないということを仮定しています。

仮定1-B. $t=1$ に関しては一時点前の状態が存在しないため、他のどの観測値や状態にも依存せず独立に決定されると仮定する。

この想定は BTYDTI2 の前提と整合的であることを確認してください。

仮定2. 「ある日の観測値」は、「その日の状態」にのみ依存して決定される。すなわち、以下が成立する。

\begin{align}
P(x_t|q_T,x_T,q_{T-1},x_{T-1},...,q_{t+1},x_{t+1},q_{t},q_{t-1},x_{t-1},...,q_1,x_1) = P(x_t|q_t)
\end{align}

例えば、「今日購買を行うかどうか」は、「今日が生存状態か、飽き状態か、死亡状態」かということのみによって決まってくるということです（BTYDTI-3 と整合的であることを確認してください）。「昨日購買した」とか、「一昨日の状態は飽き状態でした」とかの情報は、「今日が何状態か」以上の情報を一切もたらさない訳です。

　さて、以上の２つの仮定を満たすモデルを隠れマルコフモデルというわけですが、今回はもう少し簡単にするため、さらに以下の仮定を課します :

仮定3. $P(q_t|q_{t-1})$、$P(x_t | q_t)$ は時間を通じて一定である。

この仮定も BTYDTI-3 や BTYDTI-4 と整合的であることを確認してください。以上のことから、BTYDTIモデルは 仮定1-A、仮定1-B、仮定2、仮定3 を満たす隠れマルコフモデルの一種として表現できるわけです。したがって、以下ではそうした、より一般的な隠れマルコフモデルの推定方法を示すことにし、その推定方法を応用することで、BTYDTIモデルの推定を行います。

　さて、推定に入る前に、この隠れマルコフモデルをもう少し具体化させておきましょう。まず、今日状態が $i$ である人が明日状態 $j$ になる確率を $A_{i,j}$ と表現することにし、これを $i$ 行 $j$ 列成分とする行列を $A$ とします。すなわち、

\begin{align}
A = 
\begin{bmatrix}
A_{1,1} & A_{1,2} & \cdots & A_{1,K} \\
A_{2,1} & A_{2,2} & \cdots & A_{2,K} \\
\vdots & \vdots & \ddots & \vdots \\
A_{K,1} & A_{K,2} & \cdots & A_{K,K}
\end{bmatrix}
\end{align}

また、$t = 1$ の時点の状態が $i$ である確率を $\pi_i$ とし、それが列ベクトルとして格納されたものを $\boldsymbol{\pi}$ と表します。すなわち、

\begin{align}
\boldsymbol{\pi} = 
\begin{bmatrix} 
\pi_1 \\ \pi_2 \\ \vdots \\ \pi_K 
\end{bmatrix}
\end{align}

さらに、今日状態が $i$ であることを前提として、今日の実現値が $m$ となる確率を $B_{i,m}$ と表し、それを $i$ 行 $m$ 列成分として格納した行列を $B$ とします。すなわち、

\begin{align}
B = 
\begin{bmatrix}
B_{1,1} & B_{1,2} & \cdots & B_{1,M} \\
B_{2,1} & B_{2,2} & \cdots & B_{2,M} \\
\vdots & \vdots & \ddots & \vdots \\
B_{K,1} & B_{K,2} & \cdots & B_{K,M}
\end{bmatrix}
\end{align}

少々冗長ではありますが、BTYDTIモデルにおいては、$A$、$\boldsymbol{\pi}$、$B$ はそれぞれ以下のように表現されることをご確認ください。

\begin{align}
A = 
\begin{bmatrix}
1 - s_1 - s_2 & s_1 & s_2 \\
r & 1-r & 0 \\
0 & 0 & 1
\end{bmatrix}
\end{align}

\begin{align}
\boldsymbol{\pi} = 
\begin{bmatrix}
\pi \\ 1-\pi \\ 0
\end{bmatrix}
\end{align}

\begin{align}
B = 
\begin{bmatrix}
1-p & p \\
1 & 0 \\
1 & 0
\end{bmatrix}
\end{align}

3. 隠れマルコフモデルの推定

さて、推定したいパラメータをまとめて $\boldsymbol{\theta}$ と表現することにしましょう。すると、任意の $(\boldsymbol{x},\boldsymbol{q})$ の実現確率は以下のように表現できます。

\begin{align}
P(\boldsymbol{x},\boldsymbol{q}|\boldsymbol{\theta}) &= \pi_{q_1} B_{q_1,x_1} \times A_{q_1,q_2} B_{q_2,x_2} \times A_{q_2,q_3} B_{q_3,x_3} \cdots \times A_{q_{T-1},q_{T}} B_{q_{T},x_{T}} \nonumber \\
& = \pi_{q_1} \times \prod_{t=1}^{T-1} A_{q_{t},q_{t+1}} \times \prod_{t=1}^{T} B_{q_t, x_t} \tag{3-1}
\end{align}

また、自然対数をとり $\boldsymbol{\theta}$ の関数と読み替えることで、対数尤度を以下のように求めることができる :

\begin{align}
LL(\boldsymbol{\theta}|\boldsymbol{x},\boldsymbol{q}) &= \log \biggl[ \pi_{q_1} \times \prod_{t=1}^{T-1} A_{q_{t},q_{t+1}} \times \prod_{t=1}^{T} B_{q_t, x_t}  \biggr] \\
& = \log(\pi_{q_1}) + \sum_{t=1}^{T-1} \log(A_{q_{t},q_{t+1}}) + \sum_{t=1}^{T} \log(B_{q_t,x_{t}}) \tag{3-2}
\end{align}

したがって、もしも ($\boldsymbol{x}$, $\boldsymbol{q}$) の両方の系列が観測できるのであれば、(3-2) について、この値を最大化するパラメータを選択することにより、最尤推定量 $\boldsymbol{\theta}^{MLE}$ を得ることができます。しかしながら、今回のケースはそれほど簡単ではありません。なぜならば、$\boldsymbol{q}$ の系列に関しては観測できないためです。よってこのケースでは、尤度は以下のように定義しなければなりません。

\begin{align}
L(\boldsymbol{\theta} | \boldsymbol{x}) &= P(\boldsymbol{x}|\boldsymbol{\theta}) \\
& = \sum_{q \in Q} P(\boldsymbol{q}|\boldsymbol{\theta}) P(\boldsymbol{x}, \boldsymbol{q} | \boldsymbol{\theta})
\end{align}

ただし、$Q$ は起こりうる状態系列を全て集めた集合です（式変形の際には結合確率分布の計算公式を使用しています）。つまり状態が見えないため、状態の実現確率で重みづけした期待尤度を最適化する必要があるわけです。

　さて、ではこれを最適化しよう...と行きたいところですが、これはそう簡単なことではありません。なぜならば、$Q$ は $K^T$ 通り（状態の数の観測期間乗）あるためです。いちいち起こりうるすべての $\boldsymbol{q}$ について任意の $\boldsymbol{\theta}$ の時の実現確率を求めて...と計算していたのでは、到底実用可能なものとはなりません。仮にそれらを求められたとしても、期待尤度を計算する際に和を取ってしまっているため、(3-2) のようにきれいな対数尤度を計算することができず、膨大な計算コストがかかってしまいます。

そこで、 近似的であってもなるべく効率的に 計算を行う必要があります。ここで登場するのが、EMアルゴリズムの一種である Baum and Welch アルゴリズムです。大まかにいうとこのアルゴリズムは以下のステップで計算を行うことになります。

手順1 : まず、適当に $\boldsymbol{\theta}_{old}$ を定める。
手順2 : 暫定的な $\boldsymbol{\theta} _{old}$ を用いて、$P(\boldsymbol{x},\boldsymbol{q} | \boldsymbol{\theta} _{old})$ を求める。この暫定的な確率を用いて、以下のように期待対数尤度を定義する :

\begin{align} 
ELL(\boldsymbol{\theta},\boldsymbol{\theta} _{old}) = \sum_{\boldsymbol{q} \in Q} P(\boldsymbol{x}, \boldsymbol{q}|\boldsymbol{\theta} _{old}) LL(\boldsymbol{\theta}) \tag{3-4}
\end{align}

手順3 : $ELL(\boldsymbol{\theta}, \boldsymbol{\theta} _{old})$ の最大化解 $\boldsymbol{\theta} ^{*}$ を求める。

手順4 : $\boldsymbol{\theta} _{old}$ を $\boldsymbol{\theta} ^*$ に更新する。

手順5 : 手順2 ~ 手順4 を $\boldsymbol{\theta} ^*$ が収束するまで繰り返す。

　こうした手法を採用することで、膨大な計算量を削減することができ、実用可能な範囲で隠れマルコフモデルを運用することができます。少々長くなってしまったので今回はここで記事を終了し、次回は Baum and Welch アルゴリズムの具体的な計算ステップを解説していきます。