隱馬爾可夫模型(HMM)

原文地址：http://www.cnblogs.com/jacklu/p/7753471.html

本文結合了王曉剛老師的ENGG 5202 Pattern Recognition課程內容知識，和搜集的資料和自己理解的總結。

1 概述

隱馬爾可夫模型（Hidden Markov Model，HMM）是結構最簡單的貝葉斯網，這是一種著名的有向圖模型，主要用於時序數據建模（語音識別、自然語言處理等數據在時域有依賴性的問題）。

如果考慮t時刻數據依賴於0到t-1時間段的所有數據，即，在計算復雜度上是不可行的。因此Markov Model假定只依賴於最近的幾個觀測數據。

下面先從一個直觀的例子理解HMM：

假設有三個不同的骰子（6面、4面、8面），每次先從三個骰子里選一個，每個骰子選中的概率為，如下圖所示，重復上述過程，得到一串數字[1 6 3 5 2 7]。這些可觀測變量組成可觀測狀態鏈。

同時，在隱馬爾可夫模型中還有一條由隱變量組成的隱含狀態鏈，在本例中即骰子的序列。比如得到這串數字骰子的序列可能為[D6 D8 D8 D6 D4 D8]。

隱馬爾可夫模型示意圖如下所示：

圖中，箭頭表示變量之間的依賴關系。在任意時刻，觀測變量（骰子點數）僅依賴於狀態變量（哪類骰子），“觀測獨立性假設”。

同時，t時刻數據依賴於t-1時刻的數據。這就是1階馬爾可夫鏈，即系統的下一時刻的狀態僅由當前狀態決定不依賴以往的任何狀態（無記憶性），“齊次馬爾可夫性假設”。

 

0階Markov Model：

1階Markov Model：

 

2階Markov Model：

 

1階HMM

包含狀態變量（也叫latent variable，該變量是離散的、未知的、待推斷的）和觀測變量（該變量可以是離散的、也可以是連續的）,如下圖所示：

其聯合分布：

1.2 HMM中的條件獨立（在后續算法推導中非常重要）

從概率圖模型上給出條件獨立的式子非常簡單，即遮住某一節點，被分開的路徑在給定該節點時獨立。

上面六個式子，前五個式子很容易從圖模型中理解。最后一個式子可以將左邊寫成和的乘積，然后再將做分解。

 

假定每個狀態有三種取值，比如上面骰子的種類。參數如下圖所示：

初始狀態參數

狀態轉移概率,即

觀測概率（也叫emission probablity），即時刻t、狀態的概率

2 隱馬爾可夫模型三要素

以上三個參數構成隱馬爾可夫模型三要素：

狀態轉移概率矩陣A，

觀測概率矩陣B，

初始狀態概率向量

一個隱馬爾可夫模型可由來指代。

3 隱馬爾可夫模型的三個基本問題

（1） 給定模型，計算其產生觀測序列的概率, 稱作evaluation problem，比如：計算擲出點數163527的概率

（2） 給定模型和觀測序列，推斷能夠最大概率產生此觀測序列的狀態序列，即使求解，稱作decoding problem，比如：推斷擲出點數163527的骰子種類

（3） 給定觀測序列，估計模型的參數，使計算其產生觀測序列的概率最大，稱作learning problem，比如：已知骰子有幾種，不知道骰子的種類，根據多次擲出骰子的結果，反推出骰子的種類

這三個基本問題在現實應用中非常重要，例如根據觀測序列推測當前時刻最有可能出現的觀測值，這就是基本問題（1）；

在語音識別中，觀測值為語音信號，隱藏狀態為文字，根據觀測信號推斷最有可能的狀態序列，即基本問題(2)；

在大多數應用中，人工指定參數模型已變得越來越不可行，如何根據訓練樣本學得最優參數模型，就是基本問題(3)。

4 三個基本問題的解法

基於兩個條件獨立假設，隱馬爾可夫模型的這三個基本問題均能被高效求解。

4.1 基本問題（1）evaluation problem解法

4.1.1 直接計算法(概念上可行，計算上不可行)

通過列舉所有可能的長度為T的狀態序列，求各個狀態序列與觀測序列同時出現的聯合概率，然后對所有可能求和。

計算復雜度，C是狀態個數。算法不可行。

4.1.2 前向算法(t=1,一步一步向前計算)

前向概率，表示模型，時刻 t，觀測序列為且狀態為的概率。

注意求和式中有K項（Z的狀態數），計算復雜度為C*C。

通過上式可知，為了得到前向概率，可以先初始化t=1時刻的概率，然后從第一個節點開始遞推計算，每次遞推都需要計算一次c*c的的操作，因此總的算法復雜度是（C和K相同）

4.1.3 后向算法

后向概率，表示模型，時刻 t，觀測序列為且狀態為的概率。

推導過程：

通過上式可知，為了得到后向概率，可以先初始化t=T時刻的概率，然后從最后一個節點向前遞推計算，每次遞推都需要計算一次c*c的的操作，因此總的算法復雜度是（C和K相同）

算法高效的關鍵是其局部計算前向概率，根據路徑結構，如下圖所示，每次計算直接利用前一時刻計算結果，避免重復計算，減少計算量。

