Bayes theoremとHMMとViterbi Algorithm 備忘録

#概要
よく参考にさせていただいている元Google brainのYouTuberの動画がすごくわかりやすかったので、Wikipediaの例に置き換えて書き残す。

https://www.youtube.com/watch?v=kqSzLo9fenk&t=1758s

前半はHMM(Hidden Markov Model)をBayes theorem視点から理解し、後半はその行動(observation)が連なったときに使用するcomputational costの低いViterbi algorithmについて理解する。
また最後にpomegranateの使い方も少し残す。

ここで扱う[Wikipediaの例](https://en.wikipedia.org/wiki/Hidden_Markov_model#A_concrete_example)とは
- Hidden state: rainy, sunny
- Observation: walk, shop, clean
のよくあるHMMである。

またパケージは旧来の[HMMlearn](https://hmmlearn.readthedocs.io/en/latest/index.html)ではなく、より高機能な[pomegranate](https://pomegranate.readthedocs.io/en/latest/)を使用する。

- 実施期間: 2021年12月
- 環境：Ubuntu20.04LTS
- パケージ：pomegranate

##Hidden markov model (HMM)
天気(state)と行動(observation)の遷移について過去統計から確率で図示するところまで説明する。

前提として全stateの遷移は統計的にわかっているものとし、また各stateからemitするobservationも統計的にわかっているものとする。
もちろんここらがわからないといったケースも十分にあり、その場合はforward-backward algorithm(Baum-Welch algorithm)を使用すれば解決可能でありpomegranateではこちらがデフォルトになっている。

過去の統計から天気が下図のように変化することがわかっている。できるだけWikipediaの例になるようにしたので長くなってしまった。
![Screenshot from 2021-12-31 11-23-14.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/9a53ebef-323f-feb3-70be-adcd91eab060.png)
変化は４パターンしかなく、勘定すると下記は容易にわかる。
![Screenshot from 2021-12-31 11-24-47.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/a87d35cf-d960-d131-6b02-fe28c16d1d23.png)

天気、つまりstate(後々hidden state)のtransition probabilityは下図となる。
![Screenshot from 2021-12-31 11-27-34.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/01100767-199c-ed78-8f02-bd052cafeb12.png)

stateの初期状態をprior probability(事前確率)と呼び、次の３式で求めることができる。これで開始stateを確率で表現できる。
![Screenshot from 2021-12-31 12-12-05.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/d973d29a-9254-fc1f-8655-4c3ccc713f5d.png)

ここからはobservationの発生確率(emission probability)を考える。
stateがきっかけとなるobservationも統計的にわかっている。ここでは、walk, shop, cleanの3行動とする。
![Screenshot from 2021-12-31 12-13-37.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/0a75de3b-2a15-46db-411b-0f3373962d30.png)
emission probabilityも上図を勘定するだけ。
![Screenshot from 2021-12-31 12-31-47.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/e8079078-3975-ee56-81b9-a53af6b06059.png)

総合的にまとめるとモデルは下図となる。
![Screenshot from 2021-12-31 12-34-56.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/734fb23d-5913-b488-46dd-f93f6160ddc3.png)


##Bayes theorem視点でHidden stateの推論
上記までで天気->行動を図示したが、ここからは行動変化->天気変化をBayes theoremで推論する。
HMMの本体であり行動が実際にどのように推移したらそのとき天気がどのように推移していたらしいか推論する。stateは不明なのでhidden stateと呼ぶ。
この「推移していたらしい」をposterior probability(事後確率)として求める。

はじめに、observationが変化しないスポットでのhidden stateを推論してみる。
priorはすでにわかっており雨：晴れ＝0.6:0.4なので雨３日、晴れ２日の計５日間の行動はモデルより下図となる。
![Screenshot from 2021-12-31 12-57-03.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/3ad7bd4d-812c-9ab7-7493-d20b7627b8e4.png)

上図より、observationがworkのときのposterior probabilityは
![Screenshot from 2021-12-31 13-03-59.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/61b5ce71-d8d6-1069-95eb-163e1930e0e9.png)

同様にshopとcleanは下記となる。
![Screenshot from 2021-12-31 13-04-10.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/3a1fd932-13f4-a5cb-e94e-1d09f783c2bb.png)

次に２連続変化するobservationから両日のhidden stateを推論してみる。
![Screenshot from 2021-12-31 13-14-09.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/297e744c-3b98-0b67-a874-f84e29a018ce.png)

prior, transition, emissionの各確率を掛けるだけ、つまり下式の積算により最も確率の高いhidden state変化がどれなのかわかる。
![Screenshot from 2021-12-31 13-20-04.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/be289da8-c730-b737-15e9-3fa5a216f05d.png)
月曜火曜の天気の変化は4パターンあるのでモデルは全体で下図の4モデルである。
![Screenshot from 2021-12-31 13-15-40.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/c35117af-bb73-31a5-4c37-ec744ab80473.png)

つまりもっともらしいhidden stateの変化は雨->雨であったことがわかる。
２番目、３番目に大きかった確率を気にする必要がなければこのあと説明するViterbi algorithmでhidden stateの遷移は推論することができる。

##Viterbi algorithmで推論
hidden state数やobservationの変化数が増えると上述の計算方法では爆発してしまう。2xLxn^L(L:Observation長, n:Hidden state数)
そげなときは最大確率のみに注目するViterbi algorithmをつかうっちゃ。
ここでは４日間のobservationからhidden stateを推論してみる。
![Screenshot from 2021-12-31 13-32-08.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/62e9e16d-1a74-8e15-fcaf-1ac02e1f8e79.png)

手順は図示すると簡単で、それぞれの日における各hidden stateのposterior probabilityを計算し、最後にそれら4つの最大値を選ぶだけ。

1,2日目の変化は下図となる。
![Screenshot from 2022-01-03 10-23-42.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/41cbc591-3f25-db2c-13e5-269bc4193dfa.png)

こんな感じに各日、各hidden stateの最大値を書き残していく。
![Screenshot from 2021-12-31 13-49-56.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/1fc0d458-bc1f-8aec-d28d-6bd3b4116a05.png)

![Screenshot from 2021-12-31 13-50-20.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/b3fb43f8-e913-df36-211e-f34b99cdf5b9.png)

最終的に最大確率（青字）を結べばhidden stateの推論はおわり。
![Screenshot from 2021-12-31 13-52-44.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/9a26b927-68b3-bc94-4848-d05baa9f041d.png)

##pomegranateで確認
###Viterbi algorithm
pomegranateの[チュートリアル](https://github.com/jmschrei/pomegranate/blob/master/examples/hmm_rainy_sunny.ipynb)では上記で使用したWikipediaの例が使われている。

HMMのモデルを作成する。

```python
from pomegranate import *
import random
import math

random.seed(0)

model = HiddenMarkovModel(name="Rainy-Sunny")
```

rainyとsunnyのstateを作成する。このときそれぞれのstateからemitされるobservationもそのemission prob.と一緒に辞書型で定義する。

```python
rainy = State(DiscreteDistribution({'walk': 0.1, 'shop': 0.4, 'clean': 0.5}), name='Rainy')
sunny = State(DiscreteDistribution({'walk': 0.6, 'shop': 0.3, 'clean': 0.1}), name='Sunny')
```

はじめのstateをprior prob.（事前確率）を使ってモデルに定義する。

```python
model.add_transition(model.start, rainy, 0.6)
model.add_transition(model.start, sunny, 0.4)
```
transition prob.をモデルに定義する。

```python
model.add_transition(rainy, rainy, 0.7)
model.add_transition(rainy, sunny, 0.3)
model.add_transition(sunny, rainy, 0.4)
model.add_transition(sunny, sunny, 0.6)
```

bake()というコンパイルのようなファイナライズ処理を行いモデルを固める。

```python
model.bake(verbose=True)

# 可読性のため戻り値のmatrixのindexを定義する。
rainy_idx = model.states.index(rainy)   # bake後でないとindex()を認識しないから
sunny_idx = model.states.index(sunny)
```

推論させたい任意のシーケンスを設定する。ここでは手計算した前出のものを使用する。

```python
sequence = [ 'shop', 'clean', 'walk', 'walk']
```

Viterbi algorithmで推論させる。

```python
print(", ".join(state.name for i, state in model.viterbi(sequence)[1]))
```

**Rainy-Sunny-start, Rainy, Rainy, Sunny, Sunny**

手計算通りの結果が得られる。ここではstate.nameで抽出、joinしているが、viterbi()[1]の戻り値は下記であった。

```python
print(model.viterbi(sequence)[1])
```
![Screenshot from 2022-01-01 12-59-14.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/96493080-6841-ede3-8d92-0033c1761b74.png)

###同時確率
同時確率（全確率和）による推定であればシーケンス途中の詳細を知ることができる。例えばforward algorithmの場合は下図のイメージに変わる。
![Screenshot from 2022-01-03 10-30-22.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/52d70644-1c77-4df0-da45-b5af608bb6da.png)

naiveに全確率を計算させると爆発するのでforward, backward, forward-backward algorithmといったdynamic algorithmを使用する。詳細は[Wikipediaにコード付き](https://en.wikipedia.org/wiki/Forward%E2%80%93backward_algorithm)で解説あり参照のこと。

| Algorithm | [Official Doc](https://pomegranate.readthedocs.io/en/latest/HiddenMarkovModel.html)の説明 |
|:-----------|:------------|
| forward      | シーケンスを順方向に進み、各observationが各stateに整列する確率を計算する。全パスの対数確率の合計が出力される。アンダーフローのエラーを防ぐために各行を動的にスケールするため、正規化している。 |
| backward     | シーケンスを逆方向に進み、各observationが各stateに整列する確率を計算する。全パスの対数確率の合計が出力される。 |
| forward-backward  | emission matrixとtransition matrixを返す。emission matrixは全シーケンスが与えられたとき、各stateがそれらを生成する正規化された確率を返し、transition matrixはあるtransitionが使用される回数の期待値を返す。|

実際のチュートリアルではモデルのend stateも定義している。end stateがなければ永遠にHMMをぐるぐる回り続けることになり現実的ではないからだろう。雨で終わるtransition prob.と晴れで終わるそれを各0.1と定義する。そのため上で求めた4つのtransition prob.からそれぞれ0.05差し引いている。
![Screenshot from 2022-01-03 14-55-51.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/b65a8ebb-8443-8bb2-a47e-61b3ae0fd521.png)

```python
model.add_transition( rainy, rainy, 0.65 )
model.add_transition( rainy, sunny, 0.25 )
model.add_transition( sunny, rainy, 0.35 )
model.add_transition( sunny, sunny, 0.55 )
model.add_transition( rainy, model.end, 0.1 )
model.add_transition( sunny, model.end, 0.1 )
```

シーケンスが与えられたときの任意のステップが指定するhidden stateである確率は下記で得られる。

```python
print(math.e**model.forward(sequence))
```
![Screenshot from 2022-01-03 11-22-27.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/66c52736-5d17-1903-148d-8ddd64247985.png)
尚、end stateを設定しなければforwardの戻りは下記のようになり手計算と一致する。
![Screenshot from 2022-01-01 14-27-18.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/5dc01988-02e6-ba68-3dd1-57f4b20c5fb3.png)
行がシーケンスのステップ、列が各hidden stateの確率、最後から2列目は開始確率、最終列が終点確率になっている。従い、例えば最終ステップで晴れになる確率は下記で得られる。

```python
print(math.e**model.forward(sequence)[len(sequence), sunny_idx])
```
0.007304040000000002

以降、わかりやすくするためにobservationのシーケンスを長くして確認する。Viterbiの結果も併記する。

```python
sequence = [ 'walk', 'shop', 'clean', 'clean', 'walk', 'shop', 'walk', 'shop', 'clean', 'shop', 'clean', 'shop', 'clean', 'shop', 'walk', 'shop' ]
print(", ".join(state.name for i, state in model.viterbi(sequence)[1]))
```
Rainy-Sunny-start, Sunny, Rainy, Rainy, Rainy, Sunny, Sunny, Sunny, Rainy, Rainy, Rainy, Rainy, Rainy, Rainy, Rainy, Sunny, Sunny, Rainy-Sunny-end

forward-backward algorithmが出力するposterior prob.を確認する。

```python
trans, ems = model.forward_backward(sequence)
# transition matrix
print(math.e**trans)
# emission matrix
print(math.e**ems)    # model.predict_proba(sequence)の戻りと同じ
```
![Screenshot from 2022-01-03 12-01-42.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/54cd3cee-742e-71fd-190a-d129b760a6f0.png)
はじめのtransition matrixは下図のような正方行列になっていて、algorithmが推論したstate A -> state Bへの遷移回数が出力されている。
![Screenshot from 2022-01-03 12-03-54.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/275a5880-537a-66c5-5d9e-77fc652ed8d4.png)
実際にVitabiの結果を数えてみると下図となる。
![Screenshot from 2022-01-03 12-08-00.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/2bffcc34-613a-dafb-b151-2781e8433d2d.png)
Viterbiは最後のhidden stateはsunnyだが、forwardではややrainyが大きくなっている。
それ以外の大小関係はあっている。3

次のemission matrixは同様に1列目がrainy, 2列目がsunny, 各行が入力observationに対する確率を示している。まとめると下図となり、最終行のend stateへ渡るhidden stateの推論がviterbiと異なっている以外は同じ結果となっている。
![Screenshot from 2022-01-03 12-28-42.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/14473fa7-5d11-5926-2ce7-f29ece734375.png)
離散値で出力するviterbiよりも、確率で出力してくれるforward-backwardの方が使い勝手が良いだろう。グラフが大きくなったときのcomputational costが気になるところ。

backward algorithmの出力も見ておく。

```python
print(math.e**model.backward(sequence))
```

![Screenshot from 2022-01-03 14-38-11.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/9b7ed6f0-47eb-4933-da33-d489c6432db8.png)

これだけがよくわからなかった。まず全般を通してDocumentには出力されるmatrixをどのように読むのか説明が一切ない。このbackwardは特に情報が少ないので予想になるが、行はシーケンスの終端から遡るので1行目がシーケンスの最後で最終行がシーケンスの先頭なことは間違いない。1,2列目はrainy, sunnyで間違いなく、3列目がend stateで4列目がstart stateになっているようにしか見えない。チュートリアルには各ステップの任意のhidden stateの確率が出力されているとなっているが、どのようにsumをとっても1にならず、本当に確率なのか不明。forward-backwardのためのbackwardなので、単独で使用する機会はないように思う。
forward-backwardが使えれば最強のように思えるが、今後使っていって何かわかったら更新することとする。

最後に最大の事後確率をつないだhidden stateの遷移を出力、というか表示する。確率はすでにforward-backwardで計算済み。

```python
print(" ".join(state.name for i, state in model.maximum_a_posteriori(sequence)[1]))
```
**Sunny Rainy Rainy Rainy Sunny Sunny Sunny Rainy Rainy Rainy Rainy Rainy Rainy Rainy Sunny Rainy**

以上