今回は隠れマルコフモデル(Hidden Markov Model, HMM)についてです。ここではEMアルゴリズムによるパラメーターの最適化について取り上げます。 EMについてはGMMとEMアルゴリズム - 情報学関連の備忘録を参照。

概要

HMMは状態空間モデルの１つで、観測データを説明するようなマルコフ性をもつ離散的な隠れ状態を考えるモデルです。（ここでは１次のマルコフ連鎖だけを扱います。）観測データの尤度関数は $p ( {\bf X } | {\bf \theta} ) = \sum_{\bf Z} p ( {\bf X, Z} | {\bf \theta} )$ で与えられます。多くの場合では時系列データを考えるので、indexをtでおいて $\textbf{X} = ( x_1, \dots, x_t, \dots, x_T)$ とします。観測データの次元をN, 系列の長さをT, 隠れ変数の次元をKとします。 $z_t$ はbinaryの $K$ 次元ベクトルとします。どれか１つの次元 $k$ だけが1を取り、他は全て０になるような値を取ります。

f:id:endosan:20190730194224p:plain

このモデルにおける目標は、

$p ( {\bf X } | {\bf \theta} )$ を最大にするようなパラメーター ${\bf \theta}$ を求めること
その ${\bf \theta}$ を使って、 $p ( {\bf X } | {\bf \theta} )$ を最大にする隠れ状態の系列 ${\bf Z}$ を求めること

になります。

いくつかのパラメーターを導入します。

初期状態の確率 $\pi_k := p(z_{1k} = 1)$
遷移確率 $A_{jk} := p(z_{tk} = 1 | z_{t-1 j} = 1)$
出力確率のパラメーター $\phi$ 。出力確率 $p( {\bf x_t} | {\bf z_t, \phi} )$ が何かは場合による。

1. パラメーターの最適化

パラメーターの最適化をEMアルゴリズムで行います。以下のEM stepをデータの対数尤度 $p ( {\bf X } | {\bf \theta} )$ が収束するまで（または、パラメーターの変化量がある $\epsilon$ 以下になるまで）続けます。 EMはlocal maxmumにはまりがちなので、事前にこうなりそうだという予想があるならパラメーターの初期値を適切に選ぶ必要があります。

フォーワードバックワードアルゴリズム

前準備としてフォーワードバックワードアルゴリズム（HMMにおいては、Baum-Welch algorithmと呼ばれます。）について取り上げます。これを行うのは、E stepにおける計算に使用するため、かつ、データの対数尤度（収束判定に使いたい）を求めるのに使うからです。

まず、フォーワードアルゴリズムについて考えます。次のような変数を導入します。 $\alpha ( z_t ) := p(x_1, \dots, x_t, z_t | \theta)$ . この変数を丁寧に分解していくと以下のようになります。（グラフィカルモデルの方が直感的にわかりやすい） xについては分離されていることに注意。以下面倒なのでgivenなパラメーターを全て $\theta$ で書いてますが、ちゃんと書くと、 $A, \phi$ とかにしないといけないです。 \begin{align} \alpha (z_t ) &= p(x_1, \dots, x_t | z_t , \theta) p(z_t | \theta)\\ &= p(x_t | z_t, \theta) p(x_1, \dots, x_{t-1} | z_t , \theta) p(z_t | \theta)\\ &= p(x_t | z_t, \theta) p(x_1, \dots, x_{t-1}, z_t | \theta)\\ &= p(x_t | z_t, \theta) \sum_{z_{t-1}} p(x_1, \dots, x_{t-1}, z_{t-1}, z_t | \theta)\\ &= p(x_t | z_t, \theta) \sum_{z_{t-1}} p(x_1, \dots, x_{t-1}, z_t | z_{t-1}, \theta) p(z_{t-1} | \theta)\\ &= p(x_t | z_t, \theta) \sum_{z_{t-1}} p(x_1, \dots, x_{t-1} | z_{t-1}, \theta ) p( z_t | z_{t-1}, \theta) p(z_{t-1} | \theta)\\ &= p(x_t | z_t, \theta) \sum_{z_{t-1}} p(x_1, \dots, x_{t-1}, z_{t-1} |\theta) p( z_t | z_{t-1}, \theta)\\ &= p(x_t | z_t, \theta) \sum_{z_{t-1}} \alpha ( z_{t-1} ) p( z_t | z_{t-1}, \theta) \end{align} 成分ごとにみると（時刻tで隠れ状態kにいるとき、ただし正確な表記法ではない） \begin{align} \alpha ( z_{tk} =1) = p( {\bf x_t} | z_{tk}=1, {\bf \phi} ) \sum_{j=1}^{K} \alpha ( z_{t-1j} = 1) A_{jk} \end{align} この変数を使うと、観測データの尤度は \begin{align} p ( {\bf X } | {\bf \theta} ) = \sum_{z_{T}} \alpha (z_{T}) \end{align}

