一站式解決:隱馬爾可夫模型(HMM)全過程推導及實現

點擊上方「MLNLP」，選擇「星標」公眾號

重磅乾貨，第一時間送達

【導讀】隱馬爾可夫模型（Hidden Markov Model，HMM）是關於時許的概率模型，是一個生成模型，描述由一個隱藏的馬爾科夫鏈隨機生成不可觀測的狀態序列，每個狀態生成一個觀測，而由此產生一個觀測序列。

定義抄完了，下面我們從一個簡單的生成過程入手，順便引出HMM的參數。

假設有4個盒子，每個盒子裡面有不同數量的紅、白兩種顏色的球，具體如下表：

本慄子引用自《統計學習方法》

現在從這些盒子中抽取若干（  ）個球，每次抽取後記錄顏色，再放回原盒子，採樣的規則如下：

開始時，按照一個初始概率分布隨機選擇第一個盒子，這裡將第一個盒子用  表示：

將的值用變量 表示。因為有4個盒子可共選擇，所以 。然後隨機從該盒子中抽取一個球，使用表示：

將 的值用變量 表示。因為只有兩種球可供選擇，所以 。一共有4個箱子，2種球，結合前面的箱子的詳細數據，可以得到從每一個箱子取到各種顏色球的可能性，用一個表格表示：

進一步，可以用一個矩陣（稱為觀測概率矩陣，也有資料叫做發射矩陣）來表示該表

其中  表示在當前時刻給定 的條件下，給定

 表示當前的時刻，例如現在是第1時刻；然後是前面標註的初始概率分布，這個概率分布可以用一個向量（稱作初始狀態概率向量）來表示：

其中的

例如該分布是均勻分布的話，對應的向量就是

記錄抽取的球的顏色後將其放回，然後在按照如下規則選擇下一個盒子（）：

如果當前是盒子2或3，則分布以概率0.4和0.6選擇前一個或後一個盒子如果當前是盒子4，則各以0.5的概率停留在盒子4或者選擇盒子3
同樣，也可以根據以上規則做出一個表格，其中首列表示當前盒子，首行表示下一個盒子

以上，生成過程的主要流程就介紹完了，簡單概括就是：盒子，取球，盒子，取球……直到生成指定數量（T）的數據後停止。如果對這個過程還有不太理解的話，可以看看文章開頭給出的關於馬爾科夫鏈的連結。

其中N表示隱變量z的狀態數量，M表示觀測變量x可能的取值數量，在後面的討論中，用表示所有的參數。下面，根據這個慄子，寫一個數據生成代碼

import numpy as np
class HMM(object):    def __init__(self, N, M, pi=None, A=None, B=None):        self.N = N        self.M = M        self.pi = pi        self.A = A        self.B = B
    def get_data_with_distribute(self, dist):         r = np.random.rand()        for i, p in enumerate(dist):            if r < p: return i            r -= p
    def generate(self, T: int):        '''        根據給定的參數生成觀測序列        T: 指定要生成數據的數量        '''        z = self.get_data_with_distribute(self.pi)            x = self.get_data_with_distribute(self.B[z])          result = [x]        for _ in range(T-1):                    z = self.get_data_with_distribute(self.A[z])            x = self.get_data_with_distribute(self.B[z])            result.append(x)        return result
if __name__ == "__main__":    pi = np.array([.25, .25, .25, .25])    A = np.array([        [0,  1,  0, 0],        [.4, 0, .6, 0],        [0, .4, 0, .6],        [0, 0, .5, .5]])    B = np.array([        [.5, .5],        [.3, .7],        [.6, .4],        [.8, .2]])    hmm = HMM(4, 2, pi, A, B)    print(hmm.generate(10))   # 生成結果如下[0, 0, 1, 1, 1, 1, 0, 1, 0, 0]   
現在，參數介紹完了，數據生成過程也了解了，接下來就是解決HMM的基本問題了，一共有三個
不過，在討論這三個問題的相關算法之前，首先要給出兩個假設，在後面的推導過程中會不斷的用到：
齊次馬爾可夫假設：即任意時刻的狀態只依賴於前一時刻的狀態，與其他時刻的狀態無關（當然，初始時刻的狀態由參數π決定）：
觀測獨立假設：即任意時刻的觀測只依賴於該時刻的狀態，與其他無關：概率計算算法

