はじめに
python で HMM のパラメータ推定を実装しました．
教科書として『続・わかりやすいパターン認識』を使いました．
本記事の構成
- Hidden Markov Model
- サイコロ投げ
- 記号の整理
- 評価
- forward algorithm
- backward algorithm
- 推定
- Baum-Welch algorithm
- スケーリング
- 実装
- 結果
- おわりに
Hidden Markov Model
一時点前の状態に依存して現在の状態が確率的に決まるような特性を マルコフ性 と呼び，
マルコフ性を満たす過程を マルコフ過程 と呼びます．
このような事象の例として，株価や音声信号，言語などが挙げられます．
マルコフモデルにおいて，出力記号の系列のみ観測でき，状態の系列が観測できないモデルを
隠れマルコフモデル(Hidden Markov Model) と呼びます．
サイコロ投げ
HMM について，サイコロ投げを例に説明します．
$m$ 種の目 $v_1, v_2, \cdots, v_m$ を持つ $c$ 種のサイコロ $\omega_1, \omega_2, \cdots, \omega_c$ を $n$ 回投げることを考えます，
$t$ 回目に投げたサイコロの種類を $s_t$，そのサイコロを投げて観測された目を $x_t$ で表します．
$s_t$ および $x_t$ は以下のように表されます．
\begin{eqnarray}
\left\{\begin{array}{ll}
s_t \in \bigl\{\omega_1, \omega_2, \cdots, \omega_c\bigr\} \\
x_t \in \bigl\{v_1, v_2, \cdots, v_m\bigr\} \\
\end{array} \right.
\end{eqnarray}
$t$ 回目に取り出すサイコロ $s_t$ は，$(t - 1)$ 回目に取り出すサイコロ $s_{t-1}$ によって確率的に決まるとし，
その確率を $a(s_{t - 1}, s_t)$ で表します．
$(t - 1)$ 回目に取り出したサイコロが $\omega_i$ であるとき，$t$ 回目に取り出すサイコロが $\omega_j$ である確率は，
以下のように表されます．
a(\omega_i, \omega_j) = P(s_t = \omega_j \ | \ s_{t - 1} = \omega_i) \ \ \ \ \ (i, j = 1, 2, \cdots, c) \tag{1}
$t$ 回目にサイコロ $s_t$ を投げて $x_t$ が観測される確率を $b(s_t, x_t)$ で表します．
サイコロ $\omega_j$ を投げて $v_k$ が観測される確率は以下のように表されます．
b(\omega_j, v_k) = P(x_t = v_k \ | \ s_t = \omega_j) \ \ \ \ \ (j = 1, 2, \cdots, c) \ \ (k = 1, 2, \cdots, m) \tag{2}
簡単のため，$a(\omega_i, \omega_j)$ を $a_{ij}$，$b(\omega_j, v_k)$ を $b_{jk}$ と略記し，$a_{ij}$ を 遷移確率，$b_{jk}$ を 出力確率 と呼びます．
HMM では，サイコロを投げて 出力される目（= 観測結果）の情報は得られますが，
その目を出力した サイコロの種類（= 状態）の情報は得られません．
記号の整理
マルコフモデルを表現するための記号を以下に整理します．
- $n$：観測回数
- $c$：状態数
- $m$：出力記号の数
- $\omega_i$：$i$ 番目の状態
- $v_k$：$k$ 番目の出力記号
- $s_t \in \{\omega_1, \omega_2, \cdots, \omega_c\}$：時点 $t$ での状態
- $x_t \in \{v_1, v_2, \cdots, v_m\}$：時点 $t$ での観測結果
- $\boldsymbol s$：状態系列
- $\boldsymbol x$：観測記号系列
- $a_{ij}$：状態 $\omega_i$ から状態 $\omega_j$ への遷移確率
- $b_{jk}$：状態 $\omega_j$ で $v_k$ を出力する確率
- $\rho_i$：初期状態 $(t = 1)$ が $\omega_i$ である確率
- $\boldsymbol A$：$a_{ij}$ を $(i, j)$ 成分に持つ $c \times c$ の行列
- $\boldsymbol B$：$b_{jk}$ を $(j, k)$ 成分に持つ $c \times m$ の行列
- $\boldsymbol \rho = (\rho_1, \rho_2, \cdots, \rho_c)$：$\rho_i$ を成分として持つ $c$ 次元ベクトル
評価
パラメータ $\boldsymbol A, \boldsymbol B, \boldsymbol \rho$ が既知であるとき，
観測結果 $\boldsymbol x$ が得られる確率 $P(\boldsymbol x)$ を求めることを 評価 と呼びます．
\begin{align}
P(\boldsymbol x) &= \sum_{\boldsymbol s} P(\boldsymbol x, \boldsymbol s) \tag{3} \\
&= \sum_{s_1} \sum_{s_2} \cdots \sum_{s_n} P(\boldsymbol x, s_1s_2 \cdots s_n) \\
P(\boldsymbol x, \boldsymbol s) &= \prod_{t = 1}^n a(s_{t - 1}, s_t)b(s_t, x_t) \tag{4}
\end{align}
$\boldsymbol x$ および $\boldsymbol s$ の双方が同時に観測できる場合，$P(\boldsymbol x, \boldsymbol s)$ は容易に計算できます．
一方，$\boldsymbol x$ のみが観測できる場合，起こり得る全ての状態 $\boldsymbol s$ に関する周辺化の計算量は膨大になります．
効率的な計算をするアルゴリズムとして forward algorithm，backward algorithm を紹介します．
forward algorithm
観測結果 $\boldsymbol x$ が得られ，かつ $t$ 回目にサイコロ $\omega_i$ を取り出している確率は以下のようになります．
\begin{align}
& P(\boldsymbol x, s_t = \omega_i) \\
&= P(x_1x_2 \cdots x_t, s_t = \omega_i) \cdot P(x_{t + 1}x_{t + 2} \cdots x_n \ | \ s_t = \omega_i) \\
&= \alpha_t(i) \cdot \beta_t(i) \tag{5}
\end{align}
式 $(5)$ の $\alpha_t(i)$ は，$x_1x_2 \cdots x_t$ が得られ，かつ $t$ 回目にサイコロ $\omega_i$ を取り出している確率です．
$\alpha_t(i)$ は，以下のように再帰的に計算することができます．
\begin{align}
& \alpha_t(j) = \Biggl[\sum_{i = 1}^c \alpha_{t - 1}(i) a_{ij} \Biggr]b(\omega_j, x_t) \tag{6} \\
& (t = 2, 3, \cdots, n) \ \ (j = 1, 2, \cdots, c)
\end{align}
ただし，$\alpha_1(i)$ は以下のようになります．
\begin{align}
\alpha_1(i) &= P(x_1, s_1 = \omega_i) \\
&= \rho_i b(\omega_i, x_1) \ \ \ \ \ (i = 1, 2, \cdots, c) \tag{7}
\end{align}
式 $(6)$ において $t = n$ とすることで以下の式を得ます．
\alpha_n(i) = P(x_1x_2 \cdots x_n, s_n = \omega_i) \tag{8}
式 $(8)$ を $s_n$ に関して周辺化することで下式を得ます．
\begin{align}
P(\boldsymbol x) &= P(x_1x_2 \cdots x_n) \\
&= \sum_{i = 1}^c P(x_1x_2 \cdots x_n, s_n = \omega_i) \\
&= \sum_{i = 1}^c \alpha_n(i) \tag{9}
\end{align}
以上より，forward algorithm の処理手順を以下に示します．
< forward algorithm >
- 初期化
$\ \ \ \ \alpha_1(i) = \rho_i b(\omega_i, x_1) \ \ \ \ \ (i = 1, 2, \cdots, c)$ - 再帰的計算
$\ \ \ \ \alpha_t(j) = \Biggl[\sum_{i = 1}^c \alpha_{t - 1}(i) a_{ij} \Biggr]b(\omega_j, x_t)$ - 確率の計算
$\ \ \ \ P(\boldsymbol x) = \sum_{i = 1}^c \alpha_n(i)$
backward algorithm
式 $(5)$ の $\beta_t(i)$ は，$t$ 回目にサイコロ $\omega_i$ を取り出したという条件のもとで，
それ以降の観測結果が $x_{t+1}x_{t+2} \cdots x_n$ となる確率です．
$\beta_t(i)$ は，以下のように再帰的に計算することができます．
\begin{align}
& \beta_t(i) = \sum_{j = 1}^c a_{ij}b(\omega_j, x_{t + 1})\beta_{t + 1}(j) \tag{10} \\
& (t = 1, 2, \cdots, (n - 1)) \ \ (i = 1, 2, \cdots, c)
\end{align}
ただし，$\beta_n(i)$ は以下のようになります．
\beta_n(i) = 1 \ \ \ \ \ (i = 1, 2, \cdots, c) \tag{11}
式 $(10)$ において $t = 1$ とすることで以下の式を得ます．
\beta_1(i) = P(x_2x_3 \cdots x_n \ | \ s_1 = \omega_i) \tag{12}
式 $(12)$ より $P(x_1x_2 \cdots x_n, s_1 = \omega_i)$ は以下のように表せます．
\begin{align}
& P(x_1x_2 \cdots x_n, s_1 = \omega_i) \\
&= P(s_1 = \omega_i)P(x_1 \ | \ s_1 = \omega_i)P(x_2x_3 \cdots x_n \ | \ s_1 = \omega_i) \\
&= \rho_i b(\omega_i, x_1) \beta_1(i) \tag{13}
\end{align}
式 $(13)$ を $s_n$ に関して周辺化することで下式を得ます．
\begin{align}
P(\boldsymbol x) &= P(x_1x_2 \cdots x_n) \\
&= \sum_{i = 1}^c P(x_1x_2 \cdots x_n, s_1 = \omega_i) \\
&= \sum_{i = 1}^c \rho_i b(\omega_i, x_1) \beta_1(i) \tag{14}
\end{align}
以上より，backward algorithm の処理手順を以下に示します．
< backward algorithm >
- 初期化
$\ \ \ \ \beta_n(i) = 1 \ \ \ \ \ (i = 1, 2, \cdots, c)$ - 再帰的計算
$\ \ \ \ \beta_t(i) = \sum_{j = 1}^c a_{ij}b(\omega_j, x_{t + 1})\beta_{t + 1}(j) \ \ \ \ \ (t = 1, 2, \cdots, (n - 1)) \ \ (i = 1, 2, \cdots, c)$ - 確率の計算
$\ \ \ \ P(\boldsymbol x) = \sum_{i = 1}^c \rho_i b(\omega_i, x_1) \beta_1(i)$
推定
パラメータ $\boldsymbol A, \boldsymbol B, \boldsymbol \rho$ が未知であるとき，
観測結果 $\boldsymbol x$ から未知パラメータを求めることを 推定 と呼びます．
推定の解法を導くための準備として $\gamma_t(i)$ および $\xi_t(i, j)$ を以下のように定義します．
\begin{align}
\gamma_t(i) &\overset{\rm{def}}{=} P(s_t = \omega_i \ | \ \boldsymbol x) \tag{15} \\
\xi_t(i) &\overset{\rm{def}}{=} P(s_t = \omega_i, s_{t + 1} = \omega_j \ | \ \boldsymbol x) \tag{16}
\end{align}
また，$\gamma_t(i)$ および $\xi_t(i, j)$ には以下の関係が成り立ちます．
\begin{align}
\sum_{j = 1}^c \xi_t(i, j) &= \sum_{j = 1}^c P(s_t = \omega_i, s_{t + 1} = \omega_j \ | \ \boldsymbol x) \\
&= P(s_t = \omega_i \ | \ \boldsymbol x) \\
&= \gamma_t(i) \tag{17}
\end{align}
ここで，$P(\boldsymbol x, s_t = \omega_i, s_{t + 1} = \omega_j)$ は以下のように表されます．
\begin{align}
&P(\boldsymbol x, s_t = \omega_i, s_{t + 1} = \omega_j) \\
&= P(x_1 \cdots x_t, s_t = \omega_i) \cdot P(s_{t + 1} = \omega_j \ | \ s_t = \omega_i) \cdot P(x_{t + 1} \ | \ s_{t + 1} = \omega_j) \cdot P(x_{t + 2} \cdots x_n \ | \ s_{t + 1} = \omega_j) \\
&= \alpha_t(i)a_{ij}b(\omega_j, x_{t + 1})\beta_{t + 1}(j) \tag{18}
\end{align}
以上より，下式を導くことができます．
\begin{align}
\gamma_t(i) &= P(s_t = \omega_i \ | \ \boldsymbol x) \\
&= P(\boldsymbol x, s_t = \omega_i) \bigl/ P(\boldsymbol x) \\
&= \alpha_t(i)\beta_t(i) \Bigl/ \sum_{i = 1}^c \alpha_n(i) \tag{19} \\
\xi_t(i, j) &= P(s_t = \omega_i, s_{t + 1} = \omega_j \ | \ \boldsymbol x) \\
&= P(\boldsymbol x, s_t = \omega_i, s_{t + 1} = \omega_j) \bigl/ P(\boldsymbol x) \\
&= \alpha_t(i)a_{ij}b(\omega_j, x_{t + 1})\beta_{t + 1}(j) \Bigl/ \sum_{i = 1}^c \alpha_n(i) \tag{20}
\end{align}
Baum-Welch algorithm
上記の結果をもとに，最尤推定を用いてHMMのパラメータを推定します．
今回説明する推定方法を Baum-Welch algorithm と呼びます．
パラメータ $\boldsymbol A, \boldsymbol B, \boldsymbol \rho$ を以下のようにまとめます．
\boldsymbol \theta = (\boldsymbol A, \boldsymbol B, \boldsymbol \rho) \tag{21}
EMアルゴリズムを適用し，以下の $Q$ 関数を最大化します．
\begin{align}
Q(\boldsymbol \theta^0, \boldsymbol \theta) &= \sum_{\boldsymbol s} P(\boldsymbol s \ | \ \boldsymbol x; \boldsymbol \theta^0) \log P(\boldsymbol x, \boldsymbol s; \boldsymbol \theta) \\
&= \sum_{\boldsymbol s} P(\boldsymbol s \ | \ \boldsymbol x; \boldsymbol \theta^0) \bigl\{\log P(\boldsymbol s; \boldsymbol \theta) + \log P(\boldsymbol x \ | \ \boldsymbol s; \boldsymbol \theta)\bigr\} \\
&= \sum_{\boldsymbol s} P(\boldsymbol s \ | \ \boldsymbol x; \boldsymbol \theta^0) \bigl\{\log P(s_1) + \sum_{t = 1}^{n - 1} \log a(s_t, s_{t + 1}) + \sum_{t = 1}^n \log b(s_t, x_t) \bigr\} \tag{22}
\end{align}
ここで，以下のように定義します．
\begin{align}
Q(\boldsymbol \theta^0, \boldsymbol \rho) &\overset{\rm{def}}{=} \sum_{\boldsymbol s} P(\boldsymbol s \ | \ \boldsymbol x; \boldsymbol \theta^0) \log P(s_1) \tag{23} \\
Q(\boldsymbol \theta^0, \boldsymbol A) &\overset{\rm{def}}{=} \sum_{\boldsymbol s} P(\boldsymbol s \ | \ \boldsymbol x; \boldsymbol \theta^0) \sum_{t = 1}^{n - 1} \log a(s_t, s_{t + 1}) \tag{24} \\
Q(\boldsymbol \theta^0, \boldsymbol B) &\overset{\rm{def}}{=} \sum_{\boldsymbol s} P(\boldsymbol s \ | \ \boldsymbol x; \boldsymbol \theta^0) \sum_{t = 1}^n \log b(s_t, x_t) \tag{25}
\end{align}
式 $(23), (24), (25)$ をそれぞれ $\boldsymbol \rho, \boldsymbol A, \boldsymbol B$ について最大化します．
推定したパラメータは以下のようになります．（導出は割愛します）
\begin{align}
\hat{a}_{ij} &= \sum_{t = 1}^{n - 1} \xi_t(i, j) \Bigl/ \sum_{t = 1}^{n - 1} \gamma_t(i) \tag{26} \\
\hat{b}_{jk} &= \sum_{t = 1}^n \delta(x_t, v_k) \gamma_t(i) \Bigl/ \sum_{t = 1}^n \gamma_t(i) \tag{27} \\
\hat{\rho}_i &= \gamma_1(i) \tag{28}
\end{align}
Baum-Welch algorithm の処理手順を以下に示します．
< Baum-Welch algorithm >
- 初期化
$\ \ \ $ パラメータ $a_{ij}, b_{jk}, \rho_i$ に適当な初期値を設定 - 再帰的計算
$\ \ \ $ 式 $(26), (27), (28)$ より $\hat a_{ij}, \hat b_{jk}, \hat \rho_i$ を計算
$\ \ \ $ ただし $\gamma_t(i), \xi_t(i, j)$ を計算するための $\alpha_t(i), \beta_t(i)$ は forward alogorithm, backward alogorithm で計算 - パラメータの更新
$\ \ \ \ a_{ij} \leftarrow \hat a_{ij}, b_{jk} \leftarrow \hat b_{jk}, \rho_i \leftarrow \hat \rho_i$ によりパラメータを更新 - 終了判定
$\ \ \ $ 対数尤度の増分が予め設定した閾値以下ならば収束したとみなし更新終了
スケーリング
Baum-Welch algorithm では，forward algorithm, backward algorithm を使って $\alpha_t(i), \beta_t(i)$ を計算します．
系列が長くなると，アンダーフローが生じ，$\alpha_t(i), \beta_t(i)$ を計算できなくなります．
アンダーフローの対策として スケーリング を紹介します．
新たに $\hat \alpha_t(i), c_t$ を導入します．
\begin{align}
& \hat \alpha_t(i) \overset{\rm{def}}{=} P(s_t = \omega_i \ | \ x_1 \cdots x_t) = \frac{P(x_1 \cdots x_t, s_t = \omega_i)}{P(x1 \cdots x_t)} \tag{29} \\
& c_t = P(x_t \ | \ x_1 \cdots x_{t - 1}) \tag{30}
\end{align}
ここで，$P(x_1 \cdots x_t)$ を以下のように表します．
\begin{align}
P(x_1 \cdots x_t) &= P(x_1) \cdot P(x_2 \ | \ x_1) \cdot \cdots \cdot P(x_{t - 1} \ | \ x_1 \cdots x_{t - 2}) \cdot P(x_t \ | \ x_1 \cdots x_{t - 1}) \\
&= c_1 \cdot c_2 \cdot \cdots \cdot c_{t - 1} \cdot c_t \\
&= \prod_{m = 1}^t c_m \tag{31}
\end{align}
これより，以下の関係が成り立ちます．
\begin{align}
\hat \alpha_t(i) &= \frac{\alpha_t(i)}{P(x_1 \cdots x_t)} \\
&= \frac{\alpha_t(i)}{\prod_{m = 1}^t c_m} \\
&= \frac{\alpha_t(i)}{c_t \cdot \prod_{m = 1}^{t - 1} c_m} \\
&= \frac{\alpha_t(i)}{c_t \cdot P(x_1 \cdots x_{t - 1})} \tag{32}
\end{align}
式 $(6)$ は以下のように書き換えられます．
\begin{align}
\alpha_t(j) &= \Biggl[\sum_{i = 1}^c \alpha_{t - 1}(i) a_{ij} \Biggr]b(\omega_j, x_t) \\
&= \Biggl[\sum_{i = 1}^c \hat\alpha_{t - 1}(i) P(x_1 \cdots x_{t - 1}) a_{ij} \Biggr]b(\omega_j, x_t) \\
\frac{\alpha_t(j)}{P(x_1 \cdots x_{t - 1})} &= \Biggl[\sum_{i = 1}^c \hat\alpha_{t - 1}(i) a_{ij} \Biggr]b(\omega_j, x_t) \\
c_t \hat\alpha_t(j)&= \Biggl[\sum_{i = 1}^c \hat\alpha_{t - 1}(i) a_{ij} \Biggr]b(\omega_j, x_t) \tag{33}
\end{align}
式 $(33)$ より $c_t$ は以下のように計算できます．
\begin{align}
\sum_{j = 1}^c c_t \hat\alpha_t(j) &= \sum_{j = 1}^c \Biggl[\sum_{i = 1}^c \hat\alpha_{t - 1}(i) a_{ij} \Biggr]b(\omega_j, x_t) \\
c_t &= \sum_{j = 1}^c \Biggl[\sum_{i = 1}^c \hat\alpha_{t - 1}(i) a_{ij} \Biggr]b(\omega_j, x_t) \tag{34}
\end{align}
また，$\hat\beta_t(i)$ は，forward algorithm により求められた $c_t$ を使って，以下のように表されます．
c_{t + 1} \hat\beta_t(i) = \sum_{j = 1}^c a_{ij}b(\omega_j, x_{t + 1})\hat\beta_{t + 1}(j) \tag{35}
$\hat\alpha_t(i), \hat\beta_t(i)$ を持ちいて $\gamma_t(i), \xi_t(i, j)$ は以下のように書き換えられます．
\begin{align}
\gamma_t(i) &= \hat\alpha_t(i)\hat\beta_t(i) \tag{36} \\
\xi_t(i) &= \frac{\hat\alpha_t(i)a_{ij}b(\omega_j, x_{t + 1})\hat\beta_{t + 1}(j)}{c_{t + 1}} \tag{37}
\end{align}
式 $(36), (37)$ を 式 $(26), (27), (28)$ に代入することで，パラメータを更新します．
実装
python で Baum-Welch algorithm を使ったパラメータ推定を実装しました．
実装コードは ここ にあります．
サイコロ投げを題材とし，状態数を $\omega_1, \omega_2, \omega_3$ の $3$ 種（$c = 3$），
出力記号を $v_1, v_2$ の $2$ 種（$m = 2$）としています．
これは，$3$ 種の異なるサイコロを選んでは投げ，出た目の奇数($v_1$)，偶数($v_2$)を観測することに相当します．
結果
真のパラメータを以下に示します．
\boldsymbol A = \begin{pmatrix}
0.1 & 0.7 & 0.2 \\
0.2 & 0.1 & 0.7 \\
0.7 & 0.2 & 0.1
\end{pmatrix}, \ \ \
\boldsymbol B = \begin{pmatrix}
0.9 & 0.1 \\
0.6 & 0.4 \\
0.1 & 0.9
\end{pmatrix}, \ \ \
\boldsymbol \rho = \begin{pmatrix}
1.0 & 0.0 & 0.0
\end{pmatrix}
推定のための初期値は，以下のように設定しました．
\boldsymbol A = \begin{pmatrix}
0.15 & 0.60 & 0.25 \\
0.25 & 0.15 & 0.60 \\
0.60 & 0.25 & 0.15
\end{pmatrix}, \ \ \
\boldsymbol B = \begin{pmatrix}
0.5 & 0.5 \\
0.5 & 0.5 \\
0.5 & 0.5
\end{pmatrix}, \ \ \
\boldsymbol \rho = \begin{pmatrix}
1.0 & 0.0 & 0.0
\end{pmatrix}
推定結果は，以下のようになりました．
\boldsymbol A = \begin{pmatrix}
0.11918790 & 0.63863200 & 0.24218010 \\
0.13359797 & 0.07275209 & 0.79364994 \\
0.72917405 & 0.18715812 & 0.08366783
\end{pmatrix}, \ \ \
\boldsymbol B = \begin{pmatrix}
0.92503797 & 0.07496203 \\
0.56298958 & 0.43701042 \\
0.14059588 & 0.85940412
\end{pmatrix}
また，対数尤度の増大の過程は以下のようになりました．
最後に，$b_{11}, b_{21}, b_{31}$ の収束過程を以下に示します．
細い黒線は，真の値を表しています．
おわりに
python で HMM のパラメータ推定を実装することができました．
間違い等あれば，ご指摘お願いします．