In [1]:
from numpy import zeros
from functools import reduce
$ N $ を状態の数， $ t $ を時刻とする．
そして $ q_t $ を時刻 $ t $ のときの状態とする．
\begin{align} & S &=& \{ S_1, S_2, \dots, S_N \} \\ & q_t \in S \end{align}すると，本来であれば時刻 $ t $ のときに状態 $ q_t $ である確率は，それまでの時刻の情報を用いて
\begin{align} P(q_t = S_j| q_{t-1} = S_i, q_{t-2} = S_k, \dots) \end{align}と表すことができる．
（単純？）マルコフ連鎖とは， $ q_t $ の状態は $ q_{t-1} $ ，すなわち1時刻前の状態にのみ依存して決定されるという仮定をおいたモデルであり，先ほどの確率を以下のように近似する．
\begin{align} P(q_t = S_j| q_{t-1} = S_i, q_{t-2} = S_k, \dots) \approx P(q_t = S_j| q_{t-1} = S_i) \end{align}確率 $ P(q_t = S_j| q_{t-1} = S_i) $ を状態遷移確率と呼ぶ．
これは状態数 $ n $ を次元数とした正方行列を使って表すことができる．
状態遷移確率を表す行列を $ A $ として，以下のように定義する（パラメータは論文に載ってる通り）．
In [77]:
A = zeros((3, 3))
A[0, 0] = 0.4
A[0, 1] = 0.3
A[0, 2] = 0.3
A[1, 0] = 0.2
A[1, 1] = 0.6
A[1, 2] = 0.2
A[2, 0] = 0.1
A[2, 1] = 0.1
A[2, 2] = 0.8
# this is not probability (observed)
pi = [1, 1, 1]
print(A)
天気と天気の添字（天気の添字ってなんやねん）を与える関数
In [47]:
S_label = ['rainy', 'cloudy', 'sunny']
def S(label):
return S_label.index(label)
ある系列が与えられたとき，その系列が実際に生成される確率（モデルのパラメータ（遷移確率）は先程定義した $A$ を使う）
\begin{align} O = \{S_3, S_3, S_3, S_1, S_1, S_3, S_2, S_3\} \end{align}
In [72]:
O = list(map(lambda x: x-1, [3, 3, 3, 1, 1, 3, 2, 3])) # convert from 1-origin to 0-origin
確率は以下のようになる．
reduceは総積を計算している．
In [73]:
probability = pi[O[0]] * reduce(lambda x, y: x * y, [A[O[i], O[i+1]] for i in range(len(O)-1)])
print(probability)
同じ天気が続くような系列を考える．
ある天気について，それが連続するときの連続日数に関する期待値は，
\begin{align} \bar{d} &=& \sum_{d=1}^{\infty} d p_i(d) \\ &=& \sum_{d=1}^{\infty} d a_{ii}^{d-1} (1-a_{ii}) \end{align}で与えることが出来て，近似的に計算すると結果は以下の通り．
In [108]:
for i in range(len(A)):
d_bar = sum(d * pow(A[i, i], d-1) * (1 - A[i, i]) for d in range(1, 100000))
print('{weather}: {d_bar}'.format(weather=S_label[i], d_bar=d_bar))
また，無限級数の公式を使って厳密な値を求めると以下．
\begin{align} \sum_{d=1}^{\infty} d a_{ii}^{d-1} (1-a_{ii}) = \frac{1}{1-a_{ii}} \end{align}結果が一致した．
In [109]:
for i in range(len(A)):
d_bar = 1 / (1 - A[i, i])
print('{weather}: {d_bar}'.format(weather=S_label[i], d_bar=d_bar))
ここまでは，状態は観測された事象に対応していた．
pi = [1, 1, 1]
このモデルは非常に制約が強い．
そこで，状態については，物理的に観測されておらず，確率関数を返すようにする．
この物理的に観測されていない状態に対応できるよう，モデルを拡張する．
隠れ状態は，実際に観測されている状態の系列から推定する．
時刻 $ T $ までに観測された系列を表 $ \mathcal{H} $ と裏 $ \mathcal{T} $を用いて表すと以下のようになる．
\begin{align} O &=& O_1O_2O_3 \dots O_T \\ &=& \mathcal{HHTTTHTTH \dots H} \end{align}二つの部屋と，その部屋を仕切る壁を考える．片方の部屋には，一人の人間とコインがあり，もう片方の部屋には，一人の人間がいる．
コインがある部屋にいる人間は，コインを実際に投げ，「表」「裏」の結果を何らかの方法で別の部屋にいる人間に伝える．
この際，コインが無い方の部屋にいる人間は，コインが「表」「裏」のどちらだったかのみを知る．
この結果として，コインが無い方の部屋にいる人間は，コインの「表」「裏」の系列を知ることができ，コインがどのように投げられたか，コインがどのようなものであるか（コインが「表」「裏」を出す確率など）を知ることは出来ない．
このような状況において，コインの結果系列を説明するためのモデルを作ることを考える．
コインは1つだけだと仮定する．すると，コインの面に対応する2つの状態を持つモデルを考えることになる．このときは状態が観測可能なので，推定するべきなのはコインが「表」「裏」になる確率（bias）である．
コインは同様に2つあると仮定する．加えて，１回前のコイントスに使ったコインの種類によって次に使うコインの種類が決定すると仮定する．
コインが３つあると仮定する．加えて，１回前のコイントスに使ったコインの種類によって次に使うコインの種類が決定すると仮定する．
壺の例は割愛．本文読んでわからなかったら続わかりやすいパターン認識とか見ると良いかもしれません．
HMMのモデルに必要なものは何かを考える．
$ N \in \mathbb{N} $ であり，状態の数を表す．
$ M \in \mathbb{N} $ であり，観測系列の長さを表す．例えば，10回コイントスをするときは $ M = 10 $
$ N $ 次の正方行列であり，状態の遷移確率を表す．
各要素は $ a_{i j} \geq 0 $ ，および $ \sum_j a_{i j} = 1 $ を満たす．
状態 $ j $ のときの出力の確率を表す行列である．
例えば，上のコイン2つのモデルを考えると， $ 2 * 2 $ の行列となる（状態の数（コインの数 $ N $ ） * 出力ラベルの数（「表」，「裏」））．
$ A $ は1時刻前の状態に応じて現在の状態を決定するための確率を与えるが，それでは $ t=0 $ のときはどのように状態を決定すればよいのか？
それを決定するためのものが $ \pi $ であり，初期状態を制御するパラメータである．N次元のベクトルである．
また， $ \sum_i^N \pi_i = 1 $ を満たす．
値の推定が必要になるのは $ A $，$ B $，$ \pi $であり，これらをまとめたものを $ \lambda = (A, B, \pi) $ と書いておくことにする．
HMMのモデル，及びパラメータについて議論した．
実際にHMMを使うときに解かなければならない問題はなんだろうか？
問題を列挙する．なお，系列長を $ T $，観測した事象の結果を $ O = O_1 O_2 \dots O_T $ とする．
これは評価の問題であり，モデルのパラメータと観測結果が与えられている．
「モデルがどれくらい観測結果に適合しているか」をスコアリングする問題であると解釈する．
具体的に音声認識システムを考える．
$ W $ 種類の語彙が出現する言語（出現する単語列は $ L(R) $ に含まれるとする）を仮定する．
このとき， $ N $ 状態のHMMを考える．
$ P(O|\lambda) $ を計算する方法について考える．ここでも $ T $ および $ O = O_1 O_2 \dots O_T $ は与えられているものとする．
簡単に思いつく方法として，すべてのありうる状態系列についてのシュミレーションを行ってしまうことを考えよう．
状態の系列を $ Q = q_1 q_2 \dots q_t $ とする $ q_1 $ は初期状態であり $ \pi $によって決定する．
他の状態は状態遷移によって決定する．
観測結果 $ O $ が生成される確率は， $ P(O|Q, \lambda) $ で表される．
独立性を仮定しているので，以下のように展開できる．
\begin{align} P(O|Q, \lambda) &=& \prod_{t=1}^T P(O_t|q_t, \lambda) \\ &=& b_{q_1}(O_1) b_{q_2}(O_2) \dots b_{q_T}(O_T) \end{align}次に，状態系列 $ Q $ が生成される確率を考える．
\begin{align} P(Q|\lambda) = \pi_{q_1} a_{q_1 q_2} a_{q_2 q_3} \dots a_{q_{T-1} q_T} \end{align}では， $ O $ と $ Q $ の同時確率 $ P(O, Q|\lambda) $ は？
条件付き確率の定義にしたがって分解してみると，以下のようになる．
\begin{align} P(O, Q|\lambda) = P(O|Q, \lambda) P(Q|\lambda) \end{align}これで，先ほど導いた二つの確率を代入することが出来そうだ．
\begin{align} P(O, Q|\lambda) &=& P(O|Q, \lambda) P(Q|\lambda) \\ &=& b_{q_1}(O_1) b_{q_2}(O_2) \dots b_{q_T}(O_T) \pi_{q_1} a_{q_1 q_2} a_{q_2 q_3} \dots a_{q_{T-1} q_T} \end{align}求めたいのは， $ P(O|\lambda) $ である．どうすればよいか，周辺化すれば良い．周辺化及び整理をしてやると以下のようになる．
\begin{align} P(O | \lambda) &=& \sum_Q P(O, Q|\lambda) \\ &=& \sum_Q b_{q_1}(O_1) b_{q_2}(O_2) \dots b_{q_T}(O_T) \pi_{q_1} a_{q_1 q_2} a_{q_2 q_3} \dots a_{q_{T-1} q_T} \\ &=& \sum_Q \pi_{q_1} b_{q_1}(O_1) a_{q_1 q_2} b_{q_2}(O_2) a_{q_2 q_3} \dots a_{q_{T-1} q_T} b_{q_T}(O_T) \end{align}先ほどの方針にしたがって，すべてのありうる状態系列についてのシュミレーションを行ってしまうことを考える．
すると，計算量が $2T \cdot N^T $ であることがわかる．各状態ごとに $ T $ 個の系列が候補として存在する．
$ N $と$ T $が大きいと計算はまったく不可能になってしまう．どうすればよいか？
この計算は，前向き-後ろ向きアルゴリズムと呼ばれるアルゴリズムによって効率的に行える．
$ P(O| \lambda) $ を計算するため，以下のような変数（forward-variable）を定義する．
\begin{align} \alpha_t(i) = P(O_1, O_2, \dots, q_t = S_i | \lambda) \\ \end{align}この変数を， $ t = 1 $ から $ T $ まで逐次計算していくことで，効率的に $ P(O | \lambda) $を計算するのである．
この変数は，最終的に $ T \cdot N $ の行列になる．
まず，初期化処理を行う（Initialize）．
\begin{align} \alpha_1(i) = \pi_i b_i(O_1) (1 \leq i \leq N) \end{align}次に，逐次計算していく（Induction）．
\begin{align} \alpha_{t+1}(i) = \pi_i b_i(O_t) && (1 \leq j\leq N,\ 1 \leq t \leq T-1) \end{align}計算が終了したら， $ P(O | \lambda) = \sum_{i=1}^N \alpha_T(i) $ によって $ P(O | \lambda) $ を得ることができる．
本来であれば，時刻 $ t $ における $ \alpha_{t}(i) $ の計算は，時刻 $ 1 $ から時刻 $ t-1 $ までの $ \alpha_k(i) $ を計算する必要がある．
しかし，前向きアルゴリズムでは時刻 $ t-1 $ における $ \alpha_{t-1}(i) $ の計算結果を保存しておき， 時刻 $ t $ における $\alpha_{t}(i) $ の計算に利用することで繰り返し計算している処理を排除している．
これは動的計画法と呼ばれるアイデアである． $ O(N^2 \cdot T) $
同様に，以下のような変数（backward-variable）を定義する．
\begin{align} \beta_t(i) = P(O_{t+1} O_{t+2} \dots O_{T} |q_t = S_i, \lambda) \end{align}初期化処理
\begin{align} \beta_T(i) = 1 && (1 \leq i \leq N) \end{align}逐次計算
\begin{align} \beta_T(i) = \sum_{j=1}^N a_{i j} b_j(O_{t+1}) \beta_{t+1}(j) && (t = T-1, T-2, \dots, 1, 1 \leq i \leq N) \end{align}Problem2は，Problem 1と異なり，「何を持って最適とするか」を考える難しさが存在する．
まず単純に，系列の中の各状態について，最も出現しやすい出力ラベルを与えることが考えられる．
これは， $ \gamma_t(i) = \frac{\alpha_t(i) \beta_t(i)}{P(O | \lambda)} $ として， $ argmax_{1 \leq i \leq N} \gamma_t(i) $ を解くことで実現できる．
この方法では，遷移の仕方としてありえないような状態列を予測してしまうことがあり，適切ではない場合がある．
ではどうすればよりよい状態系列の予測が行えるか？
もっとも広く用いられる手法は， $ P(Q|O, \lambda) $ を最大化することで，状態系列の最適化を行うことである．
これはビタビアルゴリズムによって実現できる．
\begin{align} Q &=& \{ q_1, q_2, \dots, q_T \} \\ O &=& \{ o_1, o_2, \dots, o_T \} \\ \delta_t(i) &=& \max_{q_1, q_2, \dots, q_t-1} P(q_1, q_2, \dots, q_t = i, O_1, O_2, \dots, O_t| \lambda) \end{align}このとき， $ \delta_t(i) $ は時刻 $ t$ における状態 $ S_i $ への最短経路のスコアである．時刻をずらすと以下のようになる．
\begin{align} \delta_{t+1}(j) = [max_i \delta_t(i) a_{i j}] \cdot b_j(O_{t+1}) \end{align}このように， $ t $ を増加させ逐次計算をしていく．具体的な手続きを記述する．
初期化
\begin{align} \delta_1(1) &=& \pi_i b_i(O_1) && (1 \leq i \leq N) \\ \psi_1(i) &=& 0 \end{align}再帰計算
\begin{align} \delta_t(j) &=& max_{1 \leq i \leq N}[\delta_{t-1}(i) a_{i j} ] b_j(O_t) && (2 \leq t \leq T,\ 1 \leq j \leq N) \\ \psi_t(j) &=& argmax_{1 \leq i \leq N} [\delta_{t-1}(i) a_{i j}] && (2 \leq t \leq T,\ 1 \leq j \leq N) \end{align}計算終了
\begin{align} p^* = max_{1 \leq i \leq N} \delta_T(i) \\ q^*_T = argmax_{1 \leq i \leq N} \delta_T(i) \\ \end{align}最適パス（系列）のトラッキング
\begin{align} q^*_t = \psi_{t+1}(q^*_{t+1}) && (t = T-1, T-2, \dots, 1) \end{align}よく見ると（よく見なくても？），ビタビアルゴリズムって上で書いた前向きアルゴリズムに似てますよね．
ここまで，HMMのモデル（パラメータ）は予め得たものとして話を進めてきた．
しかし，実際の問題に取り組むとき，パラメータは予め存在しない．
したがって，HMMの学習をする必要がある（モデルがなくては何も出来ない）．
モデルの学習とは，パラメータ $ \lambda = (A, B, \pi) $ を最適化することである．
ここで言う最適化とは，観測結果を生成する確率が最大になるようにパラメータを調整することを指す．
この最適化は，解析的に行う手法はまだ知られておらず，バウム=ウェルチアルゴリズムや勾配法などの反復計算が必要な手法によって行われる．
（バウム=ウェルチアルゴリズムとは，EMアルゴリズムの一種です．HMMの学習を行うEMアルゴリズムのことをバウム=ウェルチアルゴリズムと呼ぶみたい？： 参考）
それでは，実際にパラメータの最適化に取り組みましょう．
まず，時刻 $ t $ に状態 $ S_i $ に存在し，時刻 $ t+1 $ に状態 $ S_j $ に存在する確率を表す変数を定義します．
さらに，条件付き確率の定義，前向き・後ろ向き変数の定義にしたがって分解することができ，以下のようになります．
（条件付き確率を導入するための発想ができないけど，ベイズ界隈がよくやる式変形とかなんですかね）
\begin{align} \xi_t(i, j) &=& P(q_t=S_i, q_{t+1}=S_j | O, \lambda) \\ &=& \frac{P(q_t = S_i, q_{t+1} = S_j, O | \lambda)}{P(O | \lambda)} \\ &=& \frac{\alpha_t(i) a_{i j} b_j(O_{t+1}) \beta_{t+1}(j)}{P(O | \lambda)} \\ &=& \frac{\alpha_t(i) a_{i j} b_j(O_{t+1}) \beta_{t+1}(j)}{ \sum_i^N \sum_j^N \alpha_t(i) a_{i j} b_j(O_{t+1}) \beta_{t+1}(j) } \end{align}気づく方もおられるかもしれないが，この $ \xi_t(i, j) $ を用いて先ほど定義した $ \gamma_t(i) $ を表すことができる．
$ \xi_t(i, j) $ $ j $ について周辺化してやると，
\begin{align} \gamma_t(i) = \sum_{j=1}^N \xi_t(i, j) \end{align}となる．
ちなみに，ここまでがEMアルゴリズムのEステップである．
いま， $ \gamma_t(i) $ および $ \xi_t(i, j) $には時刻 $ t $ という変数が含まれている．
最終的に，遷移確率や出力確率は時刻に依存しないはずなので，この時刻に関する変数を消去することを考えたい．
$ 1 \leq k \leq T-1 $ の範囲で $ \gamma_t(i) $ の総和を取ると，これは状態 $ S_i $ から遷移してきた回数の期待値になる．
$ 1 \leq k \leq T-1 $ の範囲で $ \xi_t(i, j) $ の総和を取ると，これは状態 $ S_i $ から状態 $ S_j $ へと遷移する回数の期待値になる．
それぞれは以下のように書ける．
\begin{align} \sum_{t=1}^{T-1} \gamma_t(i) && \mbox{ 状態 $ S_i $ から遷移してきた回数の期待値}\\ \sum_{t=1}^{T-1} \xi_t(i, j) && \mbox{ 状態 $ S_i $ から状態 $ S_j $ へと遷移する回数の期待値}\\ \end{align}上記の式を用いて， $ \bar{\lambda} = (\bar{\pi}, \bar{A}, \bar{B}) $ を求める式を与えたい．以下のようになる（最尤）．
\begin{align} \bar{\pi_i} &=& \mbox{時刻 $ t=1 $ に状態 $ S_i $ にいる回数の期待値} &=& \gamma_1(i) \\ \bar{a_{i j}} &=& \frac{\mbox{状態 $ S_i $ から状態 $ S_j $ へと遷移する回数の期待値}}{\mbox{状態 $ S_i $ から遷移してくる回数の期待値}} &=& \frac{\sum_{t=1}^{T-1} \xi_t(i, j)}{ \sum_{t=1}^{T-1} \gamma_t(i) } \\ \bar{b_j(k)} &=& \frac{\mbox{状態 $ S_j $ において出力 $ v_k $ が得られる回数の期待値}}{\mbox{状態 $ S_j $ に遷移する回数の期待値}} &=& \frac{\sum_{t=1, s.t. O_t = v_k}^{T} \gamma_t(j)}{ \sum_{t=1}^{T} \gamma_t(i) } \\ \end{align}現在のモデルを $ \lambda = (\pi, A, B) $ として，上の式をそれぞれ計算することで $ \bar{\lambda} $を得る．
このとき，つねに $ P(O | \bar{\lambda}) \geq P(O | \lambda) $ が成り立っている（パラメータは常に良い方向に更新される）．
ただし，これは大域的な最適解とは限らない．
さきほどの $ \lambda $ を得るためには，以下に示す補助関数 $ Q(\lambda | \bar{\lambda}) $ を変数 $ \bar{\lambda} $ について最大化してやれば良い（続パタのEMアルゴリズムの章とかを見たらもう少し詳細にわかりそう）．
\begin{align} Q(\lambda | \bar{\lambda}) = \sum_Q P(Q|O, \lambda) log P(O, Q| \bar{\lambda}) \end{align}\begin{align} & max_{\bar{\lambda}} Q(\lambda, | \bar{\lambda}) \Rightarrow P(O| \bar{\lambda}) \geq P(O | \lambda) & \\ & \therefore \bar{\lambda} = argmin_{\bar{\lambda}} Q(\lambda, | \bar{\lambda}) & \end{align}あと一息，ここで $ \pi, A, B $ に関する性質について考えよう．確率の定義から，
\begin{align} \sum_{i=1}^N \pi_i &=& 1 \\ \sum_{j=1}^N a_{i j} &=& 1 && (1 \leq i \leq N) \\ \sum_{k=1}^M b_j(k) &=& 1 && (1 \leq j \leq N) \\ \end{align}が言える．
あとは，先ほど定義したQ関数を，今得た３つの等式を制約式としてラグランジュの未定定数法を用いることで， $ \lambda $ を得ることができる．
\begin{align} \pi_i &=& \frac{\pi_i \frac{\partial{P}}{\partial{\pi_i}}}{\sum_{k=1}^N \pi_k \frac{\partial{P}}{\partial{\pi_k}} } \\ a_{i j} &=& \frac{a_{i j} \frac{\partial{P}}{\partial{a_{i j}}}}{\sum_{k=1}^N a_{i k} \frac{\partial{P}}{\partial{a_{i k}}} } \\ b_j(k) &=& \frac{b_j(k) \frac{\partial{P}}{\partial{b_{j}(k)}}}{\sum_{l=1}^N b_j(l) \frac{\partial{P}}{\partial{b_j(l)}} } \\ \end{align}この３つの等式は，上で定義した数え上げを用いた等式に変形できる．
これでパラメータ更新のアルゴリズムが導出できた．