Matlab統計學工具箱之（隱）馬爾可夫模型：Markov Models

原文：https://chunqiu.blog.ustc.edu.cn/?p=454

目錄

馬爾科夫鏈

隱馬爾科夫模型（HMM）

HMM簡介

分析隱馬爾科夫模型

產生測試序列

估計狀態序列

估計狀態轉移矩陣和輸出矩陣

估計后驗狀態概率

改變初始狀態分布

 

馬爾科夫鏈

         如下圖的馬爾科夫模型狀態轉移矩陣

 

         圖中的方框代表你要建模過程的可能狀態，箭頭表示狀態間的轉移情況。箭頭上的標簽表示轉移概率。在過程中的每一步，模型會產生一個輸出（決定於所處的狀態）然后轉移到另一個狀態。馬爾科夫模型的一個重要特性是下一個狀態只與當前的狀態有關，與過去的轉移無關。

         例如，對於一個順序拋擲硬幣的過程，存在兩個狀態：正和反。最近一次的投擲決定了模型當前的狀態，下一次投擲決定了如何轉移到下一狀態。如果硬幣是公平的，則轉移概率都是二分之一，輸出就是當前的狀態。在更復雜的模型中，每一個狀態可以以一個隨機過程產生輸出，如可以通過投擲骰子來決定輸出。

         馬爾科夫鏈是離散狀態馬爾科夫模型的數學描述。馬爾科夫鏈如下定義：

l  一系列狀態{1,2,…,M}；

l  一個MXM的狀態轉移矩陣T，其中Tij表示從狀態i到j的狀態轉移概率。矩陣每行的和為1，因為它表示從一個給定狀態到其他各狀態的所有轉移概率的和。

l  一系列可能的輸出{s1,s2,…,sN}，默認時可以當做{1,2,…,N}，N為可能輸出的個數，當然你可以選擇不同的數據集或符號來表示；

l  一個MXN的輸出矩陣E，其中Eik表示模型在狀態i時的輸出sk的概率。

馬爾科夫鏈的過程：從第0步的初始狀態i0開始；然后以概率T1i1轉移到狀態i1，並以概率Ei1k1輸出sk1；最后，以概率T1i1 Ei1k1 Ti1i2 Ei2k2 … Tir-1ir Eirk觀察到狀態i1i2…ir，同時輸出sk1sk2…skr。

 

隱馬爾科夫模型（HMM）

HMM簡介

         隱馬爾科夫模型是你只觀察到輸出，但不知道模型產生輸出所經歷的狀態。隱馬爾科夫模型分析就是從觀察的數據恢復相關的狀態序列。

         例如，考慮一個有2個狀態6個可能輸出的馬爾科夫模型，模型中用到

n  一個紅骰子，6個面，分別標號1~6；

n  一個綠骰子，12個面，5個面標號2~6，其他面標記1；

n  一個非均衡的紅硬幣，以0.9的概率擲出正，0.1的概率擲出反；

n  一個非均衡的綠硬幣，以0.95的概率擲出正，0.05的概率擲出反。

模型以下面的規則產生一系列來自{1,2,3,4,5,6}的數字序列：

n  先擲紅骰子，記下骰子的輸出；

n  然后擲紅硬幣：

l  如果是正，則擲紅骰子記下結果；

l  否則擲綠骰子記結果；

n  接下來的下一步，擲與上一步所用骰子同顏色的硬幣，如果是正，則擲同顏色的骰子；否則擲另一個骰子。

模型的狀態轉移圖如下

 

         這個模型中，通過擲與狀態同顏色的骰子來決定輸出，通過擲與狀態同顏色的硬幣來決定轉移到哪一狀態。

         狀態轉移矩陣為

 

         輸出矩陣為

 

         這個模型是非隱的，因為你從硬幣、骰子的顏色可以知道狀態序列。如果假設有人產生了一系列輸出，但不給你展示硬幣和骰子，所有你能看到的只是輸出序列。如果你開始看到1出現的次數比其他都多，你可能會假設模型處於綠狀態，但你不能確定，因為你不能看到所擲骰子的顏色。

         隱馬爾科夫模型提出了如下問題：

1)         給定輸出序列，求最可能的狀態序列？