バックワードアルゴリズムについても同じように考えます。次のような変数を導入します。 $\beta (z_t) := p(x_{t+1}, \dots, x_{T} | z_t)$ . 上と同じように分解していくと、 \begin{align} \beta (z_t) = \sum_{z_{t+1}} \beta (z_{t+1}) p(x_{t+1} | z_{t+1}, \theta) p(z_{t+1} | z_{t} , \theta) \end{align} $\beta$ の初期値についてはあまり自明ではないですが、 $\beta(z_{T}) = 1$ とおくことになります。これは \begin{align} \alpha (z_{T}) \beta (z_{T}) &= p(z_{T}, \textbf{X})\\ p(z_{T}, \textbf{X}) \beta (z_{T}) &= p(z_{T}, \textbf{X})\\ \beta(z_{T}) &= 1 \end{align} からわかります。

実際の計算の際には、この $\alpha, \beta$ は実際にはとても小さい値になるのでうまくスケーリングし直す必要があります。また今度取り上げます。

E step

データの対数尤度を最大にするような隠れ変数の分布 $q(Z)$ は次で与えられます。（KL divergenceをゼロにするようなものからすぐにわかる） \begin{align} q({\bf Z}) = p( {\bf Z} | {\bf X , \theta}) \end{align} この分布を利用して隠れ変数に関する２つの期待値を計算しておきます。GMMで計算したresponsibilityのようなものです。

時刻tで隠れ変数kである期待値 $\sum_{z_{t}} z_{tk} p( {\bf Z} | {\bf X , \theta}) = \sum_{z_{t}} z_{tk} p(z_{tk} | {\bf X , \theta}) = p(z_{tk} = 1| {\bf X , \theta})$
時刻t-1で隠れ変数jであり、時刻tで隠れ変数kである期待値 $\sum_{z_{t-1j}, z_{t}} z_{t-1j} z_{tk} p(z_{t-1j}, z_{tk} | {\bf X , \theta}) = p(z_{t-1j} = 1, z_{tk} = 1| {\bf X , \theta})$

議論を楽にするために便利な変数 $\gamma, \xi$ を導入します。 \begin{align} \gamma (z_{tk}) &:= p(z_{tk} = 1 | {\bf X , \theta})\\ \xi (z_{t-1j}, z_{tk}) &:= p(z_{t-1j}=1, z_{tk} = 1 | {\bf X , \theta}) \end{align} この変数は上で計算した $\alpha, \beta$ を使って計算できます。 \begin{align} \gamma(z_{tk}) &= \frac{\alpha(z_{tk}=1) \beta (z_{tk} = 1)}{p({\bf X} | \theta)} = \sum_{z_{t}} z_{tk} \frac{\alpha(z_{t}) \beta (z_{t})}{p({\bf X} | \theta)} \\ \xi (z_{t-1j}, z_{tk}) &= \frac{\alpha(z_{t-1j}=1) A_{jk} p(x_{t} | z_{tk} = 1, \phi) \beta (z_{tk} = 1)}{p({\bf X} | \theta)} = \sum_{z_{t-1j}, z_{t}} z_{t-1j} z_{tk} \frac{\alpha(z_{t-1}) p(z_t | z_{t-1}) p(x_{t} | z_{t}, \phi) \beta (z_{t})}{p({\bf X} | \theta)} \end{align}

