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

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的警告，還沒有去掉這個警告，查了一下說不影響結果，后面會去解決這個問題，下面貼出我的代碼

# -*- coding: utf-8 -*-
"""
Created on Thu Feb 16 19:28:39 2017
2017-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)