2)         給定輸出序列，你如何估計模型的轉移概率和輸出概率？

3)         模型得出給定序列的先驗概率？

4)         模型在序列某一位置是特定狀態的后驗概率？

分析隱馬爾科夫模型

         與隱馬爾科夫模型有關的統計工具箱函數如下

l  Hmmgenerate——從一個馬爾科夫模型產生狀態、輸出序列；

l  Hmmestimate——從輸出序列和已知的狀態序列計算轉移、輸出序列概率的極大似然估計；

l  Hmmtrain——從輸出序列計算轉移、輸出序列概率的極大似然估計；

l  Hmmviterbi——計算隱馬爾科夫模型的最大可能狀態序列；

l  Hmmdecode——計算輸出序列的狀態后驗概率。

下面來講如何使用這些函數分析隱馬爾科夫模型。

產生測試序列

         下面的命令創建狀態轉移矩陣和輸出矩陣

TRANS = [.9 .1; .05 .95;];

 

EMIS = [1/6, 1/6, 1/6, 1/6, 1/6, 1/6;...

7/12, 1/12, 1/12, 1/12, 1/12, 1/12];

        要產生一組狀態和輸出的隨機序列，使用函數hmmgenerate

[seq,states] = hmmgenerate(1000,TRANS,EMIS);

輸出參數seq為輸出序列，states為狀態序列。

         Hmmgenerate函數第0步開始於狀態1，然后第一步轉移到狀態i1，而i1作為states的第一個元素。如果要改變初始狀態，看后面的改變初始狀態分布一節。

估計狀態序列

         給定狀態轉移矩陣TRANS和輸出矩陣EMIS，函數hmmviterbi使用Viterbi算法計算產生給定輸出序列seq的極大可能狀態序列

likelystates = hmmviterbi(seq, TRANS, EMIS);

likelystates與seq長度相同。

         為了測試hmmviterbi的精確度，計算實際狀態序列states與估計狀態序列likelystates吻合的百分比

sum(states==likelystates)/1000

ans =

   0.8200

         在這個例子中，估計的狀態和實際吻合達82%

估計狀態轉移矩陣和輸出矩陣

函數hmmestimate和hmtrain可以用來根據給定的輸出序列估計狀態轉移矩陣和輸出矩陣。

使用hmmestimate

         此函數需要你知道產生輸出序列的狀態序列。基於輸出序列和狀態序列，如下估計轉移矩陣和輸出矩陣

[TRANS_EST, EMIS_EST] = hmmestimate(seq, states)

 

TRANS_EST =

0.8989    0.1011

0.0585    0.9415

 

EMIS_EST =

0.1721    0.1721    0.1749    0.1612    0.1803    0.1393

0.5836    0.0741    0.0804    0.0789    0.0726    0.1104

         與原始的序列對比

TRANS

TRANS =

0.9000    0.1000

0.0500    0.9500

 

EMIS

EMIS =

0.1667    0.1667    0.1667    0.1667    0.1667    0.1667

0.5833    0.0833    0.0833    0.0833    0.0833    0.0833

使用hmmtrain

         如果你不知道狀態序列，但知道TRANS和EMIS矩陣的初始猜測，可以使用hmmtrain函數估計。

         假設已知TRANS和EMIS矩陣的初始猜測如下

TRANS_GUESS = [.85 .15; .1 .9];

EMIS_GUESS = [.17 .16 .17 .16 .17 .17;.6 .08 .08 .08 .08 08];

下面做矩陣估計

[TRANS_EST2, EMIS_EST2] = hmmtrain(seq, TRANS_GUESS, EMIS_GUESS)

 

TRANS_EST2 =

0.2286    0.7714

0.0032    0.9968

 

EMIS_EST2 =

0.1436    0.2348    0.1837    0.1963    0.2350    0.0066

0.4355    0.1089    0.1144    0.1082    0.1109    0.1220

         Hmmtrain函數使用迭代算法調整TRANS_GUESS和EMIS_GUESS，使得基於調整的矩陣產生更加類似於觀察序列seq的輸出。當兩次迭代中矩陣的變化不大時迭代結束。

         如果如果算法在默認的最大迭代次數100次以內仍沒結束，將返回迭代的最后結果並產生一個警告。這時可以通過增加最大迭代次數來增大迭代的次數

