隱馬爾可夫模型（HMM）學習筆記一

        學習了李航的《統計學習方法》中隱馬爾可夫模型(Hidden Markov Model, HMM)，這裏把本身對HMM的理解進行總結(大部分是書本原文，O(∩_∩)O哈哈~，主要是想利用python將其實現一遍，這樣印象深入一點兒)，並利用python將書本上的例子運行一遍。HMM是可用於標註問題的統計學習模型，描述由隱藏的馬爾科夫鏈隨機生成觀測序列的過程，屬於生成模型。HMM在語音識別、天然語言處理、生物信息、模式識別等領域都有着普遍的應用。

 

1、隱馬爾可夫模型的定義

       HMM由初始機率分佈、狀態轉移機率分佈A以及觀測機率分佈B三部分肯定。HMM模型能夠用三元符號表示，即：。A、B、稱爲HMM的三要素。

       設Q是全部可能的狀態的集合，V是全部可能的觀測的集合。其中，N是可能的狀態數，M是可能的觀測數。這裏觀測狀態v是可見的，狀態q是不可見的。記I是長度爲T的狀態序列，O是對應的觀測序列，。

　　一、則可定義狀態轉移機率矩陣B以下：

其中是在時刻t處於狀態 q_i 的條件下在時刻t+1轉移到狀態 q_j 的機率。

　　二、觀測機率矩陣B定義以下：

其中是在時刻t處於狀態 q_j 的條件下生成觀測 v_k 的機率。

　　三、初始狀態機率向量以下：

其中是時刻t=1處於狀態qj的機率。

 

2、HMM的兩個假設

　　一、齊次馬爾科夫性假設。即隱藏的馬爾科夫鏈在任意時刻 t 的狀態只與前一時刻的狀態有關，與時刻 t 也無關。

。

　　二、觀測獨立性假設。即任意時刻的觀測只和該時刻的馬爾科夫鏈的狀態有關，與其餘觀測即狀態無關。

。

 

例題：例子採用書本上的盒子和球模型。

　　分析：這裏一共四個盒子，則表示全部可能的狀態有四種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、隱馬爾可夫模型的3個基本問題

　　 一、機率計算問題。給定模型和觀測序列,計算在模型下觀測序列O出現的機率。

　　二、學習問題。己知觀測序列,估計模型參數，使得在該模型下觀測序列機率最大。即用極大似然估計的方法估計參數。

　　三、預測問題，也稱爲解碼（decoding）問題。已知模型和觀測序列，求對給定觀測序列條件機率最大的狀態序列。即給定觀測序列，求最有可能的對應的狀態序列。

 

　　問題1 機率計算問題

　　若是直接採用直接計算法的話，須要列舉全部長度爲T的狀態序列 I，而後求出 I 與觀測序列O的聯合機率，最後再對全部可能的狀態序列求和，獲得序列O在HMM模型下生成的機率。這樣時間複雜度是O(TN^T)階的，這種算法觀測序列長度很長的話計算量太大，因此採用前向、後向算法計算機率。

 

　 （1）前向算法

 　　首先定義(前向機率）給定隱馬爾可夫模型，定義到時刻t爲止的觀測序列爲且狀態爲的機率爲前向機率，記做

　　輸入：隱馬爾可夫模型，觀測序列;

　　輸出：觀測序列機率。

　　(1)初值

　　(2)遞推對 

　　(3)終止

　　因爲到了時間T，一共有N種狀態發轉移到最後那個觀測，因此最終的結果要將這些機率加起來。而每次遞推都是在前一次的基礎上作的，因此只需累加O(T)次，因此整體複雜度是O(N²T)。

 

　　分析：前向機率矩陣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)遞推。對

 

 

a_ij表示狀態i轉移到j的機率，b_j表示發射O_t+1，表示j後面的序列對應的後向機率。

 

　　(3)終止

 

 

最後的求和是由於，在第一個時間點上有N種後向機率都能輸出從2到T的觀測序列，因此乘上輸出O₁的機率後求和獲得最終結果。

 

　　分析：前向機率矩陣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)