Python實現HMM（隱馬爾可夫模型）

標籤：sha   文字   多個   細節   次數   圖結構   簡單的   代碼   test   

1. 前言

隱馬爾科夫HMM模型是一類重要的機器學習方法，其主要用於序列資料的分析，廣泛應用於語音辨識、文本翻譯、序列預測、中文分詞等多個領域。雖然近年來，由於RNN等深度學習方法的發展，HMM模型逐漸層得不怎麼流行了，但並不意味著完全退出應用領域，甚至在一些輕量級的任務中仍有應用。本系列部落格將詳細剖析隱馬爾科夫鏈HMM模型，同以往網路上絕大多數教程不同，本系列部落格將更深入地分析HMM，不僅包括估計序列隱狀態的維特比演算法（HMM解碼問題）、前向後向演算法等，而且還著重的分析HMM的EM訓練過程，並將所有的過程都通過數學公式進行推導。

由於隱馬爾科夫HMM模型是一類非常複雜的模型，其中包含了大量機率統計的數學知識，因此網路上多數部落格一般都採用舉例等比較通俗的方式來介紹HMM，這麼做會讓初學者很快明白HMM的原理，但要丟失了大量細節，讓初學者處於一種似懂非懂的狀態。而本文並沒有考慮用非常通俗的文字描述HMM，還是考慮通過詳細的數學公式來一步步引導初學者掌握HMM的思想。另外，本文重點分析了HMM的EM訓練過程，這是網路上其他教程所沒有的，而個人認為相比於維特比演算法、前向後向演算法，HMM的EM訓練過程雖然更為複雜，但是一旦掌握這個訓練過程，那麼對於通用的鏈狀圖結構的推導、EM演算法和網路訓練的理解都會非常大的協助。另外通過總結HMM的數學原理，也能非常方便將數學公式改寫成代碼。

最後，本文提供了一個簡單版本的隱馬爾科夫鏈HMM的Python代碼，包含了高斯模型的HMM和離散HMM兩種情況，代碼中包含了HMM的訓練、預測、解碼等全部過程，核心代碼總共只有200~300行代碼，非常簡單！個人代碼水平比較渣=_=||，大家按照我的教程，應該都可以寫出更魯棒性更有高效的代碼，附上Github地址：https://github.com/tostq/Easy_HMM

覺得好，就點星哦！

覺得好，就點星哦！

覺得好，就點星哦！

重要的事要說三遍！！！！

為了方便大家學習，我將整個HMM代碼完善成整個學習項目，其中包括

hmm.py：HMM核心代碼，包含了一個HMM基類，一個高斯HMM模型及一個離散HMM模型

DiscreteHMM_test.py及GaussianHMM_test.py：利用unnitest來測試我們的HMM，同時引入了一個經典HMM庫hmmlearn作為對照組

Dice_01.py：利用離散HMM模型來解決丟色子問題的例子

Wordseg_02.py：解決中文分詞問題的例子

Stock_03.py：解決股票預測問題的例子

 

2. 隱馬爾科夫鏈HMM的背景

首先，已知一組序列 ：

我們從這組序列中推匯出產生這組序列的函數，假設函數參數為 ，其表示為

即使得序列X發生機率最大的函數參數，要解決上式，最簡單的考慮是將序列 的每個資料都視為獨立的，比如建立一個神經網路。然後這種考慮會隨著序列增長，而導致參數爆炸式增長。因此可以假設當前序列資料只與其前一資料值相關，即所謂的一階馬爾科夫鏈：

有一階馬爾科夫鏈，也會有二階馬爾科夫鏈（即當前資料值取決於其前兩個資料值）

當前本文不對二階馬爾科夫鏈進行深入分析了，著重考慮一階馬爾科夫鏈，現在根據一階馬爾科夫鏈的假設，我們有：

因此要解一階馬爾科夫鏈，其關鍵在於求資料（以下稱觀測值）之間轉換函式 ，如果假設轉換函式同序列中位置 （時間）無關，我們就能根據轉換函式而求出整個序列的機率：

然而，如果觀測值x的狀態非常多（特別極端的情況是連續資料），轉換函式會變成一個非常大的矩陣，如果x的狀態有K個，那麼轉換函式就會是一個K*(K-1)個參數，而且對於連續變數觀測值更是困難。

為了降低馬爾科夫鏈的轉換函式的參數量，我們引入了一個包含較少狀態的隱狀態值，將觀測值的馬爾科夫鏈轉換為隱狀態的馬爾科夫鏈（即為隱馬爾科夫鏈HMM）

其包含了一個重要假設：當前觀測值只由當前隱狀態所決定。這麼做的一個重要好處是，隱狀態值的狀態遠小於觀測值的狀態，因此隱藏狀態的轉換函式 的參數更少。

此時我們要決定的問題是：

