隱馬可夫模型的三大經典問題

Sampling and Monte Carlo (Special)： 3 classic questions of a hidden markov model

隱馬可夫模型 Hidden Markov Model

回顧上面幾篇提到的馬可夫鏈，本質上是對一個連續過程的模擬，並假設每次狀態變化的機率僅與上次的狀態有關。

然而，現實生活中我們可觀測到的馬可夫鏈往往只有輸出狀態，而不知其隱含的轉移機率與當前狀態。舉幾個例子 : Siri 接收了使用者輸入的一段語音以後，反饋一段語音答案，其中使用者輸入的語音會被切割成字段，並編碼成包含初始狀態與轉移機率的一段馬可夫鏈，最後輸出結果，而我們通常只能觀察到輸出結果，無法得知後端是怎麼編碼與轉移狀態的，如果我們能夠用某些方法，從輸出值觀察到後端發生的事情，我們就可以玩壞 Siri …?

當然，現在的語音識別有許多更精確的演算法在背後支援，但背後的概念都與上面提及的概念擺脫不了關係，這樣的算法過程我們稱做隱馬可夫模型，因此理解隱馬可夫模型對後續做自然語言處理，甚至到強化學習的概念，有至關重要的作用。因為另一種層面來說，隱馬可夫模型也可視為一種包含環境改變來選擇輸出策略的架構。

首先，我們先視覺化看一下隱馬可夫模型 :

Source : Hidden Markov Model (Wikipedia)

我們透過數學的語言來描述，在隱馬可夫模型中，有幾個重要的名詞 :

1. x(t) 是隱藏狀態或隱藏變數，基本上觀察者無法直接觀察到該變數，因此這是我們 "假想" 出來的，代表實際上有某種決策因子在影響你的結果。後續為了方便閱讀我們用 h(t) 代替，表示 hidden。
2. y(t) 是觀測狀態或觀測變數，是我們實際觀察到的情況，比方說投硬幣，連續投了五次正正反正正，這五次就是我們觀察到的狀態，而每次拋硬幣時手的力道大小與方向，空氣的風速等就是隱藏狀態。後續為了方便閱讀以 o(t) 代替，表示 observed。

此外，我們對這個模型做了兩個假設 :

一、齊次馬可夫鏈假設 : 也就是任意時刻的隱藏狀態只跟上一時刻的隱藏狀態有關，以剛剛的投硬幣假設來說，就是假設你每次出去的力道，會跟你上次的力道有關，舉個例子來說像是你跟人打賭能投幾次正面，發現上次太小力出現反面時，可能你下意識會加重力道。(Bad Assumption Though)

二、觀測獨立假設 : 觀測狀態只依賴於當前的隱藏狀態，也就是說當有某種機率函數將觀測狀態轉換為輸出狀態，該機率函數不會受其他時刻影響

因此，一個隱馬可夫模型的產出過程如下 :

1. 設定初始隱含狀態並輸出初始觀察狀態
2. 根據隱狀態轉移機率，轉移至下一狀態
3. 根據該狀態的輸出機率函數，輸出當前觀測狀態

這個概念相對簡單，前面的章節也有談到過一些，我就不贅述了。如果想了解如何實現隱馬可夫模型，可以參考台師大的演算法筆記 :

演算法筆記 - Hidden Markov Model

而今天要談及的就是，我們如何在給定其中幾者的情況下，對未知的參數做出最佳的估計，這衍生出了馬可夫的三大經典問題，分別是 :

1. 給定模型參數 lambda = (A, B, I)，其中A為隱藏狀態下轉移的機率矩陣、B為觀測狀態生成機率矩陣、I 則為初始機率分布，當我們有觀察序列 O = (o1, o2, … , ot)，就能評估在該參數組合下，此序列出現的概率。其應用是我們能夠對給定的模型與輸出去評估其機率值，藉以針對幾種可能的發生情況做期望值評估，比如我們事先知道幾種可能的輸出狀況，每一種輸出狀況對應到不同的策略與效益成本損失，因此我們可以評估所有狀況後，得出未來整體的期望損失以及了解如何結合不同策略。
2. 當然承上述，我們也想知道在給定模型參數與觀測序列後，反向推估隱含序列，相比前兩個問題不涉及對隱含序列的推估，這種算法的應用面又更廣泛一些，因為涉及對觀察現象背後成因的推估。比如根據症狀了解病人健康狀態的可能性、根據供貨銷售量推估供應鏈是否充足等。
3. 最後一個問題則是給定觀察序列以後，我們如何設定模型參數來讓該序列發生的條件機率最大。通常這可以用來反向設計模型，讓其達成我們想要的效果，比如我們已知一組投擲硬幣的正反面結果，反向去估計該硬幣出現正反面的機率有多少。

