第１３章では系列データを扱うモデルとして隠れマルコフモデルが紹介されています。今回はその隠れマルコフモデル(HMM:Hidden Markov Model)での最尤推定を実装します。パラメータを固定してフォーワード-バックワードアルゴリズムを使って隠れ状態を推定するということは何度かしたことがあるのですが、パラメータまで推定するということはしたことなかったのでいい機会なのでやってみることにしました。ちなみに、HMMはカルマンフィルタや粒子フィルタなどの基礎になっています。

隠れマルコフモデル

隠れマルコフモデルでは系列データ（例えば、コインを百回投げて出た裏表）を考えます。そして、それらの観測された系列データの背後には潜在変数（例えば、裏（表）が出やすいコインが使われた）があるとしています。
隠れマルコフモデルのグラフィカルモデルは下の図(PRML図13.5より)のようになっています。

$\{x_i\}$を観測変数（コインの裏表）、$\{z_i\}$を潜在変数（コインの偏り）としています。潜在変数はひとつ前の潜在変数にのみ依存し（例：表が出やすいコインが使われた後は、また同じ表に偏ったコインが使われる）、観測変数はそのときの潜在変数にのみ依存しています（例：表に偏ったコインが使われたら、表が出やすい）。このように潜在変数がひとつ前の潜在変数にのみ依存している場合、一次マルコフ連鎖といいます。

数式を用いてもう少し細かく見ていきます。
潜在変数$z_n$が$z_{n-1}$にのみ依存していて、その関係を遷移確率行列$A$で表します。行列$A$のそれぞれの成分は潜在変数の遷移確率を表します。成分$A_{ij}$は$z_{n-1}$が状態$i$（$z_{n-1,i}=1$）のときに$z_n$の状態が$j$（$z_{n,j}=1$）になる確率となっています。ただし、潜在変数には一対K符号化法が用いられているとします。これらを数式にして書き下すと、

p(z_n|z_{n-1},A) = \prod_{i=1}^K\prod_{j=1}^KA_{ij}^{z_{n-1,i},z_{nj}}

となります。
ただし、一番最初の潜在変数$z_1$にはその前の潜在変数が存在しないので、K次元の確率ベクトル$\pi$を用いて

p(z_1|\pi) = \prod_{i=1}^K\pi_i^{z_{1i}}

とします。

観測変数$x$がD個の離散的な状態をとる離散変数（例：コインの裏表なら２つの状態がある離散変数）とした場合、観測変数の条件付き分布は

p(x_n|z_n,\mu) = \prod_{i=1}^D\prod_{k=1}^K\mu_{ik}^{x_{ni}z_{nk}}

となります。ただし、$x$にも一対K符号化法が用いられています。

HMMの最尤推定

上で紹介した隠れマルコフモデルでは３つのパラメータ$\pi,A,\mu$が登場してきました。ここでは、観測された系列データ$X = \{x_n\}_{n=1}^N$からそれら３つのパラメータを最尤推定します。これら３つのパラメータをひとまとめにして$\theta$と表すと、最大化すべき尤度関数は

p(X|\theta) = \sum_Zp(X,Z|\theta)

となります。ただし、$Z=\{z_n\}_{n=1}^N$。Zの場合の数は、$z$の状態数KのN乗となって（例えば：コインの偏りが２通りあってコインを１００回投げた場合、$2^{100}\approx10^{30}$通り）、直接計算できません。そこで、EMアルゴリズムの枠組みで一次の隠れマルコフモデルであることを活かして計算量を減らす方法を使います。

EMアルゴリズムでは、以下の式の計算と$\theta$についての最大化を繰り返します。

Q(\theta,\theta^{old}) = \sum_Zp(Z|X,\theta^{old})\ln p(X,Z|\theta)

Mステップ

$\gamma(z_n)$を潜在変数$z_n$の周辺事後分布、$\xi(z_{n-1},z_n)$を連続した潜在変数の同時事後分布とします。