即在所有可能隱藏狀態序列情況下，求使得序列 發生機率最大的函數參數 。

這裡我們再總結下：

隱馬爾科夫鏈HMM三個重要假設：

1. 當前觀測值只由當前隱藏狀態確定，而與其他隱藏狀態或觀測值無關（隱藏狀態假設）

2. 當前隱藏狀態由其前一個隱藏狀態決定（一階馬爾科夫假設）

3. 隱藏狀態之間的轉換函式機率不隨時間變化（轉換函式穩定性假設）

隱馬爾科夫鏈HMM所要解決的問題：

在所有可能隱藏狀態序列情況下，求使得當前序列X產生機率最大的函數參數θ。

 

 

代碼

http://blog.csdn.net/sinat_36005594/article/details/69568538

 

前幾天用MATLAB實現了HMM的代碼，這次用python寫了一遍，依據仍然是李航博士的《統計學習方法》

由於第一次用python，所以代碼可能會有許多缺陷，但是所有代碼都用書中的例題進行了測試，結果正確。

這裡想說一下python，在編寫HMM過程中參看了之前寫的MATLAB程式，發現他們有太多相似的地方，用到了numpy庫，在python代碼過程中最讓我頭疼的是數組角標，和MATLAB矩陣角標從1開始不同，numpy庫數組角標都是從0開始，而且數組的維數也需要謹慎，一不小心就會出現too many indices for array的錯誤。程式中最後是維特比演算法，在運行過程中出現了__main__:190: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future的警告，還沒有去掉這個警告，查了一下說不影響結果，後面會去解決這個問題，下面貼出My Code

