隱馬爾可夫模型差不多是學(xué)習(xí)中遇到的最難的模型了，本節(jié)通過對《統(tǒng)計(jì)學(xué)習(xí)方法》進(jìn)行學(xué)習(xí)并結(jié)合網(wǎng)上筆記，用Python代碼實(shí)現(xiàn)了隱馬模型觀測序列的生成、前向后向算法、Baum-Welch無監(jiān)督訓(xùn)練、維特比算法。比較清晰的了解了隱馬爾可夫模型，其實(shí)在實(shí)際運(yùn)用中我們只需要調(diào)用庫就一行代碼解決問題。

在調(diào)用分詞函數(shù)時(shí)就采用的隱馬爾可夫模型HMM原理來進(jìn)行切詞。

具體原理解釋介紹如下：

1 隱馬爾可夫模型定義

隱馬爾可夫模型是關(guān)于時(shí)序的概率模型，描述由一個(gè)隱藏的馬爾可夫鏈隨機(jī)生成不可觀測的狀態(tài)隨機(jī)序列，再由各個(gè)狀態(tài)生成一個(gè)觀測而產(chǎn)生觀測隨機(jī)序列的過程。隱藏的馬爾可夫鏈隨機(jī)生成的狀態(tài)的序列，稱為狀態(tài)序列（state sequence);每個(gè)狀態(tài)生成一個(gè)觀測，而由此產(chǎn)生的觀測的隨機(jī)序列，稱為觀測序列（observation sequence)。序列的每一個(gè)位置又可以看作是一個(gè)時(shí)刻。

什么叫隱馬爾科夫鏈呢？簡單說來，就是把時(shí)間線看做一條鏈，每個(gè)節(jié)點(diǎn)只取決于前N個(gè)節(jié)點(diǎn)。就好比你打開朋友圈通過你朋友的最近幾條狀態(tài)就可以判斷他接下來要干嘛。

接下來引入一些符號(hào)來表述定義：

設(shè)Q是所有可能的狀態(tài)的集合，V是所有可能的觀測的集合。

其中，N是可能的狀態(tài)數(shù)，M是可能的觀測數(shù)。

狀態(tài)q不可見，觀測v可見。應(yīng)用到詞性標(biāo)注系統(tǒng)，詞就是v，詞性就是q。

I是長度為T的狀態(tài)序列，O是對應(yīng)的觀測序列。

這可以理解為相當(dāng)于給定了一個(gè)詞（O）+詞性（I）的訓(xùn)練集或者說是中文分詞系統(tǒng)里的詞典。

定義A為狀態(tài)轉(zhuǎn)移概率矩陣：

其中，

是在時(shí)刻t處于狀態(tài)qi的條件下在時(shí)刻t+1轉(zhuǎn)移到狀態(tài)qj的概率。

這實(shí)際在表述一個(gè)一階的HMM，所作的假設(shè)是每個(gè)狀態(tài)只跟前一個(gè)狀態(tài)有關(guān)。

B是觀測概率矩陣:

其中，

是在時(shí)刻t處于狀態(tài)qj的條件下生成觀測vk的概率。

π是初始狀態(tài)概率向量：

其中，

是時(shí)刻t=1處于狀態(tài)qj的概率。

隱馬爾可夫模型由初始狀態(tài)概率向量π、狀態(tài)轉(zhuǎn)移概率矩陣A和觀測概率矩陣B決定。π和A決定狀態(tài)序列，B決定觀測序列。因此，隱馬爾可夫模型可以用三元符號(hào)表示，即

括號(hào)里的三個(gè)變量稱為隱馬爾可夫模型的三要素。加上一個(gè)具體的狀態(tài)集合Q和觀測序列V，則構(gòu)成了HMM的五元組。

狀態(tài)轉(zhuǎn)移概率矩陣A與初始狀態(tài)概率向量π確定了隱藏的馬爾可夫鏈，生成不可觀測的狀態(tài)序列。觀測概率矩陣B確定了如何從狀態(tài)生成觀測，與狀態(tài)序列綜合確定了如何產(chǎn)生觀測序列。

從定義可知，隱馬爾可夫模型作了兩個(gè)基本假設(shè)：