\begin{align}
\gamma(z_n) &= p(z_n|X,\theta^{old})\\
\xi(z_{n-1},z_n) &= p(z_{n-1},z_n|\theta^{old})
\end{align}

これらを用いて、最大化すべき式を書き直すと、

Q(\theta,\theta^{old}) = \sum_{k=1}^K\gamma(z_{1k})\ln\pi_k + \sum_{n=2}^N\sum_{j=1}^K\sum_{k=1}^K\xi(z_{n-1,j},z_{n,k})\ln A_{jk} + \sum_{n=1}^N\sum_{k=1}^K \gamma(z_{nk})\ln p(x_n|z_n,\mu)

これをパラメータ$\pi,A,\mu$について最大化すると、そのときのパラメータはラグランジュ乗数法を用いて

\begin{align}
\pi_k &= {\gamma(z_{1k})\over\sum_{k=1}^K\gamma(z_{1j})}\\
A_{jk} &= {\sum_{n=2}^N\xi(z_{n-1,j},z_{nk})\over\sum_{l=1}^K\sum_{n=2}^N\xi(z_{n-1,j},z_{nl})}\\
\mu_{ik} &= {\sum_{n=1}^N\gamma(z_{nk})x_{ni}\over\sum_{n=1}^N\gamma(z_{nk})}
\end{align}

と求まります。

Eステップ

このステップでは、Mステップで用いた潜在変数についての２つの事後確率分布$\gamma(z_n),\xi(z_{n-1},z_n)$を計算します。その方法はフォーワード-バックワードアルゴリズムと言われていて、名前の通り系列データに関して前からと後ろからの２つの計算をすることで効率的に事後確率分布を計算できます。

フォーワード

$\alpha(z_1)=p(x_1,z_1)=p(x_1|z_1)p(z_1)$として、$\alpha(z_n)=p(x_1,\dots,x_n,z_n)$を再帰的に計算します。

\begin{align}
\alpha(z_n) &= p(x_1,\dots,x_n,z_n)\\
&= \sum_{z_{n-1}}p(x_1,\dots,x_n,z_n,z_{n-1})\\
&= \sum_{z_{n-1}}p(x_1,\dots,x_n|z_n,z_{n-1})p(z_n|z_{n-1})p(z_{n-1})\\
&= \sum_{z_{n-1}}p(x_1,\dots,x_{n-1}|z_{n-1})p(x_n|z_n)p(z_n|z_{n-1})p(z_{n-1})\\
&= p(x_n|z_n)\sum_{z_{n-1}}\alpha(z_{n-1})p(z_n|z_{n-1})
\end{align}

バックワード

$\beta(z_N)=p(x_N|z_N)$として、$\beta(z_n)=p(x_{n+1},\dots,x_N|z_n)$を再帰的に計算します。

\begin{align}
\beta(z_n)&=p(x_{n+1},\dots,x_N|z_n)\\
 &= \sum_{z_{n+1}}p(z_{n+1},x_{n+1},\dots,x_N|z_n)\\
&= \sum_{z_{n+1}}p(x_{n+1},\dots,x_N|z_n,z_{n+1})p(z_{n+1}|z_n)\\
&= \sum_{z_{n+1}}p(x_{n+2},\dots,x_N|z_{n+1})p(x_{n+1}|z_{n+1})p(z_{n+1}|z_n)\\
&= \sum_{z_{n+1}}\beta(z_{n+1})p(x_{n+1}|z_{n+1})p(z_{n+1}|z_n)
\end{align}

事後確率分布

上のフォーワードとバックワードより

\begin{align}
\gamma(z_n) &= {\alpha(z_n)\beta(z_n)\over p(X)}\\
\xi(z_{n-1},z_n) &= {\alpha(z_{n-1})p(x_n|z_n)p(z_n|z_{n-1})\beta(z_n)\over p(X)}
\end{align}

のように事後確率分布が求まります。

コード

import

いつもどおりnumpyとmatplotlibを使います。

import matplotlib.pyplot as plt
import numpy as np