學習了李航的《統計學習方法》中隱馬爾可夫模型(Hidden Markov Model, HMM)，這里把自己對HMM的理解進行總結(大部分是書本原文，O(∩_∩)O哈哈~，主要是想利用python將其實現一遍，這樣印象深刻一點兒)，並利用python將書本上的例子運行一遍。HMM是可用於標注問題的統計學習模型，描述由隱藏的馬爾科夫鏈隨機生成觀測序列的過程，屬於生成模型。HMM在語音識別、自然語言處理、生物信息、模式識別等領域都有着廣泛的應用。
一、隱馬爾可夫模型的定義
HMM由初始概率分布
、狀態轉移概率分布A以及觀測概率分布B三部分確定。HMM模型
可以用三元符號表示，即：
。A、B、稱為HMM的三要素。
設Q是所有可能的狀態的集合，V是所有可能的觀測的集合。
其中，N是可能的狀態數，M是可能的觀測數。這里觀測狀態v是可見的，狀態q是不可見的。記I是長度為T的狀態序列，O是對應的觀測序列，
。
1、則可定義狀態轉移概率矩陣B如下：
其中
是在時刻t處於狀態 qi 的條件下在時刻t+1轉移到狀態 qj 的概率。
2、觀測概率矩陣B定義如下：
其中
是在時刻t處於狀態 qj 的條件下生成觀測 vk 的概率。
3、初始狀態概率向量如下：
其中
是時刻t=1處於狀態qj的概率。
二、HMM的兩個假設
1、齊次馬爾科夫性假設。即隱藏的馬爾科夫鏈在任意時刻 t 的狀態只與前一時刻的狀態有關，與時刻 t 也無關。
。
2、觀測獨立性假設。即任意時刻的觀測只和該時刻的馬爾科夫鏈的狀態有關，與其他觀測即狀態無關。
。
例題：例子采用書本上的盒子和球模型。
分析：這里一共四個盒子，則表示所有可能的狀態有四種Q={盒子1，盒子2，盒子3，盒子4}，N=4。球的顏色只有紅白兩種，即所有可能的觀測為V={紅，白}，M=2。
用python中的列表分別表示可能的狀態集合Q和觀測集合V，用字典存儲初始概率分布、狀態轉移概率分布和觀測概率分布。
# 所有可能的狀態集合 states = ['盒子1', '盒子2', '盒子3', '盒子4'] # 所有可能的觀測集合 observations = ['紅球', '白球'] # 初始狀態概率分布pi Pi = {'盒子1':0.25, '盒子2':0.25, '盒子3':0.25, '盒子4':0.25} # 狀態轉移概率A A = {'盒子1':{'盒子1':0, '盒子2':1, '盒子3':0, '盒子4':0}, '盒子2':{'盒子1':0.4, '盒子2':0, '盒子3':0.6, '盒子4':0}, '盒子3':{'盒子1':0, '盒子2':0.4, '盒子3':0, '盒子4':0.6}, '盒子4':{'盒子1':0, '盒子2':0, '盒子3':0.5, '盒子4':0.5} } # 觀測概率分布B B = {'盒子1':{'紅球':0.5, '白球':0.5}, '盒子2':{'紅球':0.3, '白球':0.7}, '盒子3':{'紅球':0.6, '白球':0.4}, '盒子4':{'紅球':0.8, '白球':0.2} }
為了后面方便處理，這里引入python的numpy庫將Pi、A和B轉換為矩陣和向量形式。
import numpy as np # 將初始狀態概率轉換成向量形式 Pi = np.array([Pi[x] for x in Pi]) print('初始狀態pi:', Pi) # 將初始概率分布轉換為矩陣形式 trans_A = np.array([A[x][y] for x in A for y in A[x]]) A = trans_A.reshape((len(A), -1)) print('狀態轉移概率分步A:', A) # 將觀測概率分布轉換成矩陣形式 trans_B = np.array([B[x][y] for x in B for y in B[x]]) B = trans_B.reshape((len(B), -1)) print('觀測概率分布B:', B)
觀測序列的生成過程
根據隱馬爾可夫模型定義，可以將一個長度為T的觀測序列
的生成過程描述如下：
算法(觀測序列的生成）
輸入：隱馬爾可夫模型
，觀測序列長度
輸出：觀測序列。
(1)按照初始狀態分布
產生狀態
(2)令t=1
(3)按照狀態的觀測概率分布
生成
(4)按照狀態的狀態轉移概率分布
產生狀態
令
；如果
則轉步(3)；否則，終止。
分析：對於按照某個概率分布生成狀態，這里采用python中的np.random.multinomial函數按照多項式分布，生成數據，然后再利用np.where得到滿足條件的元素對應的下標。定義觀測序列生成函數generation_observation。
# 定義生成觀測序列的函數 def generation_observation(pi, A, B, T): def random_res(probs): """ 1.np.random.multinomial: 按照多項式分布，生成數據 >>> np.random.multinomial(20, [1/6.]*6, size=2) array([[3, 4, 3, 3, 4, 3], [2, 4, 3, 4, 0, 7]]) For the first run, we threw 3 times 1, 4 times 2, etc. For the second, we threw 2 times 1, 4 times 2, etc. 2.np.where: >>> x = np.arange(9.).reshape(3, 3) >>> np.where( x > 5 ) (array([2, 2, 2]), array([0, 1, 2])) """ return np.where(np.random.multinomial(1,probs) == 1)[0][0] observation_data = np.zeros(T, dtype=int) # 初始化生成的觀測序列 observation_states = np.zeros(T, dtype=int) # 觀測序列對應的狀態 observation_states[0] = random_res(pi) # 根據初始狀態分布產生狀態1 observation_data[0] = random_res(B[observation_states[0]]) # 由狀態1生成對應的觀測值 for i in range(1, T): observation_states[i] = random_res(A[observation_states[i-1]]) # 由前一個狀態轉變到下一個狀態 observation_data[i] = random_res(B[observation_states[i]]) return observation_data, observation_states observation_data, observation_states = generation_observation(Pi, A, B, 10) print(observation_data, observation_states) print('觀測序列:', [observations[i] for i in observation_data]) print('狀態序列:', [states[i] for i in observation_states])
三、隱馬爾可夫模型的3個基本問題
1、概率計算問題。給定模型和觀測序列,計算在模型
下觀測序列O出現的概率
。
2、學習問題。己知觀測序列,估計模型參數，使得在該模型下觀測序列概率
最大。即用極大似然估計的方法估計參數。
3、預測問題，也稱為解碼（decoding）問題。已知模型和觀測序列，求對給定觀測序列條件概率
最大的狀態序列
。即給定觀測序列，求最有可能的對應的狀態序列。
問題1 概率計算問題
如果直接采用直接計算法的話，需要列舉所有長度為T的狀態序列 I，然后求出 I 與觀測序列O的聯合概率，最后再對所有可能的狀態序列求和，得到序列O在HMM模型下生成的概率。這樣時間復雜度是O(TNT)階的，這種算法觀測序列長度很長的話計算量太大，所以采用前向、后向算法計算概率。
（1）前向算法
首先定義(前向概率）給定隱馬爾可夫模型
，定義到時刻t為止的觀測序列為
且狀態為
的概率為前向概率，記作
輸入：隱馬爾可夫模型，觀測序列
;
輸出：觀測序列概率
。
(1)初值
(2)遞推對
(3)終止
由於到了時間T，一共有N種狀態發轉移到最后那個觀測，所以最終的結果要將這些概率加起來。而每次遞推都是在前一次的基礎上做的，所以只需累加O(T)次，所以總體復雜度是O(N2T)。
分析：前向概率矩陣alpha是N*T維的，可以先用零初始化，之后再按照前向算法一步一步對其值進行替換。定義前向算法方法forward_compute。再利用例題10.2對其進行驗證。
def forward_compute(pi, A, B, o_seq): T = len(o_seq) N = A.shape[0] alpha = np.zeros((N, T)) # 初始化alpha矩陣 N*T維 # print(alpha) # print(B[:, o_seq[0]]) alpha[:, 0] = pi*B[:, o_seq[0]] # step1 初值 for t in range(1, T): alpha[:, t] = alpha[:, t-1].dot(A) * B[:, o_seq[t]] # step2 遞推 # print(alpha) return np.sum(alpha[:, T-1]), alpha # step3 終止 對時刻T的alpha求和
# 測試 例題10.2 A = np.array([[0.5, 0.2, 0.3], [0.3, 0.5, 0.2], [0.2, 0.3, 0.5]]) Pi = np.array([0.2, 0.4, 0.4]) B = np.array([[0.5, 0.5], [0.4, 0.6], [0.7, 0.3]]) o_seq = [0, 1, 0] # 紅、白、紅 prob, alpha = forward_compute(Pi, A, B, o_seq) print('觀測概率為：', prob) print('alpha:', alpha)
（2）后向算法
定義后向概率為給定隱馬爾可夫模型,定義在時刻t狀態為的條件下，從t+1到T的部分觀測序列為
的概率為后向概率，記作
可以用遞推的方法求得后向概率
及觀測序列概率
。
算法(觀測序列概率的后向算法）
輸入：隱馬爾可夫模型,觀測序列:
輸出：觀測序列概率
(1)初值
根據定義，從T+1到T的部分觀測序列其實不存在，所以硬性規定這個值是1。
(2)遞推。對
aij表示狀態i轉移到j的概率，bj表示發射Ot+1，
表示j后面的序列對應的后向概率。
(3)終止
最后的求和是因為，在第一個時間點上有N種后向概率都能輸出從2到T的觀測序列，所以乘上輸出O1的概率后求和得到最終結果。
分析：前向概率矩陣beta也是N*T維的，可以先用零初始化，之后再按照前向算法一步一步對其值進行替換。定義前向算法方法backward_compute。再利用例題10.2對其進行驗證。
def backword_compute(pi, A, B, o_seq): T = len(o_seq) beta = np.zeros((A.shape[0], T)) # 初始化beta矩陣 beta[:, T-1] = np.ones(A.shape[0]) # step1 初值 for t in range(T-2, -1, -1): beta[:, t] = np.sum(A*B[:, o_seq[t+1]]*beta[:, t+1], axis=1) # step2 遞推 res_prop = np.sum(pi*B[:, o_seq[0]]*beta[:,0], axis=0) # 求觀測序列概率 # print(beta) return res_prop, beta
# 測試 例題10.2 A = np.array([[0.5, 0.2, 0.3], [0.3, 0.5, 0.2], [0.2, 0.3, 0.5]]) Pi = np.array([0.2, 0.4, 0.4]) B = np.array([[0.5, 0.5], [0.4, 0.6], [0.7, 0.3]]) o_seq = [0, 1, 0] # 紅、白、紅 prob, beta = backword_compute(Pi, A, B, o_seq) print('觀測概率：', prob) print('beta:', beta)