M step

Q関数（完全データ対数尤度の期待値）を微分してゼロとおいて最大になるパラメーター $\theta$ を求めます。この時、計算した期待値 $\gamma, \xi$ は固定で、E stepで計算した際のパラメータを $\theta_{old}$ とします。 Q関数の具体的な形は \begin{align} Q(\theta_{old}, \theta) &= \sum_{Z} p(Z | X, \theta_{old}) \log p (X, Z | \theta)\\ &= \sum_{Z} p(Z | X, \theta_{old}) \log p (X| Z, \theta) p (Z | \theta)\\ &= \sum_{Z} p(Z | X, \theta_{old}) \log \left( \prod_{t=1}^{T} p (x_{t} | z_{t} , \theta) p(z_{1} | \theta) \prod_{t=1}^{T-1} p(z_{t+1} | z_{t}, \theta) \right) \\ &= \sum_{Z} p(Z | X, \theta_{old}) \log p(z_{1} | \theta) + \sum_{t=1}^{T-1} \sum_{Z} p(Z | X, \theta_{old}) \log p(z_{t+1} | z_{t}, \theta) + \sum_{t=1}^{T} \sum_{Z} p(Z | X, \theta_{old}) \log p (x_{t} | z_{t} , \theta) \\ & = \sum_{z_{1}} p(z_{1} | X, \theta_{old}) \log p(z_{1} | \pi) + \sum_{t=1}^{T-1} \sum_{z_{t+1}, z_{t}} p(z_{t+1}, z_{t} | X, \theta_{old}) \log p(z_{t+1} | z_{t}, \textbf{A}) + \sum_{t=1}^{T} \sum_{z_{t}} p(z_{t} | X, \theta_{old}) \log p (x_{t} | z_{t} , \phi)\\ & = \sum_{k=1}^{K} p(z_{1k}=1 | X, \theta_{old}) \log \pi_{k} + \sum_{t=1}^{T-1} \sum_{j,k=1}^{K} p(z_{tj} =1, z_{t+1k}=1 | X, \theta_{old}) \log A_{jk} + \sum_{t=1}^{T} \sum_{k=1}^{K} p(z_{tk}=1 | X, \theta_{old}) \log p (x_{t} | z_{tk} , \phi)\\ & = \sum_{k=1}^{K} \gamma (z_{1k}) \log \pi_{k} + \sum_{t=1}^{T-1} \sum_{j,k=1}^{K} \xi (z_{tj}, z_{t+1k}) \log A_{jk} + \sum_{t=1}^{T} \sum_{k=1}^{K} \gamma (z_{tk}) \log p (x_{t} | z_{tk} , \phi)\\ \end{align} また、いくつかの制約条件があるのでそれを付け加えてラグランジュ最適化を行います。目的関数 $F$ は以下の形になります。 \begin{align} F(\theta, \alpha, \{\beta_{j} \}) = Q(\theta_{old}, \theta) + \alpha (\sum_{k=1}^{K} \pi_{k} - 1) + \sum_{j=1}^{K} \beta_j (\sum_{k=1}^{K} A_{jk} -1) \end{align} これをそれぞれの変数で微分する。出力確率についてはどのような分布であるかによる。 \begin{align} \frac{\partial F}{\partial \pi_k} = 0, \frac{\partial F}{\partial A_{jk}} = , \frac{\partial F}{\partial \alpha} = 0, \frac{\partial F}{\partial \beta_{j}} = 0 \end{align} 計算すると、 \begin{align} \pi_{k} = \frac{\gamma (z_{1k})}{\sum_{k=1}^{K} \gamma (z_{1k}) }, A_{jk} = \frac{ \sum_{t=1}^{T-1} \xi (z_{tj}, z_{t+1k})}{ \sum_{k=1}^{K} \sum_{t=1}^{T-1} \xi (z_{tj}, z_{t+1k})} \end{align} この結果に対する解釈を与えると、 \begin{align} \pi_{k} = \frac{\textrm{時刻１で状態kにいる確率}}{\textrm{規格化因子}}, A_{jk} = \frac{\textrm{時刻tで状態j, 時刻t+1で状態kにいる確率、の全ての時刻での和}}{\textrm{規格化因子}} \end{align} $A_{jk}$ の分子については、遷移はマルコフ過程で遷移行列が時刻で変わらないとしているので、和をとっていることは自然なことです。