hmmtrain(seq,TRANS_GUESS,EMIS_GUESS,'maxiterations',maxiter)

maxiter是新設的最大迭代次數。

還可以改變默認的最小容忍值

hmmtrain(seq, TRANS_GUESS, EMIS_GUESS, 'tolerance', tol)

tol為新設的最小容忍值，增大它可以讓算法更快的停止，但結果將不精確。

         兩個因素將降低hmmtrain函數估計結果的可靠性：

l  算法收斂到一個局部極大值，但此極大值並不是真正的狀態轉移矩陣與輸出矩陣的位置。對此，可以使用多個不同的初始猜測值來估計，使用最好的結果；

l  Seq序列太短，以至於不足以訓練矩陣，對此，可以使用更長的seq序列。

估計后驗狀態概率

輸出序列seq的后驗狀態概率是模型給出seq序列時所處狀態的后驗概率。可以使用hmmdecode函數計算后驗狀態概率：

PSTATES = hmmdecode(seq,TRANS,EMIS)

         隨着序列長度的增加，序列的概率會趨向於0，所以一個足夠長的序列的概率將小於計算機的最小精度，為了避免此問題，hmmdecode還返回了概率的對數值。

改變初始狀態分布

         默認地，統計工具箱的隱馬爾科夫模型函數是從狀態1開始的。換句話說，初始狀態的分布集中在狀態1處，其他位置為0。如果想給初始狀態賦予一個概率分布p=[p1,p2,…,pM]，可以如下實現

1. 創建一個M+1XM+1的輔助狀態轉移矩陣如下

 

其中的T是真正的狀態轉移矩陣，第一列包含M+1個0，向量p的和為1。

1. 創建M+1XN的輔助輸出矩陣如下

 

         如果狀態轉移矩陣與輸出矩陣式TRANS和EMIS，可以使用如下命令創建輔助矩陣

TRANS_HAT = [0 p; zeros(size(TRANS,1),1) TRANS];

 EMIS_HAT = [zeros(1,size(EMIS,2)); EMIS];

 

原文：http://www.aiseminar.cn/bbs/forum.php?mod=viewthread&tid=596

此文講述的內容在Matlab 7.0、7.5（R2007b）中均有――馬爾可夫工具箱，主要內容如下。

簡介：馬爾可夫處理是隨機處理的一個典型例子――此種處理根據特定的概率產生隨機輸出或狀態序列。馬爾可夫處理的特別之處在於它的無記憶性――他的下一個狀態僅依賴他的當前狀態，不考慮導致他們的歷史。馬爾可夫處理的模型在實際應用中使用非常廣泛，從每日股票價格到染色體中的基因位置都有應用。

馬爾可夫鏈

馬爾可夫模型用狀態圖可視化描述如下。
<ignore_js_op> 

在圖中，矩形代表你要描述的模型在處理中可能出現的狀態，箭頭描述了狀態之間的轉換。每個箭頭上的標簽表明了該轉換出現的概率。在處理的每一步，模型都可能根據當前狀態產生一種output或emission，然后做一個到另一狀態的轉換。馬爾可夫模型的一個重要特點是：他的下個狀態僅僅依賴當前狀態，而與導致其成為當前狀態的歷史變換無關。

馬爾可夫鏈是馬爾可夫模型的一組離散狀態集合的數學描述形式。馬爾可夫鏈特征歸納如下：

1. 一個狀態的集合{1, 2, ..., M}
2. 一個M * M的轉移矩陣T，(i, j)位置的數據是從狀態i轉到狀態j的概率。T的每一行值的和必然是1，因為這是從一個給定狀態轉移到其他所有可能狀態的概率之和。
3. 一個可能的輸出（output）或發布（emissions）的集合{S1, S2, ..., SN}。默認情況下，發布的集合是{1, 2, ..., N}，這里N是可能的發布的個數，當然，你也可以選擇一個不同的數字或符號的集合。
4. 一個M * N的發布矩陣E，(i, k)入口給出了從狀態i得到發布的標志Sk的概率。

