隱馬爾科夫模型（HMM）學習筆記二

時間 2019-12-08

標籤模型 hmm 學習筆記简体版

原文原文鏈接

　　這裏接着學習筆記一中的問題2，說實話問題2中的Baum-Welch算法編程時矩陣轉換有點燒腦，開始編寫一直不對（編程還不熟練hh），後面在紙上仔細推了一遍，由特例慢慢改寫才運行成功，因此代碼裏面好多處都有print。html

　　筆記一中對於問題1（機率計算問題）採用了前向或後向算法，根據前向和後向算法能夠獲得一些後面要用到的機率與指望值。算法

1、問題2 學習問題已知觀測序列，估計模型參數，使得在該模型下觀測序列機率最大

　　隱馬爾可夫模型的學習，根據訓練數據除包括觀測序列O外是否包括了對應的狀態序列 I 分爲監督學習和非監督學習。若是包括了狀態序列 I ，則能夠直接採用極大似然估計來估計隱馬爾可夫模型的初始狀態機率Pi，狀態轉移機率A和觀測機率B（通常統計頻數）。可是通常沒有對應的狀態序列 I ，若是人工標註訓練數據的話代價過高，因此大多時候採用非監督學習方法------Baum-Welch算法。　編程

　　Baum-Welch算法

　　假設給定訓練數據只包含S個長度爲T的觀測序列而沒有對應的狀態序列，目標是學習隱馬爾可夫模型的參數。咱們將觀測序列數據看做觀測數據O，狀態序列數據看做不可觀測的隱數據I，那麼隱馬爾可夫模型事實上是一個含有隱變量的機率模型app

它的參數學習能夠由EM算法實現（EM算法能夠參考博客，博主寫得通俗易懂）。函數

　　一、肯定徹底數據的對數似然函數學習

　　全部觀測數據寫成,全部隱數據寫成，徹底數據是。徹底數據的對數似然函數是。測試

　　二、EM算法的E步：spa

　　求Q函數.net

其中，是隱馬爾可夫模型參數的當前估計值，是要極大化的隱馬爾可夫模型參數。（Q函數的標準定義是：，式子內部實際上是條件機率，其中的對應；其與無關，因此省略掉了。）code

因而函數能夠寫成：

式中求和都是對全部訓練數據的序列總長度T進行的。這個式子是將代入後（取對數後變成加法，便於求解），將初始機率、轉移機率、觀測機率這三部分乘積的對數拆分爲對數之和，因此有三項。

　　三、EM算法的M步:極大化Q函數求模型參數，因爲要極大化的參數在Q函數表達式中單獨地出如今3個項中，因此只需對各項分別極大化。

　　第1項能夠寫成：

注意到知足約束條件利用拉格朗日乘子法，寫出拉格朗日函數：

對其求偏導數並令結果爲0

獲得

這個求導是很簡單的，求和項中非i的項對πi求導都是0，logπ的導數是1/π，γ那邊求導就剩下πi本身對本身求導，也就是γ。等式兩邊同時乘以πi就獲得了上式。

對i求和獲得γ：

代入中獲得：

第2項能夠寫成：

相似第1項，應用具備約束條件的拉格朗日乘子法能夠求出

第3項爲:

一樣用拉格朗日乘子法，約束條件是。注意，只有在對時對的偏導數纔不爲0,以表示。求得

　　Baum-Welch模型參數估計公式

　　將這三個式子中的各機率分別簡寫以下：

　　則可將相應的公式寫成：

這三個表達式就是Baum-Welch算法（Baum-Welch algorithm)，它是EM算法在隱馬爾可夫模型學習中的具體實現，由Baum和Welch提出。