以下，我們就針對這三個問題來描述 :

前向後向算法 Forward-Backward Algorithm

第一個問題是給定參數組合與觀察序列的情況下，推估該序列在此模型中出現的條件機率，當然，以電腦科學的角度來說，這可以被暴力(Bruteforce)推算出來，只要給定初始狀態以及參數組合(轉移機率)，我們可以用這種方式去計算所有可能的狀態，然後給出某一特定序列狀態條件機率的估計值。假設轉移狀態只會有兩種，也就是兩種機率，那麼 :

當然，學過演算法的同學就看的出來，這至少是個指數級別複雜度的演算法，計算複雜度是 O(T*N^T)，如果序列 T 足夠長或是隱藏狀態數 N 足夠多，基本上算個三天三夜你也得不到結果，因此我們得找一個簡潔的辦法。

但在此之前，我們必須稍微複習一下演算法裏頭提過的動態規劃算法 : 實際上動態規劃與其說是一種算法，不如說是一種思維，當我們處理一個大問題的時候，有時候我們會發現解決該問題並不一定要從頭解起，而是先解決該問題前面的一些子問題，將該子問題的狀態儲存起來，就可以省去從頭推算的時間，舉例來說我們可以看下圖 :

實際上要推導出F1的期望機率，並不需要你給定每個初始狀態可能的值並從頭推導起，只要知道F3與F2的期望機率就行，而反之推導就會得到F3來自F4, F5，F4僅來自於F5，換言之F3來自F5以及F5平方的加總機率，順著這樣的思路推導下去，給定初始狀態選取任一值的機率與轉移機率，我們就可以用O(N)的時間推算出期望機率，而不需要將每條路徑走一遍，大大增加了效率，然而，卻會占用掉一部分記憶體空間已儲存各節點的期望機率。

這個跟我們要解決的問題有甚麼關係呢? 首先，我們把剛才暴力解的方法用數學的語言補充一下 ，首先先定義一些參數:

其中，我們要求解的對象為 : P(O, H|λ)，如果依照傳統的想法，我們可以通過 A, B, I 三個參數，並且知道任意一組隱含序列 H，就可以通過分別求出 P(H|λ) 與 P(O|H, λ)，這樣就可以計算聯合分布來獲得 P(O, H|λ)，接著只需要列舉出所有的隱含序列 H，理論上就可以知道所以 O 出現的機率。

當然，這方法聽起來十分蠢，因此我們需要透過剛才提及到的動態規劃方法，通過逆推的方式將局部問題推廣到全域問題上。首先，我們定義在時刻T的時候，如果隱藏狀態為 j，那麼我們可以得到機率 :

這東西很簡單，給定的參數 λ 自然不必多說，這意思是說在給定參數的時刻t下，觀測到隱含狀態為 j，觀測序列為 o1, o2, …, ot 的機率，我們稱這東西叫做是前向機率 (Forward Probabilities)。那麼，只要我們將這個機率乘上狀態轉移的機率，我們就可以得到下一時刻 t+1 時，狀態為 k 的機率，記為 :

那麼此時所有在時刻 t 往下轉移到新隱藏狀態 k 的機率有多少呢? 很簡單，因為馬可夫互相獨立的性質因素，我們可以簡單將其加總得到 :

這東西就是時刻 t 觀察到狀態為 o1, o2, …, ot，且時刻 t+1 為新隱含狀態 k 的機率，並且也由於這個互相獨立的性質，我們可以推廣到時刻 t+1 時，觀察到 o1, o2, …, ot, ot+1，且其隱藏狀態為 k 的機率 :

比照這個模式，我們可以從 t = 1 開始往後推廣，如果要計算在時刻 T 的時候觀測序列為 o1, o2, …, ot，則只要將所有隱藏序列的值加總計算即可 :