(1)齊次馬爾可夫性假設(shè)，即假設(shè)隱藏的馬爾可夫鏈在任意時(shí)刻t的狀態(tài)只依賴于其前一時(shí)刻的狀態(tài)，與其他時(shí)刻的狀態(tài)及觀測無關(guān)。

(2)觀測獨(dú)立性假設(shè)，即假設(shè)任意時(shí)刻的觀測只依賴于該時(shí)刻的馬爾可夫鏈的狀態(tài)，與其他觀測及狀態(tài)無關(guān)。

2 觀測序列的生成過程

根據(jù)隱馬爾可夫模型定義，將一個(gè)長度為T的觀測序列O的生成過程描述如下：

觀測序列的生成算法：

輸入：隱馬爾可夫模型

，觀測序列長度

輸出：觀測序列

(1)按照初始狀態(tài)分布

產(chǎn)生狀態(tài)

(2)令t=1

(3)按照狀態(tài)

的觀測概率分布生成

(4)按照狀態(tài)

的狀態(tài)轉(zhuǎn)移概率分布產(chǎn)生狀態(tài)

令t=t+1；如果t<>

2.1 觀測序列生成Python實(shí)現(xiàn)

算法首先初始化兩個(gè)長度為T的向量，接著按照初始狀態(tài)分布pi生成第一個(gè)狀態(tài)，取出狀態(tài)對應(yīng)的觀測的概率分布，生成一個(gè)觀測。

接下來都是按前一個(gè)狀態(tài)取出狀態(tài)轉(zhuǎn)移概率分布生成狀態(tài)，再取出狀態(tài)對應(yīng)的觀測的概率分布生成一個(gè)觀測。重復(fù)該步驟就得到長度為T的觀測和狀態(tài)向量。

代碼如下：

import numpy as np
class HMM():

   def __init__(self, A, B, pi):
       self.A = A
       self.B = B
       self.pi = pi
   def simulate(self, T):
       # draw_from接受一個(gè)概率分布，然后生成該分布下的一個(gè)樣本。    
       def draw_from(probs):
           return np.where(np.random.multinomial(1, probs) == 1)[0][0]
       observations = np.zeros(T, dtype=int)
       states = np.zeros(T, dtype=int)    
       states[0] = draw_from(self.pi)
       observations[0] = draw_from(self.B[states[0], :])
       for t in range(1, T):
           states[t] = draw_from(self.A[states[t-1], :])
           observations[t] = draw_from(self.B[states[t], :])
       return states, observations


3 隱馬爾可夫模型的3個(gè)基本問題

(1)   概率計(jì)算問題。給定模型

和觀測序列

,計(jì)算在模型

下觀測序列O出現(xiàn)的概率

。

(2)   學(xué)習(xí)問題。己知觀測序列

,估計(jì)模型

參數(shù)，使得在該模型下觀測序列概率

最大。即用極大似然估計(jì)的方法估計(jì)參數(shù)。

(3)   預(yù)測問題，也稱為解碼（decoding)問題。已知模型

和觀測序列

，求對給定觀測序列條件概率

最大的狀態(tài)序列

。即給定觀測序列，求最有可能的對應(yīng)的狀態(tài)序列。

下面各節(jié)將逐一介紹這些基本問題的解法。

4 概率計(jì)算算法

這節(jié)介紹計(jì)算觀測序列概率的前向（forward)與后向（backward)算法，前向后向算法就是求第一個(gè)狀態(tài)的前向概率或最后一個(gè)狀態(tài)的后向概率，然后向后或向前遞推即可。

4.1 前向算法

前向概率給定隱馬爾可夫模型

，定義到時(shí)刻t為止的觀測序列為且狀態(tài)為的概率為前向概率，記作

可以遞推地求得前向概率

及觀測序列概率。

觀測序列概率的前向算法：

輸入：隱馬爾可夫模型

，觀測序列;

輸出：觀測序列概率

。

(1)初值