出力確率についてもいくつかの場合で議論しておきます。

観測変数がPoisson分布に従うとき

出力確率は以下のような形になる。 $\phi = \{ \mu_{k} \}_{k=1, \dots, K}$ \begin{align} p (x_{t} | z_{tk} , \mu_{k} ) = \frac{ \exp( - \mu_{k}) \mu_{k}^{x_{k}} }{ x_{k}! } \end{align} 最適化する目的関数の $\{ \mu_{k} \}_{k=1, \dots, K}$ に関するところだけに注目すると、 \begin{align} F(\theta, \alpha, \{\beta_{j} \}) = \sum_{t=1}^{T} \sum_{k=1}^{K} \gamma (z_{tk}) \left( -\mu_{k} + x_{k} \log \mu_{k} \right) + \textrm{const.} \end{align} 微分する。 \begin{align} \frac{\partial F}{\partial \mu_{k} } &= 0\\ \therefore \ \mu_{k} &= \frac{\sum_{t=1}^{T} x_{k} \gamma (z_{tk})}{ \sum_{t=1}^{T} \gamma (z_{tk})} \end{align}

観測変数がGauss分布に従う時

出力確率は以下のような形になる。 $\phi = \{ \mu_{k},\Sigma_{k} \}_{k=1, \dots, K}$ $\displaystyle \begin{align} p (x_{t} | z_{tk} , \mu_{k},\Sigma_{k} ) = {\it N} ( x_{t} | \mu_{k},\Sigma_{k}) = \frac{1}{(\sqrt{2 \pi})^{N} \sqrt{|\Sigma_{k}|}} \exp\left( - \frac{1}{2} \left(x-\mu_{k}\right)^{T} \Sigma^{-1} \left(x-\mu_{k}\right) \right) \end{align}$ 最適化する目的関数の $\{ \mu_{k},\Sigma_{k} \}_{k=1, \dots, K}$ に関するところだけに注目すると、

$\begin{aligned} F(\theta, \alpha, \{\beta_{j} \}) = \sum_{t=1}^{T} \sum_{k=1}^{K} \gamma (z_{tk}) \left( -\frac{1}{2}\log |\Sigma_{k} | - \frac{1}{2} \left(x-\mu_{k}\right)^{T} \Sigma^{-1} \left(x-\mu_{k}\right) \right) + \textrm{const.} \end{aligned}$

微分すると、

$\begin{aligned} \frac{\partial F}{\partial \mu_{k} } &= 0, \frac{\partial F}{\partial \Sigma_{k} } = 0\\ \therefore \ \mu_{k} &= \frac{\sum_{t=1}^{T} x_{k} \gamma (z_{tk})}{ \sum_{t=1}^{T} \gamma (z_{tk})}, \Sigma_{k} = \frac{\sum_{t=1}^{T} \gamma (z_{tk}) \left(x-\mu_{k}\right) \left(x-\mu_{k}\right)^{T}}{ \sum_{t=1}^{T} \gamma (z_{tk})} \end{aligned}$