整個流程如果不好理解，可以參考 Stanford 自然語言處理課程中放上的示意圖，我將底下的機率從 α 改寫成 P，以免容易與轉移矩陣中的 A 混淆。

接下來我們很快速透過程式碼來實戰看看 :

程式碼我就不一行一行解釋了，應該不難懂，簡單來說我創造了兩個 Class，一個產生隨機機率矩陣，一個是算法計算，使用者要輸入的僅是觀察序列、狀態可能的數量與序列時間長度即可，由 34~45 行可以看出，計算複雜度在 O(TN**2)。以下來說說後向算法 :

後向算法的概念很簡單，與前向幾乎一模一樣，差別只是在於動態規劃的路徑是由末端往前倒推的，這邊先直接說明公式 :

可以看見公式與剛才的前向算法長得十分相似，這是因為整個後向算法的思路只是由後往前推演，還記得剛才的前向算法我們提到，對於當前時刻 t 且觀測序列 o1, ..., ot 時，可透過所有對應隱藏狀態的機率相加得到，對應的後向算法則是先有了 t+1 時間點的隱藏狀態機率後，把所有可能產生的隱藏狀態的機率加總，往回推到時刻 t 且隱藏狀態為特定值時的機率。

這邊就不再逐步推導，請大家自行嘗試；此外，結合兩種算法的思想，我們也可以推算出單一狀態 t 時隱藏狀態為 i 的機率，或是狀態 t 時隱藏狀態 i，且狀態 t+1 時隱藏狀態為 j 的機率 :

隱藏序列在時刻 t 時為 i 的機率

很容易理解，給定觀測序列 O 與參數 λ 下隱藏狀態的時刻 t，其隱藏狀態為 i 的機率可以透過條件機率改寫，並透過機率矩陣相乘的方式得到。

隱藏序列在時刻 t 時為 i，且時刻 t+1 為 j 的機率 (P’ 為後向機率，P為前向)

維特比算法 Viterbi Algorithm

在 Stanford 自然語言課程的講義中有一個很有趣的比喻 : 在一個冰淇淋銷售員的業績中，給定一組觀察到的 (3, 1, 3)，你可以去猜測出最佳的隱藏狀態天氣序列為 (熱天、熱天、熱天)。Viterbi 就是一個解碼的過程，與上面介紹的 Forward-Backward 有點像，我們手上都有參數組合與觀測序列，且同樣是基於動態規劃去求解的演算法，只是不同於專注在觀測序列本身出現的機率，在這裡我們想要反向解碼最佳的可能隱藏序列。

當然，有人可能會好奇說，剛才透過前向後向算法不是已經可以推估出在時刻 t 之下隱藏狀態為 i 的機率了嗎 ? 但一來預測單一條件並不能保證整體是最佳輸出序列，二來就算是最佳情況，即轉移機率都良好的定義(大於0小於1)的情況下，也需要繼續進行迭代讓算法複雜度進一步提升。

因此，我們需要來重新研究一個可以同時考慮整體的算法。在這裡維特比因為是 Viterbi，我統一用 v(t) 來做為符號。首先，我們先從局部著手(動態算法的典型思維)，在時刻 t 隱藏狀態為 j 的機率為 :

上面式子的第二行為 Recursive Expression，讀者可自行演算一下。式子代表甚麼意思呢? 意思是說每一次時刻 t 隱藏狀態為 j 的機率，都是從每一條時刻 t-1，狀態為 j，並且乘上 j 到 i 轉移機率中選取一條最大的路徑，得到該機率後乘上該隱藏狀態轉移到 Ot 的機率得來，這樣就得到了在所有隱藏路徑中得到下一個狀態 j 的最大路徑機率值。(對這邊有點拗口，要慢慢讀)

接下來，這邊要用到一個觀念叫做 Backpointer，其定義是在時刻 t 且轉移狀態為 i 的情況下，那麼從 t = 1 ~ t-1 中所有可能的轉移狀態路徑中，找出最佳的那一個，因此可以寫成 :

注意以上是 argmax，因此只需要返回 index 即可，之所以這樣做的原因是因為，Backpointer 是返回在我們得到上一個時刻的最佳可能路徑時，該路徑最後一時點的隱藏狀態，而非該機率值。