馬爾可夫鏈在第0步，從一個初始狀態i0開始。接着，此鏈按照T(1, i1)概率轉移到狀態i1，且按概率E(i1, k1)概率發布一個輸出S(k1)。因此，在前r步，狀態序列(i1, i2, i3, ..., ir)和發布序列(Sk1, Sk2, ..., Skr)的可能的觀測結果是
T(1, i1)E(i, k1), T(i1, i2)E(i2, k2), ..., T(ir-1, ir)E(ir, k)


隱馬爾可夫模型

簡介：隱馬爾可夫模型相對馬爾可夫模型的不同之處在於，你觀測到一組發布序列，但是卻不知道模型通過什么樣的狀態序列得到這樣的發布。隱馬爾可夫模型分析就是要從觀測數據恢復出這一狀態序列。

一個例子：考慮一個擁有2個狀態和6個位置發布的馬爾可夫模型。模型描述如下：

1. 一個紅色骰子，有6個面，標記為1到6。
2. 一個綠色骰子，有12個面，有5個標記為2到6，余下的全標記為1。
3. 一個加權的紅硬幣，正面的概率是0.9，背面的概率是0.1。
4. 一個加權的綠硬幣，正面的概率為0.95，背面的概率為0.05。

這個模型按照下面的規則從集合{1, 2, 3, 4, 5, 6}產生一個數字序列：

1. 從滾動紅骰子開始，記下出現的數字，作為發布結果。
2. 投紅硬幣：
如果結果是正面，滾動紅骰子，記下結果；
如果結果是背面，滾動綠骰子，記下結果。
3. 在下面的每一步，你都拋和前一步所投骰子相同顏色的硬幣。如果硬幣出現正面，滾和硬幣相同顏色的骰子。如果硬幣出現反面，改為投另種顏色的骰子。

這個模型的狀態圖如下，有兩個狀態，紅和綠：
<ignore_js_op> 
通過投擲和狀態一樣顏色的骰子來解決輸出（發布），通過拋同樣顏色的硬幣來決定狀態的轉移。

狀態轉移概率矩陣為：
>> T = [0.9 0.1; 0.05 0.95]
T =
    0.9000    0.1000
    0.0500    0.9500

輸出或發布概率矩陣為：
>> E = [1/6 1/6 1/6 1/6 1/6 1/6; 7/12 1/12 1/12 1/12 1/12 1/12]
E =
    0.1667    0.1667    0.1667    0.1667    0.1667    0.1667
    0.5833    0.0833    0.0833    0.0833    0.0833    0.0833

這個模型並不是隱藏的，因為我們從硬幣和骰子的顏色已經知道狀態序列。假設，有其他人產生了一個發布結果，而沒有向你展示硬幣和骰子，你能看到的只有結果。當你看到1比其他數字多時，你也許猜測這個模型是在綠色狀態，但是因為你不能看到被投骰子的顏色，所以你並不能確定。

隱馬爾可夫模型提出了如下問題：

1. 給定一個輸出序列，最有可能的狀態過程是什么？
2. 給定一個輸出序列，如何估計模型的轉移和發布概率？
3. 如何求這個模型產生一個給定序列的先驗概率？
4. 如何求這個模型產生一個給定序列的后驗概率？


隱馬爾可夫模型分析

統計工具箱包含5個與隱馬爾可夫模型相關的函數：
hmmgenerate 從一個馬爾可夫模型產生一個狀態序列和輸出序列；
hmmestimate 計算轉移和輸出的極大似然估計；
hmmtrain 從一個輸出序列計算轉移和輸出概率的極大似然估計；
hmmviterbi 計算一個隱馬爾可夫模型最可能的狀態變化過程；
hmmdecode 計算一個給定輸出序列的后驗狀態概率。

下面部分介紹如何使用這些函數來分析隱馬爾可夫模型。

1. 產生一個測試序列

下面代碼產生上面簡介中模型的轉移和輸出矩陣：
TRANS = [.9 .1; .05 .95;];
EMIS = [1/6, 1/6, 1/6, 1/6, 1/6, 1/6;...
7/12, 1/12, 1/12, 1/12, 1/12, 1/12];

要從模型產生一個隨機的狀態序列和輸出序列，使用hmmgenerate：
[seq,states] = hmmgenerate(1000,TRANS,EMIS);