計算量

フォーワードバックワードアルゴリズムの時間計算量 $O(K^2 T)$ で、空間計算量は $O(K)$ になります。全ての場合をただ計算していくと、時間計算量は $O(K^T)$ になるのでかなり少なくなっています。

EMについては収束がいつ頃起こるわからないのでなんとも言えません。ただ、何千stepも必要になっていたら初期値を選び直した方がいいかもしれません。問題にもよりますが、経験的には数十stepで収束します。

2. 隠れ変数の系列の最大化

viterbiアルゴリズムを紹介します。要するに動的計画法のことです。 $p(\textbf{Z}, \textbf{X} | \theta )$ を最大化します。次のような変数を導入します。 \begin{align} \delta_{t} (z_{t}) = \max_{z_{1}, \dots, z_{t-1}} p(z_{1}, \dots, z_{t}, x_{1}, \dots, x_{t} | \theta ) \end{align} これは時刻t-1までは最適なパスになっている時に、時刻tで $z_{t}$ をとった時の観測変数との同時確率です。引数が行列だと見にくいので、 $\delta_{t} (k) := \delta_{t} (z_{tk}=1)$ と書くことにします。

時刻t+1に進む時を考えると \begin{align} \delta_{t+1} (z_{t+1} ) = [ \max_{z_{t}} \delta_{t} (z_{t}) p(z_{t+1} | z_{t} , \textbf{A}) ] p (x_{t+1} | z_{t+1} ,\phi) \end{align} 成分で見ると \begin{align} \delta_{t+1} (k) = [ \max_{j} \delta_{t} (j) A_{jk} ] p (x_{t+1} | z_{t+1k} ,\phi) \end{align} tで最適な $z_{tk}=1$ を選んだ時の前の時刻の $z_{t-1j}=1$ を記録しておく関数 $\psi_t (k) := j$ を定義する。

初期化 \begin{align} \delta_{1}(k) &= \pi_{k} p (x_{1} | z_{1k},\phi)\\ \psi_{1} &= 0 \end{align}
繰り返し \begin{align} \delta_{t+1} (k) &= [ \max_{j} \delta_{t} (j) A_{jk} ] p (x_{t+1} | z_{t+1k} ,\phi)\\ \psi_{t+1} (k) &= arg\max_{j} \delta_{t} (j) A_{jk} \end{align}
終了 \begin{align} P^{*} &= \max_{k} \delta_{T} (k)\\ k_{T}^{*} &= arg\max_{k} \delta_{T} (k) \end{align}

あとは隠れ変数を $\psi$ で遡れば求めたい系列が得られます。

おまけ: Q関数の最大化がデータ尤度の上昇になることについて

Q関数 $Q(\theta_{old}, \theta) = \sum_{Z} p(Z | X, \theta_{old}) \log p (X, Z | \theta)$ に対して、 \begin{align} Q(\theta_{old}, \theta) - Q(\theta_{old}, \theta_{old}) &= \sum_{Z} p(Z | X, \theta_{old}) \log \frac{p (X, Z | \theta)}{p (X, Z | \theta_{old})}\\ & \le (\sum_{Z} p(Z | X, \theta_{old})) \log \frac{\sum_{Z} p (X, Z | \theta)}{\sum_{Z} p (X, Z | \theta_{old})}\\ & = \log \frac{p (X | \theta)}{p (X | \theta_{old})} \end{align} （途中でJensenの不等式を利用しています。）

よって、 $Q(\theta_{old}, \theta) \ge Q(\theta_{old}, \theta_{old})$ なら、 $p (X | \theta) \ge p (X | \theta_{old})$ . \begin{align} \max_{\theta} Q(\theta_{old}, \theta) \to p (X | \theta) \ge p (X | \theta_{old}) \end{align}