以上兩個都是定義，接下來就要進入到算法部分。由於文字會有點拗口，這邊用圖片呈現以幫助理解 :

Source : Wikipedia

在上面的圖片中說明了一個病人的隱藏狀態(健康或發燒)，以及其觀測狀態(頭暈、發寒與健康)，想像一個病人走進診間說 : 我第一天感到正常，第二天開始發寒，第三天感到頭暈，那我們是如何計算 Viterbi 來說明他這三天的實際健康狀態呢? (請參考實際案例佐演算法步驟服用)

1. 初始化:給定觀測序列 o1 的情況下，得到所有隱藏狀態的可能機率 v1(i)，並將起始 Backpointer 設為 0
2. Recursion:計算時刻 t = 2, 3, ..., T-1 的所有最佳可能路徑機率值，並且對每一時刻計算 Backpointer 取得其上一時刻的所有隱藏狀態
3. 至時刻 t = T，我們可以開始倒推得到時刻 T 時的前一時點最佳隱藏狀態，以及該狀態出現的機率值
4. 一路從 t = T 倒推至 t = T-1, T-2, ..., 1，即得到最佳隱藏序列

第一步、首先，初始化所有狀態，這代表當沒有任何資訊的時候，我們只知道一個病人說自己正常的時候，其真實健康狀態的所有可能機率，比如根據過去資料判斷，這樣的病人六成健康，四成為生病，Backpointer 設定為 0 因為此時我們並不知道他來之前的真實狀態為何。

第二步、我們就開始計算啦！時刻 t = 2 的時候病人說自己發寒，發寒有兩種情況，一種是健康但是發寒，機率 0.4，一種是生病且發寒，機率 0.6，這樣我們的路徑可能有四種 : 健康且健康發寒、健康且生病發寒、生病且健康發寒、生病且生病發寒。根據計算我們可以得出 t = 2 最佳路徑為第一天健康且第二天健康但是發寒，此時我們一併紀錄 Backpointer，也就是如果 t = 2 的現在他觀察倒是健康的，那麼 t = 1 最有可能他也是健康的；如果觀察到生病，那麼 t = 1 最有可能他是健康的。根據這樣的思路往下遞迴計算到 t = T 的時候終止。

第三步、也就是 t = T 的當下，我們觀察到他說自己頭暈。這時候因為我們已經用動態規劃計算了 t = 1~T 的所有可能最佳路徑機率與狀態，我們可以往前倒推，由於 t = T 觀察到頭暈的時候，我們可以得到一組最有可能的路徑狀態 Ht，且當前最有可能的隱藏狀態根據 Backpointer，我們知道他是生病的。

第四步、有了Ht之後我們就可以倒推，在Ht中的每一個節點利用Backpointer可以得知在那當下最有可能的隱藏狀態，於是得到了最終結論。該病人可能是 {健康、健康、生病} 的真實健康狀態。

接下來，又進入到演算法實況，請參考 :

鮑姆威爾奇算法 Baum-Welch Algorithm

接下來進入到 HMM 的最後一個問題，即反向求解轉移參數，這個問題比前兩者更為複雜，因為我們想處理的角度通常是站在沒有隱藏序列的情況下去進行，這也更符合現實的情況，以上述那個病人的例子而言，即我有這個病人的闡述狀況後，我能否去得知 :

1. 健康與生病狀態下，產生各種症狀(頭痛、發寒、正常)的機率
2. 各種狀態中在下一期的轉移機率是多少? 比如一個正常的人，通常有多少機率下一期生病
3. 一般來到診所的病人中，開始處與健康與生病的機率各是多少

該算法用到了最大期望算法的原理，如果對該算法不熟的讀者建議參考 :

機器學習: EM 演算法(Expectation-Maximization Algorithm, EM)、高斯混合模型(Gaussian Mixture Model, GMM)和GMM-EM詳細推導

簡單來說，EM 演算法的精髓在於透過透過 E-step 與 M-step 兩者交替求解，不斷逼近函數可能的下界(E-Step)，並在給定下界的情況下對其進行最佳化(M-step)，藉此逼近目標函數。