輸出中，seq是輸出序列，states是狀態序列。hmmgenerate在第0步從狀態1開始，在第一步轉移到狀態i1
，並返回i1作為狀態的第一個入口。

2. 估計狀態序列

給定了轉移和輸出矩陣TRANS和EMIS，函數hmmviterbi使用Viterbi算法計算模型給定輸出序列seq最有可能
通過的狀態序列：
likelystates = hmmviterbi(seq, TRANS, EMIS);

likelystates是和seq一樣長的序列。計算hmmvertibi的精度如下：
sum(states == likelystates) / length(states)
ans =
    0.8680

3. 估計轉移和輸出矩陣
函數hmmestimate和hmmtrain用於估計給定輸出序列seq的轉移和輸出矩陣TRANS和EMIS。

使用hmmestimate
[TRANS_EST, EMIS_EST] = hmmestimate(seq, states)
TRANS_EST =
    0.9065    0.0935
    0.0406    0.9594
EMIS_EST =
    0.1452    0.1516    0.1581    0.1968    0.1581    0.1903
    0.5841    0.0754    0.0986    0.0812    0.0841    0.0768
由上面使用方式可知，hmmestimate函數需要事先知道了得到輸出序列seq，以及得到此結果的狀態變化序
列。

使用hmmtrain
如果不知道狀態序列，但是知道TRANS和EMIS的初始猜測，那就可以使用hmmtrain來估計TRANS和EMIS。

假設已知如下初始猜測：
TRANS_GUESS = [.85 .15; .1 .9];
EMIS_GUESS = [.17 .16 .17 .16 .17 .17;.6 .08 .08 .08 .08 08];

TRANS和EMIS的估計如下：
[TRANS_EST2, EMIS_EST2] = hmmtrain(seq, TRANS_GUESS, EMIS_GUESS)
TRANS_EST2 =
    0.9207    0.0793
    0.0370    0.9630
EMIS_EST2 =
    0.1792    0.1437    0.1436    0.1855    0.1509    0.1971
    0.5774    0.0775    0.1042    0.0840    0.0859    0.0710

hmmtrain使用迭代算法來不斷修改TRANS_GUESS和EMIS_GUESS，使得每一步修改得到的矩陣都更加可能產生觀測序列seq。當前后兩個兩次迭代矩陣的變化在一個小的容錯范圍內時，迭代停止。如果算法無法達到容錯的范圍，則迭代到達一定次數時就會停止，並返回一個警告提示。默認的最大迭代次數為100。

如果算法達不到目標誤差范圍，則可以通過增加迭代次數和/或加大容錯誤差值來使其獲得較合適結果：
改變迭代次數maxiter：hmmtrain(seq,TRANS_GUESS,EMIS_GUESS,'maxiterations',maxiter)
改變容錯誤差tol：hmmtrain(seq, TRANS_GUESS, EMIS_GUESS, 'tolerance', tol)

影響hmmtrain輸出的矩陣可靠性的兩點因素：
（1）算法收斂於局部極值，這點可以使用不同的初始猜測矩陣來嘗試解決；
（2）序列seq太短而無法很好的訓練矩陣，可以嘗試使用較長的序列。

4. 估計后驗狀態概率（不太理解）
一個輸出序列seq的后驗狀態概率是在特定狀態下的模型產生在seq中一個輸出的條件概率。假定seq已經給出，你可以使用hmmdecode得到后驗狀態概率。

PSTATES = hmmdecode(seq,TRANS,EMIS)

輸出為一個M * N的矩陣。M是狀態的個數，L是seq的長度。PSTATES(i, j)是模型在狀態i時，產生seq第j個輸出的條件概率。

參考文獻：
[1] Matlab R2007b幫助文檔，http://www.aiseminar.cn/bbs，jink2005簡譯！

 

 

----------------------------------------------------------------------------------------------------------------

統計工具箱包含5個與隱馬爾可夫模型相關的函數：
hmmgenerate 從一個馬爾可夫模型產生一個狀態序列和輸出序列；
hmmestimate 計算轉移和輸出的極大似然估計；
hmmtrain 從一個輸出序列計算轉移和輸出概率的極大似然估計；
hmmviterbi 計算一個隱馬爾可夫模型最可能的狀態變化過程；
hmmdecode 計算一個給定輸出序列的后驗狀態概率。