　　算法 (Baum-Welch算法）

　　輸入：觀測數據

　　輸出：隱馬爾可夫模型參數。

　　(1)初始化。對,選取，獲得模型。

　　(2)遞推。對

右端各值按觀測和模型計算。

　　(3)終止。獲得模型參數。

　　分析：先根據前向算法得出alpha，根據後向算法獲得beta，以後再根據本文最開始的公式計算gamma和sigamma矩陣，這裏注意gamma是N*N維的，而sigamma是T-1*N*N維的，最後執行Baum-Welch算法。代碼裏面有相應的註釋，這裏我也卡了很久，能夠先本身到紙上推一下。

def Baum_Welch(pi, A, B, obs_seq, error=0.005):
    switch = 0  # 判斷是否中止迭代
    if not switch:
        N = A.shape[0]  # 可能的狀態個數
        M = B.shape[1]  # 可能的觀測結果個數
        T = len(obs_seq)  # 觀測序列長度
        newB = np.zeros((N, M)) # 初始化觀測機率
        alpha = forward_compute(pi, A, B, obs_seq)[1]   # 前向算法獲得的alpha矩陣  N*T維
        print('alpha:', alpha)
        beta = backword_compute(pi, A, B, obs_seq)[1]   # 後向算法獲得的beta矩陣  N*T維
        print('beta:', beta)
        gamma = np.zeros((N, T))  # gamma_t_i 表示t時刻在狀態q_i的機率 N*N維
        sigamma = np.zeros((T-1, N, N))  # sigamma_t_i_j t時刻處於q_i,t+1時刻處於q_j的機率  (T-1)*N*N維
 
        for t in range(T-1):
            gamma[:, t] = (alpha[:, t]*beta[:, t]) / (alpha[:, t].dot(beta[:, t]))  # 求出gamma矩陣
            denom = np.dot(np.dot(alpha[:,t].T, A) * B[:,obs_seq[t+1]].T, beta[:,t+1])  # sigmma矩陣分母
            for i in range(N):
                molecule = alpha[i, t]* A[i,:] * B[:,obs_seq[t+1]]*beta[:,t+1]  # 分子
                sigamma[t,i,:] = molecule / denom   # 求sigamma
        
        gamma[:, T-1] = (alpha[:, T-1]*beta[:, T-1]) / (alpha[:, T-1].dot(beta[:, T-1]))  # 因爲sigamma直到T-1,因此gamma最後要添加一列
        print('gamma:', gamma)
        print('sigamma:', sigamma)
        # print('-----------')
        # print(np.sum(sigamma, axis=0))
        # print('*********')
        # print(np.sum(gamma[:,:T-1], axis=1))
        # print(np.sum(gamma[:,:T-1], axis=1).reshape(N, -1))
        #更新
        newA = np.sum(sigamma, axis=0)/(np.sum(gamma[:,:T-1], axis=1).reshape(N, -1))
        # print(newA)
        # print(obs_seq==0)
        # print(np.sum(gamma, axis=1)
        for m in range(M):
            newB[:,m] = np.sum(gamma.T[obs_seq==m], axis=0)/(np.sum(gamma, axis=1))
        # print(newB)
        newpi = gamma[:, 0] 
    
        # 檢查是否知足要求
        if np.max(abs(pi - newpi)) < error and np.max(abs(A - newA)) < error and np.max(abs(B - newB)) < error:
                switch = 1
        pi, A, B = newpi, newA, newB
    return A, B, pi

　　以後帶入相關數據進行測試。

A = np.array([[0.1, 0.3, 0.6],[0.3, 0.5, 0.2], [0.3, 0.5, 0.2]])
B = np.array([[0.5, 0.5],[0.5, 0.5], [0.5,0.5]])
pi = np.array([0.2, 0.4, 0.4])

observations_data = np.array([0, 1, 0, 0, 1])

newA, newB, newpi = Baum_Welch(pi, A, B, observations_data)
print("newA: ", newA)
print("newB: ", newB)
print("newpi: ", newpi)

2、問題3 預測問題已知模型和觀測序列，求給定序列最大機率下的狀態

　　預測算法有兩種：近似算法和維特比算法。

　　1.近似算法

　　基本思想：在每一個時刻t選擇在該時刻最有可能出現的狀態 i^*_t ，從而獲得一個狀態序列 I^* = (i^*₁, i^*₂, ..., i^*_T)。

　　對於給定隱馬爾可夫模型和觀測序列O，在 t 時刻處於狀態 qi 的機率gamma(i) 爲：

在每一時刻 t 最有可能的狀態 i*t 爲：

從而獲得狀態序列 I^* = (i^*₁, i^*₂, ..., i^*_T)。

這個代碼很簡單，直接對gamma按行求最大值就行。

'''
gamma: [[0.2     0.26    0.248   0.2504  0.44992]
        [0.4     0.46    0.448   0.2992  0.24992]
        [0.4     0.28    0.304   0.4504  0.30016]]
'''
def envolution_state(gamma):
    states = []
    for i in range(gamma.shape[1]):
        max_index = np.where(gamma==np.max(gamma[:,i]))[0][0]
        states.append(max_index)
    return states

print('*********test*********')
gamma = np.array([[0.2 ,0.26,0.248,0.2504,0.44992], [0.4,0.46,0.448,0.2992,0.24992], [0.4,0.28,0.304,0.4504,0.30016]])
print('state:', envolution_state(gamma))

　　二、維特比算法

　　基本思想：利用動態規劃求解隱馬爾科夫模型預測問題，即用動態規劃求解機率最大路徑（最優路徑），一條路徑對應一個狀態序列。

　　算法(維特比算法）

　　輸入：模型和觀測;

　　輸出：最優路徑。

　　⑴初始化

　　(2)遞推。對

　　(3)終止

　　(4)最優路徑回溯。對

求得最優路徑。

def vetebi(pi, A, B, obs_seq):
    N = A.shape[0]  # 狀態個數
    T = len(obs_seq) # 觀測序列長度
    
    deta = np.zeros((N, T))  # 初始化deta矩陣
    fia = np.zeros((N, T), dtype=int)  # 初始化fia矩陣
    path = np.zeros(T, dtype=int)  # 初始化最優路徑
    
    deta[:, 0] = pi*B[:, obs_seq[0]]  # 計算deta_1
    for t in range(1, T):
        matrix = (deta[:, t-1]*A.T)*B[:, obs_seq[t]].reshape(N, -1)  # 計算乘法，deta按列取其中最大的
        deta[:, t] = np.max(matrix, axis=1)  # 給deta賦值
        fia[:, t] = np.argmax(matrix, axis=1)  # 給fia賦值，用於求最優路徑
    print(deta)
    print(fia)
    
    # 最大機率
    P_max = np.max(deta[:, T-1]) 
    # print(P_max)
    
    # 求最優路徑
    path[T-1] = np.max(fia[:, T-1])
    for t in range(T-2, -1, -1):
        path[t] = fia[:, t+1][path[t+1]]
    # print(path)
    return path+1

　　利用書本上例題進行測試。

# test
A = np.array([[0.5, 0.2, 0.3], [0.3, 0.5, 0.2], [0.2, 0.3, 0.5]])
B = np.array([[0.5, 0.5], [0.4, 0.6], [0.7, 0.3]])
pi = np.array([0.2, 0.4, 0.4])
obs_seq = np.array([0, 1, 0])
path = vetebi(pi, A, B, obs_seq)
print('the best path is :',path)