前向概率的定義中一共限定了兩個(gè)條件，一是到當(dāng)前為止的觀測序列，另一個(gè)是當(dāng)前的狀態(tài)。所以初值的計(jì)算也有兩項(xiàng)（觀測和狀態(tài)），一項(xiàng)是初始狀態(tài)概率，另一項(xiàng)是發(fā)生到當(dāng)前觀測的概率。

(2)遞推對

每次遞推同樣由兩部分構(gòu)成，大括號(hào)中是當(dāng)前狀態(tài)為i且觀測序列的前t個(gè)符合要求的概率，括號(hào)外的是狀態(tài)i發(fā)生觀測t+1的概率。

(3)終止

由于到了時(shí)間T，一共有N種狀態(tài)發(fā)生了最后那個(gè)觀測，所以最終的結(jié)果要將這些概率加起來。

4.2 前向算法Python實(shí)現(xiàn)

代碼如下：

1def forward(self, obs_seq):
2    '''前向算法'''
3    N = self.A.shape[0]
4    T = len(obs_seq)
5    F = np.zeros((N, T))
6    F[:, 0] = self.pi * self.B[:, obs_seq[0]]
7    for t in range(1, T):
8        for n in range(N):
9            F[n, t] = np.dot(F[:, t-1], (self.A[:, n])) * self.B[n, obs_seq[t]]
10    return F


4.3 后向算法

后向算法相當(dāng)于前向算法的反向遞推，具體可查看書178頁。

4.4 后向算法Python實(shí)現(xiàn)

代碼如下：

1def backward(self, obs_seq):
2    '''后向算法'''
3    N = self.A.shape[0]
4    T = len(obs_seq)
5    M = np.zeros((N, T))
6    M[:, T-1] = 1
7    # 或者M(jìn)[:, -1:] = 1，列上表示最后一行
8    for t in reversed(range(T-1)):
9        for n in range(N):
10            M[n, t] = np.dot(self.A[n, :], M[:, t+1]) * self.B[n, obs_seq[t+1]]
11    return M


5 學(xué)習(xí)問題

Baum-Welch算法也就是EM算法，己知觀測序列

,估計(jì)模型參數(shù)，使得在該模型下觀測序列概率最大。即用極大似然估計(jì)的方法估計(jì)參數(shù)。而如何解決這個(gè)估計(jì)問題這里采用Baum-Welch算法。 

關(guān)于算法的詳細(xì)解釋可參考：https://blog.csdn.net/u014688145/article/details/53046765?locationNum=7&fps=1

5.1 Baum-Welch算法Python實(shí)現(xiàn)

代碼如下：

1def EM(self, observation, criterion=0.05):
2    '''EM算法進(jìn)行參數(shù)學(xué)習(xí)'''
3    n_state = self.A.shape[0]
4    n_sample = len(observation)
5    done = 1
6    while done:
7        Alpha = self.forward(observation)
8        Beta = self.backward(observation)
9        xi = np.zeros((n_state, n_state, n_sample-1))
10        gamma = np.zeros((n_state, n_sample))
11        for t in range(n_sample-1):
12            denom = np.sum(np.dot(Alpha[:, t].T, self.A) * self.B[:, observation[t+1]].T * Beta[:, t+1].T)
13            sum_gamma1 = np.sum(Alpha[:, t] * Beta[:, t])
14            for i in range(n_state):
15                numer = Alpha[i, t] * self.A[i, :] * self.B[:, observation[t+1]].T * Beta[:, t+1].T
16                xi[i, :, t] = numer/denom
17            gamma[i, t] = Alpha[i, t] * Beta[i, t] / sum_gamma1
18        last_col = Alpha[:, n_sample-1].T * Beta[:, n_sample-1]
19        gamma[:, n_sample-1] = last_col / np.sum(last_col)
20        # 更新參數(shù)
21        n_pi = gamma[:, 0]
22        n_A = np.sum(xi, 2) / np.sum(gamma[:, :-1], 1)
23        n_B = np.copy(self.B)
24        num_level = self.B.shape[1]
25        sum_gamma = 0
26        a = 0
27        for lev in range(num_level):
28            for h in range(n_state):
29                for t in range(n_sample):
30                    sum_gamma = sum_gamma + gamma[h, t]
31                    if observation[t] == lev:
32                        a = a + gamma[h, t]
33                n_B[h, lev] = a / sum_gamma
34                a = 0
35        # 檢查閾值
36        if np.max(np.abs(self.pi-n_pi)) < criterion="">and np.max(np.abs(self.B-n_B)) < criterion="">
37                and np.max(np.abs(self.A-n_A)) <>criterion:
38            done = 0
39        self.A, self.B, self.pi = n_A, n_B, n_pi
40    return self.A, self.B, self.pi