下面部分介紹如何使用這些函數來分析隱馬爾可夫模型。

1. 產生一個測試序列

下面代碼產生上面簡介中模型的轉移和輸出矩陣：
TRANS = [.9 .1; .05 .95;];
EMIS = [1/6, 1/6, 1/6, 1/6, 1/6, 1/6;...
7/12, 1/12, 1/12, 1/12, 1/12, 1/12];

要從模型產生一個隨機的狀態序列和輸出序列，使用hmmgenerate：
[seq,states] = hmmgenerate(1000,TRANS,EMIS);

輸出中，seq是輸出序列，states是狀態序列。hmmgenerate在第0步從狀態1開始，在第一步轉移到狀態i1
，並返回i1作為狀態的第一個入口。

2. 估計狀態序列

給定了轉移和輸出矩陣TRANS和EMIS，函數hmmviterbi使用Viterbi算法計算模型給定輸出序列seq最有可能
通過的狀態序列：
likelystates = hmmviterbi(seq, TRANS, EMIS);

likelystates是和seq一樣長的序列。計算hmmvertibi的精度如下：
sum(states == likelystates) / length(states)
ans =
    0.8680

3. 估計轉移和輸出矩陣
函數hmmestimate和hmmtrain用於估計給定輸出序列seq的轉移和輸出矩陣TRANS和EMIS。

使用hmmestimate
[TRANS_EST, EMIS_EST] = hmmestimate(seq, states)
TRANS_EST =
    0.9065    0.0935
    0.0406    0.9594
EMIS_EST =
    0.1452    0.1516    0.1581    0.1968    0.1581    0.1903
    0.5841    0.0754    0.0986    0.0812    0.0841    0.0768
由上面使用方式可知，hmmestimate函數需要事先知道了得到輸出序列seq，以及得到此結果的狀態變化序
列。

使用hmmtrain
如果不知道狀態序列，但是知道TRANS和EMIS的初始猜測，那就可以使用hmmtrain來估計TRANS和EMIS。

假設已知如下初始猜測：
TRANS_GUESS = [.85 .15; .1 .9];
EMIS_GUESS = [.17 .16 .17 .16 .17 .17;.6 .08 .08 .08 .08 08];

TRANS和EMIS的估計如下：
[TRANS_EST2, EMIS_EST2] = hmmtrain(seq, TRANS_GUESS, EMIS_GUESS)
TRANS_EST2 =
    0.9207    0.0793
    0.0370    0.9630
EMIS_EST2 =
    0.1792    0.1437    0.1436    0.1855    0.1509    0.1971
    0.5774    0.0775    0.1042    0.0840    0.0859    0.0710

hmmtrain使用迭代算法來不斷修改TRANS_GUESS和EMIS_GUESS，使得每一步修改得到的矩陣都更加可能產生觀測序列seq。當前后兩個兩次迭代矩陣的變化在一個小的容錯范圍內時，迭代停止。如果算法無法達到容錯的范圍，則迭代到達一定次數時就會停止，並返回一個警告提示。默認的最大迭代次數為100。

如果算法達不到目標誤差范圍，則可以通過增加迭代次數和/或加大容錯誤差值來使其獲得較合適結果：
改變迭代次數maxiter：hmmtrain(seq,TRANS_GUESS,EMIS_GUESS,'maxiterations',maxiter)
改變容錯誤差tol：hmmtrain(seq, TRANS_GUESS, EMIS_GUESS, 'tolerance', tol)

影響hmmtrain輸出的矩陣可靠性的兩點因素：
（1）算法收斂於局部極值，這點可以使用不同的初始猜測矩陣來嘗試解決；
（2）序列seq太短而無法很好的訓練矩陣，可以嘗試使用較長的序列。

4. 估計后驗狀態概率（不太理解）
一個輸出序列seq的后驗狀態概率是在特定狀態下的模型產生在seq中一個輸出的條件概率。假定seq已經給出，你可以使用hmmdecode得到后驗狀態概率。

PSTATES = hmmdecode(seq,TRANS,EMIS)

輸出為一個M * N的矩陣。M是狀態的個數，L是seq的長度。PSTATES(i, j)是模型在狀態i時，產生seq第j個輸出的條件概率。