利用前向概率和后向概率可以計算：

整個觀測序列的概率：

給定觀測序列，t時刻的狀態后驗概率：

給定觀測序列，t時刻從某一狀態，在t+1時刻轉換成新的狀態的后驗概率：

4.2 基本問題（2）decoding problem解法

4.2.1 近似算法

選擇每一時刻最有可能出現的狀態，即根據上述計算t時刻的狀態后驗概率，選擇概率最大的狀態，從而得到一個狀態序列。這個方法計算簡單，此方法但是不能保證整個狀態序列的出現概率最大。因為可能兩個相鄰的狀態轉移概率為0，即實際上不可能發生這種狀態轉換。

4.2.2 Viterbi算法

使用動態規划求解概率最大（最優）路徑。t=1時刻開始，遞推地計算在時刻t狀態為i的各條部分路徑的最大概率，直到計算到時刻T，狀態為i的各條路徑的最大概率，時刻T的最大概率即為最優路徑的概率，最優路徑的節點也同時得到。

如果還不明白，看一下李航《統計學習方法》的186-187頁的例題就能明白算法的原理。

考慮一個表格數據結構，存儲着t時刻時，狀態為j的能夠產生觀測序列的最大概率值。

t=1時，

t>1時

          

維特比算法：

使用記錄求解的狀態序列。

維特比算法圖示：

狀態[3 3 3]極為概率最大路徑。

4.3 基本問題（3）解法

4.3.1 監督學習方法

給定T個長度相同的（觀測序列，狀態序列）作為訓練集，使用極大似然估計法來估計模型參數。

轉移概率 的估計:樣本中t時刻處於狀態i，t+1時刻轉移到狀態j的頻數為，則

觀測概率和初始狀態概率的估計類似。

4.3.2 Baum-Welch算法

使用EM算法得到模型參數估計式

EM算法是常用的估計參數隱變量的利器，它是一種迭代方法，基本思想是：

（1） 選擇模型參數初始值；

（2） (E步)根據給定的觀測數據和模型參數，求隱變量的期望；

（3） (M步)根據已得隱變量期望和觀測數據，對模型參數做極大似然估計，得到新的模型參數，重復第二步。

作業題：

第一問用於理解HMM產生數據的過程，第二問用於理解維特比算法。

自己寫的答案（運行結果如上圖）：

import numpy
import random

def random_pick(pick_list,probability_list):
    x=random.uniform(0,1)
    cumulative_probability=0.0
    for item,item_probability in zip(pick_list,probability_list):
        cumulative_probability+=item_probability
        if x < cumulative_probability: break
    return item

Zt=[1,2]
Pi=[0.6,0.4]
Xt=[1,2,3]
a1j=[0.7, 0.3]
a2j=[0.4, 0.6]
b1j=[0.1, 0.4, 0.5]
b2j=[0.6, 0.3, 0.1]
x=[-1 for n in range(10)]
z=[-1 for n in range(10)]
#for function test
#temp_counter = 0
#for i in range(100):
#    if random_pick(Zt,Pi) == 1: temp_counter+=1
#print(temp_counter)
##for function test

z[0] = random_pick(Zt,Pi)
for i in range(10):
    if z[i] == 1:
        x[i] = random_pick(Xt, b1j)
        if i < 9: z[i+1] = random_pick(Zt, a1j)
    else:
        x[i]= random_pick(Xt, b2j)
        if i < 9: z[i+1]= random_pick(Zt, a2j)
print(z)
print(x)

bp=[-1 for n in range(10)]
Fi=[[-1,-1] for n in range(10)]
for i in range(2):
    if i == 0:
        Fi[0][0]=Pi[i]*b1j[x[0]-1]
        print('Fi00',Fi[0][0])
    elif i == 1:
        Fi[0][1]=Pi[i]*b2j[x[0]-1]
        print('Fi01',Fi[0][1])
if Fi[0][0] < Fi[0][1]:
    bp[0] = 1
else:
    bp[0] = 0
    
for t in range(9):
    for j in range(2):
        if j == 0:
            if bp[t] == 0: 
                 Fi[t+1][j]=Fi[t][0]*a1j[0]*b1j[x[t+1]-1]
            elif bp[t] == 1:
                Fi[t+1][j]=Fi[t][1]*a2j[0]*b1j[x[t+1]-1]
        if j == 1:
            if bp[t] == 0: 
                Fi[t+1][j]=Fi[t][0]*a1j[1]*b2j[x[t+1]-1]
            elif bp[t] == 1:
                Fi[t+1][j]=Fi[t][1]*a2j[1]*b2j[x[t+1]-1]
    print('Fit0',Fi[t+1][0])
    print('Fit1',Fi[t+1][1])
    if Fi[t+1][0] < Fi[t+1][1]:
        bp[t+1] = 1
    else:
        bp[t+1] = 0

print(bp)#z=bp+1

 

 

參考資料：

CUHK 王曉剛老師的ENGG 5202 Pattern Recognition 課堂講義

《機器學習》周志華

《統計學習方法》李航

如何用簡單易懂的例子解釋隱馬爾可夫模型https://www.zhihu.com/question/20962240

《Pattern Classification》

《PRML》