# -*- coding: utf-8 -*-"""Created on Thu Feb 16 19:28:39 20172017-4-2    ForwardBackwardAlg函數功能：實現前向演算法    理論依據：李航《統計學習方法》2017-4-5    修改了ForwardBackwardAlg函數名稱為ForwardAlgo以及輸出的alpha數組形式    完成了BackwardAlgo函數功能：後向演算法    以及函數FBAlgoAppli：計算在觀測序列和模型參數確定的情況下，    某一個隱含狀態對應相應的觀測狀態的機率2017-4-6    完成BaumWelchAlgo函數一次迭代2017-4-7    實現維特比演算法@author: sgp"""import numpy as np#輸入格式如下：#A = np.array([[.5,.2,.3],[.3,.5,.2],[.2,.3,.5]])#B = np.array([[.5,.5],[.4,.6],[.7,.3]])#Pi = np.array([[.2,.4,.4]])#O = np.array([[1,2,1]])#應用ndarray在數組之間進行相互運算時，一定要確保數組維數相同！#比如：#In[93]:m = np.array([1,2,3,4])#In[94]:m#Out[94]: array([1, 2, 3, 4])#In[95]:m.shape#Out[95]: (4,)#這裡表示的是一維數組#In[96]:m = np.array([[1,2,3,4]])#In[97]:m#Out[97]: array([[1, 2, 3, 4]])#In[98]:m.shape#Out[98]: (1, 4)#而這裡表示的就是二維數組#注意In[93]和In[96]的區別,多一對中括弧！！#N = A.shape[0]為數組A的行數， H = O.shape[1]為數組O的列數#在下列各函數中，alpha數組和beta數組均為N*H二維數組，也就是橫向座標是時間，縱向是狀態def ForwardAlgo(A,B,Pi,O):    N = A.shape[0]#數組A的行數    M = A.shape[1]#數組A的列數    H = O.shape[1]#數組O的列數    sum_alpha_1 = np.zeros((M,N))    alpha = np.zeros((N,H))    r = np.zeros((1,N))    alpha_1 = np.multiply(Pi[0,:], B[:,O[0,0]-1])        alpha[:,0] = np.array(alpha_1).reshape(1,N)#alpha_1是一維數組，在使用np.multiply的時候需要升級到二維數組。#錯誤是IndexError: too many indices for array    for h in range(1,H):        for i in range(N):            for j in range(M):                sum_alpha_1[i,j] = alpha[j,h-1] * A[j,i]            r = sum_alpha_1.sum(1).reshape(1,N)#同理，將數組升級為二維數組            alpha[i,h] = r[0,i] * B[i,O[0,h]-1]    #print("alpha矩陣: \n %r" % alpha)        p = alpha.sum(0).reshape(1,H)    P = p[0,H-1]    #print("觀測機率: \n %r" % P)    #return alpha    return alpha, P    def BackwardAlgo(A,B,Pi,O):    N = A.shape[0]#數組A的行數    M = A.shape[1]#數組A的列數    H = O.shape[1]#數組O的列數        #beta = np.zeros((N,H))    sum_beta = np.zeros((1,N))    beta = np.zeros((N,H))    beta[:,H-1] = 1    p_beta = np.zeros((1,N))        for h in range(H-1,0,-1):        for i in range(N):            for j in range(M):                sum_beta[0,j] = A[i,j] * B[j,O[0,h]-1] * beta[j,h]            beta[i,h-1] = sum_beta.sum(1)    #print("beta矩陣: \n %r" % beta)    for i in range(N):        p_beta[0,i] = Pi[0,i] * B[i,O[0,0]-1] * beta[i,0]    p = p_beta.sum(1).reshape(1,1)    #print("觀測機率: \n %r" % p[0,0])    return beta, p[0,0]  def FBAlgoAppli(A,B,Pi,O,I):    #計算在觀測序列和模型參數確定的情況下，某一個隱含狀態對應相應的觀測狀態的機率    #例題參考李航《統計學習方法》P189習題10.2    #輸入格式：    #I為二維數組，存放所求機率P(it = qi,O|lambda)中it和qi的角標t和i，即P=[t,i]    alpha,p1 = ForwardAlgo(A,B,Pi,O)    beta,p2 = BackwardAlgo(A,B,Pi,O)    p = alpha[I[0,1]-1,I[0,0]-1] * beta[I[0,1]-1,I[0,0]-1] / p1    return p    def GetGamma(A,B,Pi,O):    N = A.shape[0]#數組A的行數    H = O.shape[1]#數組O的列數    Gamma = np.zeros((N,H))    alpha,p1 = ForwardAlgo(A,B,Pi,O)    beta,p2 = BackwardAlgo(A,B,Pi,O)    for h in range(H):        for i in range(N):            Gamma[i,h] = alpha[i,h] * beta[i,h] / p1    return Gamma    def GetXi(A,B,Pi,O):    N = A.shape[0]#數組A的行數    M = A.shape[1]#數組A的列數    H = O.shape[1]#數組O的列數    Xi = np.zeros((H-1,N,M))    alpha,p1 = ForwardAlgo(A,B,Pi,O)    beta,p2 = BackwardAlgo(A,B,Pi,O)    for h in range(H-1):        for i in range(N):            for j in range(M):                Xi[h,i,j] = alpha[i,h] * A[i,j] * B[j,O[0,h+1]-1] * beta[j,h+1] / p1    #print("Xi矩陣: \n %r" % Xi)    return Xi    def BaumWelchAlgo(A,B,Pi,O):    N = A.shape[0]#數組A的行數    M = A.shape[1]#數組A的列數    Y = B.shape[1]#數組B的列數    H = O.shape[1]#數組O的列數    c = 0    Gamma = GetGamma(A,B,Pi,O)    Xi = GetXi(A,B,Pi,O)    Xi_1 = Xi.sum(0)    a = np.zeros((N,M))    b = np.zeros((M,Y))    pi = np.zeros((1,N))    a_1 = np.subtract(Gamma.sum(1),Gamma[:,H-1]).reshape(1,N)    for i in range(N):        for j in range(M):            a[i,j] = Xi_1[i,j] / a_1[0,i]    #print(a)    for y in range(Y):        for j in range(M):            for h in range(H):                if O[0,h]-1 == y:                    c = c + Gamma[j,h]            gamma = Gamma.sum(1).reshape(1,N)            b[j,y] = c / gamma[0,j]            c = 0    #print(b)    for i in range(N):        pi[0,i] = Gamma[i,0]    #print(pi)    return a,b,pi    def BaumWelchAlgo_n(A,B,Pi,O,n):#計算迭代次數為n的BaumWelch演算法    for i in range(n):        A,B,Pi = BaumWelchAlgo(A,B,Pi,O)    return A,B,Pi    def viterbi(A,B,Pi,O):    N = A.shape[0]#數組A的行數    M = A.shape[1]#數組A的列數    H = O.shape[1]#數組O的列數    Delta = np.zeros((M,H))    Psi = np.zeros((M,H))    Delta_1 = np.zeros((N,1))    I = np.zeros((1,H))        for i in range(N):        Delta[i,0] = Pi[0,i] * B[i,O[0,0]-1]            for h in range(1,H):        for j in range(M):            for i in range(N):                Delta_1[i,0] = Delta[i,h-1] * A[i,j] * B[j,O[0,h]-1]            Delta[j,h] = np.amax(Delta_1)            Psi[j,h] = np.argmax(Delta_1)+1    print("Delta矩陣: \n %r" % Delta)    print("Psi矩陣: \n %r" % Psi)    P_best = np.amax(Delta[:,H-1])    psi = np.argmax(Delta[:,H-1])    I[0,H-1] = psi + 1    for h in range(H-1,0,-1):        I[0,h-1] = Psi[I[0,h]-1,h]    print("最優路徑機率: \n %r" % P_best)    print("最優路徑: \n %r" % I)

 

Python實現HMM（隱馬爾可夫模型）