概率計算算法現在，來看第一個問題，關於概率的計算，由於存在隱變量，所以X的邊際概率需要將所有的聯合概率  加和得到：
由於給出了T個觀測數據，所以相應的狀態也有T 個：
將（1）式中的 展開得到：
即使不考慮內部的計算，這起碼也是 階的計算量，所以需要更有效的算法，下面介紹兩種：前向算法和後向算法。
現在定義一個前向概率 ，它t時刻的狀態以及1,2,...t時刻的觀測在給定參數下的聯合概率：
也就是下圖中標記的那一部分
其中  表示由狀態 生成給定觀測數據的概率，例如設t時刻觀測數據  ，有
接著，根據（2）式，還可以得到：
由此公式，遍歷的取值求和，可以得到  的邊際概率
首先，來看上式中紅色部分，根據觀測獨立假設（下文再引用該假設時稱作假設2）：
然後是藍色部分，根據齊次馬爾可夫假設（下文再引用該假設時稱作假設1）
以上就是前向算法的推導，下面根據一個慄子來寫代碼，假設前面抽了五個球，分別是：紅、紅、白、白、紅，求概率
class HMM(object):    def evaluate(self, X):        '''        根據給定的參數計算條件概率        X: 觀測數據        '''        alpha = self.pi * self.B[:,X[0]]        for x in X[1:]:                                                            alpha = np.sum(self.A * alpha.reshape(-1,1) * self.B[:,x].reshape(1,-1), axis=0)        return alpha.sum()
print(hmm.evaluate([0,0,1,1,0]))   
上面的注釋中的代碼是按照公式來寫的，可以看出，時間複雜度降為了 ，比之前至少 的起步價已經好太多了.


另外說一點，如果對於前向算法還有印象的話，你會發現：上面的定義。實際上，對於任意時刻t，存在以下等式
觀察上式，藍色部分自然就是 。而紅色部分，根據假設2，都是無關（即互相獨立），可以省去，所以這部分最終變為：
    