放在 Baum-Welch 演算法，EM 算法的目標為求解一個聯合分布 P(O,H|λ)基於條件機率P(I|O,λ’)的期望值，其中λ’為當前模型參數，接著不斷透過 M-Step 更新模型參數直到收斂；直觀意涵是，當我們手頭資訊僅有觀測序列 O，以及一個被初始化的參數λ’時，我可以先透過 E-Step 找出該條件下的聯合分配，根據該分配去更新模型參數直到無法更新為止。

以上是簡單的表達式，接下來就要進入到實際計算的部分，首先，我們知道M-Step 可以改寫成 :

其中，此時的觀測數據與訓練數據由於轉變成 :

因此我們可以用如下方式改寫聯合分布 P(O, H|λ) :

將此式帶入上方 M-Step，得到 :

接下來就是大家都熟知的過程啦！只要透過MLE的思路，對A, B, I 分別求其導數就可以得到估計參數了，這邊我就不慢慢寫出來，直接寫求導結果(絕對不是因為我太久沒推微積分)，有興趣的讀者可以自己推導看看 :

其中 γ, δ 可以透過前向後向算法得出。至此，整個推導過程結束，而要實現 B-W 算法，可以按照以下的步驟寫程式碼 :

1. 初始化所有參數 A, B, I
2. 透過前向後向算法計算所有 γ, δ，其中 γ = [d, t, i]；δ = [d, t, i, j]
3. 更新模型參數 A, B, I 直到收斂為止 (收斂可以自行定義，或直接設定迭代最大次數)；如果不收斂，則回到第二步重新計算後更新

接下來我們一樣進入到程式碼環節，但這邊我要偷懶一下，直接貼上 Package 連結 :

hmmlearn/hmmlearn

這個 Package 是從 Scikit-Learn 移出來獨立的，可靠性相對好一些，為甚麼我不自己實作呢 ? 因為我寫著寫著，發現以 B-W 算法需要的資料結構完整性、算法與迭代複雜性等因素，如果不透過多層 Class 去包裝實在很難寫出讓人看得懂又可以簡單 Implement 的 Code (或是我 Coding 技術太差 QQ)，我在 Medium 上也沒有搜尋到有人將其完整寫出，故證明不是我懶惰不想花時間寫而是真的很麻煩(威~)

圖示化請參考 Wikipedia，我認為寫得很清楚。

Baum-Welch algorithm

結尾

又是幾百年沒寫，想單純扯些廢話。

文章內的這些東西實際上是台大應數所某堂課的回家作業，當時我依稀記得教授只講完了隱馬可夫鏈，讓課上同學盯著黑板看半天，某些學過演算法的同學在下面恍然大悟說可以DFS, DP把路徑走完，他笑著點了點頭。就讓我們回去把這三大算法給搞出來……

其實這也就算了，本來自己找材料研究然後做一遍就是研究中的重要精神，沒什麼可抱怨的，不過記得我當時連甚麼隱馬可夫都沒聽過，教授出作業的時候也是直接給了一組模擬資料和情境敘述，壓根沒說這算法的名稱叫甚麼，你還一時半會下不了關鍵字，一度看旁邊研究所的同學都老神在在的時候（而且同系都會互相分享作業！），我真的對著自己的智商懷疑了很久(笑)。不過好處也是，我當時就真的把這些算法的核心思維也看了一遍，其中每步運算牽涉到的基礎定理也逐一複習，雖然每天都熬夜到三點多寫程式，最後作業也是做了出來，頗有一番收穫就是了。(不過 B-W 算法我就完全不管 Code 完整性了，寫了超多個腳本分別運算哈哈，寫了大概兩百多行，不過這也是時間趕以及當時對資料結構的架設靈活度不夠高的因素)

另外就是，教授看到學生名單的時候不知道有沒有滿頭問號 : 「咦? 你這國企的跑來湊甚麼熱鬧?」他大概覺得我沒被勸退也是很神奇，不過管院每年似乎都有這麼一些人走上奇怪的路線就是了，當時同班我記得也有一個財金系的，只能說不管往哪條路走大概都不至於孤單，會有一些奇葩陪你XD；順道一提中研院來的教授好像都有一種，我隨便講反正有人會聽懂的上課哲學，100行的東西可以濃縮成3行講完。