6 預(yù)測算法

6.1 維特比算法

輸入：模型

和觀測;

輸出：最優(yōu)路徑

。

⑴初始化

(2)遞推。對

(3)終止

(4)最優(yōu)路徑回溯。對

求得最優(yōu)路徑

。

6.2 維特比算法Python實(shí)現(xiàn)

維特比算法是一種動(dòng)態(tài)規(guī)劃算法用于尋找最有可能產(chǎn)生觀測事件序列的-維特比路徑-隱含狀態(tài)序列。維特比算法其實(shí)就是多步驟每步多選擇模型的最優(yōu)選擇問題，其在每一步的所有選擇都保存了前續(xù)所有步驟到當(dāng)前步驟當(dāng)前選擇的最小總代價(jià)（或者最大價(jià)值）以及當(dāng)前代價(jià)的情況下前繼步驟的選擇。依次計(jì)算完所有步驟后，通過回溯的方法找到最優(yōu)選擇路徑。

代碼如下：

1def viterbi(self, obs_seq):
2    '''viterbi算法預(yù)測狀態(tài)序列'''
3    N = self.A.shape[0]
4    T = len(obs_seq)
5    P = np.zeros((N, T))
6    prev_point = np.zeros((N, T))
7    prev_point[:, 0] = 0
8    P[:, 0] = self.pi * self.B[:, obs_seq[0]]
9    for t in range(1, T):
10        for n in range(N):
11            P[n, t] = np.max(P[:, t - 1] * self.A[:, n]) * self.B[n, obs_seq[t]]
12            prev_point[n, t] = np.argmax(P[:, t - 1] * self.A[:, n] * self.B[n, obs_seq[t]])
13    return P, prev_point


6.3 最優(yōu)路徑生成Python實(shí)現(xiàn)

最優(yōu)路徑其實(shí)也是維特比算法的一部分，當(dāng)已經(jīng)確定了T時(shí)刻的最優(yōu)狀態(tài)i，接下來通過回溯的方式找到最優(yōu)路徑。

代碼如下：

1def build_path(self, obs_seq):
2    '''return the optimal path'''
3    P, prev_point = self.viterbi(obs_seq)
4    T = len(obs_seq)
5    opt_path = np.zeros(T)
6    last_state = np.argmax(P[:, T-1])
7    opt_path[T-1] = last_state
8    for t in reversed(range(T-1)):
9        opt_path[t] = prev_point[int(opt_path[t+1]), t+1]
10    last_path = reversed(opt_path)
11    return last_path


6.4 實(shí)踐

用《統(tǒng)計(jì)學(xué)習(xí)方法》中的案例數(shù)據(jù)進(jìn)行測試，代碼如下：

1if __name__ == '__main__':
2    # 用《統(tǒng)計(jì)學(xué)習(xí)方法》中的案例進(jìn)行測試
3    A = np.array([[0.5, 0.2, 0.3], [0.3, 0.5, 0.2], [0.2, 0.3, 0.5]])
4    B = np.array([[0.5, 0.5], [0.4, 0.6], [0.7, 0.3]])
5    pi = np.array([0.2, 0.4, 0.4])
6    test1 = HMM(A, B, pi)
7    obs_seq = [0, 1, 0]
8    print(test1.forward(obs_seq))
9    print(test1.backward(obs_seq))
10    print(test1.viterbi(obs_seq))
11    print(test1.build_path(obs_seq))
12    print(test1.EM(obs_seq))


參考：

1. 李航. 統(tǒng)計(jì)學(xué)習(xí)方法[M]. 北京：清華大學(xué)出版社，2012

       碼農(nóng)場.隱馬爾可夫模型：http://www.hankcs.com/ml/hidden-markov-model.html