def evaluate_backward(self, X):        beta = np.ones(self.N)        for x in X[:0:-1]:            beta_next = np.empty(self.N)            for i in range(self.N):                beta_next[i] = np.sum(self.A[i,:] * self.B[:,x] * beta)            beta = beta_next        return np.sum(beta * self.pi * self.B[:,X[0]])
和前向算法差不多，而且是照著公式寫的，就不寫注釋了，還是使用前面的慄子，跑了一下發現結果是一樣的。我想，同時寫兩個BUG得出同一個結果的概率應該很小很小吧
學習算法現在，概率計算的問題就解決了，接著來看第二個問題，參數學習，這裡需要用到EM算法，不熟悉的可以參考一下：https://zhuanlan.zhihu.com/p/85236423
然後，對Q函數中的每一項進行化簡，首先是第一項，用到了齊次馬爾可夫假設：
其中， 就是當前參數下觀測數據的概率，就是第一個問題所求解的。另外，利用第一個問題中定義的前向概率和後向概率，有：
最終得到：
同樣，將上式化簡，另外為了在後面方便引用，將該式設為一個函數f
可以看到，一共是  個相似的項，我們提一個（紅色部分）出來化簡，看看能不能找到通項公式
因為  是一個概率分布的矩陣，例如前面的慄子，每一行的和等於1將該函數對矩陣A的每一個元素求（偏）導並令導數為0：
將兩邊同時乘上  ，得到
注意一下上面的下標t與上標中（t+1）它們是不同的，由於變量比較多，各種ijk比較多，所以這裡需要注意一下。然後利用的約束，代入（6）式，得到：
然後化簡：
代入（7）式，得到
（8）式中，分母部分前面已經解決了，下面來看分子部分，進行化簡
注意，上面的化簡中，  。然後紅色和藍色部分的化簡用到了前面前面提過的兩個假設，將條件中不被依賴的變量去掉了。最後代入（8）式得到：
M是矩陣B的列數，前面已經定義過的，構造拉格朗日函數：
化簡得到（9）式的分母
class HMM(object):    def fit(self, X):        '''        根據給定觀測序列反推參數        '''                self.pi = np.random.sample(self.N)        self.A = np.ones((self.N,self.N)) / self.N        self.B = np.ones((self.N,self.M)) / self.M        self.pi = self.pi / self.pi.sum()        T = len(X)        for _ in range(50):                        alpha, beta = self.get_something(X)            gamma = alpha * beta
            for i in range(self.N):                for j in range(self.N):                    self.A[i,j] = np.sum(alpha[:-1,i]*beta[1:,j]*self.A[i,j]*self.B[j,X[1:]]) / gamma[:-1,i].sum()
            for j in range(self.N):                for k in range(self.M):                    self.B[j,k] = np.sum(gamma[:,j]*(X == k)) / gamma[:,j].sum()                        self.pi = gamma[0] / gamma[-1].sum()

    def get_something(self, X):        '''        根據給定數據與參數，計算所有時刻的前向概率和後向概率        '''        T = len(X)        alpha = np.zeros((T,self.N))        alpha[0,:] = self.pi * self.B[:,X[0]]        for i in range(T-1):            x = X[i+1]            alpha[i+1,:] = np.sum(self.A * alpha[i].reshape(-1,1) * self.B[:,x].reshape(1,-1), axis=0)
        beta = np.ones((T,self.N))        for j in range(T-1,0,-1):            for i in range(self.N):                beta[j-1,i] = np.sum(self.A[i,:] * self.B[:,X[j]] * beta[j])
        return alpha, beta
if __name__ == "__main__":    import matplotlib.pyplot as plt    def triangle_data(T):           data = []        for x in range(T):            x = x % 6            data.append(x if x <= 3 else 6-x)        return data        data = np.array(triangle_data(30))    hmm = HMM(10, 4)    hmm.fit(data)                   gen_obs = hmm.generate(30)      x = np.arange(30)    plt.scatter(x, gen_obs, marker='*', color='r')    plt.plot(x, data, color='g')    plt.show()
上面的代碼，使用最開始的慄子無法收斂，或者收斂到坑裡（公式和書上《統計學習方法》是一樣的），但是使用別人的例子又能很好的工作。調了一晚上後我覺得還是把它貼上來算了，希望大神發現了問題所在能告知一下。
預測算法最後一個問題了，解決這個問題的算法叫做維特比（Viterbi）算法。實際上它是一個動態規劃求解最優路徑的算法，這裡的最優路徑不過就是對應成最大概率而已，比前面兩個問題容易解決得多。直接上例子，如下圖所示：
如果  ，那麼最優路徑（索引）自然就是
接著，假設還有  ，末端的計算自然還是一樣問題在於從 如何計算最大概率
然後，又因為我們所求的是路徑，所以還要記錄最大概率所對應的索引值
class HMM(object):    def decode(self, X):        T = len(X)        x = X[0]        delta = self.pi * self.B[:,x]        varphi = np.zeros((T, self.N), dtype=int)        path = [0] * T        for i in range(1, T):            delta = delta.reshape(-1,1)                 tmp = delta * self.A            varphi[i,:] = np.argmax(tmp, axis=0)            delta = np.max(tmp, axis=0) * self.B[:,X[i]]        path[-1] = np.argmax(delta)                for i in range(T-1,0,-1):            path[i-1] = varphi[i,path[i]]        return path
推薦閱讀：
關於句子表徵的學習筆記
常用的 Normalization 方法：BN、LN、IN、GN
PaddlePaddle實戰NLP經典模型 BiGRU + CRF詳解