隱馬爾科夫模型（HMM）一前向與後向演算法

01-24

在該篇文章中講了隱馬爾科夫模型（HMM）一基本模型與三個基本問題 - 知乎專欄隱馬爾科夫鏈的基本概念和最經典的三個問題，這篇文章總結一下隱馬爾科夫鏈（HMM）中的前向與後向演算法，首先給出這倆個演算法是為了解決HMM的第一個基本問題。先回憶一下第一個問題：

第一個問題是求，給定模型的情況下，求某種觀測序列出現的概率。
比如，給定的HMM模型參數已知，求出三天觀察是(Dizzy,Cold,Normal)的概率是多少？對應的HMM模型參數已知的意思，就是說的A(trainsition_probability),B(emission_probability),pi矩陣是已經知道的。

相關條件如下圖所示：

由上圖所示，也就是說，可以寫成如下代碼：

trainsition_probability = [[0.7,0.3],[0.4,0.6]] emission_probability = [[0.5,0.4,0.1],[0.1,0.3,0.6]] pi = [0.6,0.4]

在第一個問題中，我們需要求解出三天觀察是(Dizzy,Cold,Normal)的概率是多少？

這裡為了演示簡單，我只求解出倆天觀察為(Dizzy,Cold)的概率是多少！

這個問題太好求解了，最暴力的方法就是將路徑全部遍歷一遍。下面儘可能通俗易懂的說明一下：

首先畫出時間序列狀態圖如下：

下面，我詳細走一遍一條路徑的暴力演算法，這樣既可以避開公式的晦澀，也不失正確性。其它路徑完全類似

第一天為Healthy的概率為：0.6

在第一天為Healthy的基礎上，觀察為Dizzy的概率為：P（Dizzy|Healthy）=0.6P(Healthy->Dizzy)=0.60.1=0.06

然後求出在第一天為Healthy的基礎上，並且第一天表現為Dizzy的前提下，第二天也為Healthy的概率為：

P(Healthy|Healthy,Dizzy)=0.6p(Healthy->Healthy) = P(Dizzy|healthy)07 = 0.060.7

上面求完的時候，代表下圖中的紅線已經轉移完了。

好，那麼當在前面基礎上，第二天觀察為Cold的概率為：
P(Cold|(Healthy,Dizzy),(Healthy)) = P(Healthy|Healthy,Dizzy)0.4 = 0.060.70.4

現在我們已經完成一條路徑的完整結果了。

就是在第一天隱含狀態為Healthy和第二天隱含狀態為Healthy的基礎上，觀察序列為Dizzy，Cold的概率為

P（Dizzy,Cold|Healthy,Healthy） = 0.060.70.4=0.0168

那麼同理，我們就可以求出其它三條路徑。

（1）在第一天隱含狀態為Healthy和第二天隱含狀態為Fever的基礎上，觀察序列為Dizzy，Cold的概率

（2）在第一天隱含狀態為Fever和第二天隱含狀態為Healthy的基礎上，觀察序列為Dizzy，Cold的概率

（3）在第一天隱含狀態為Fever和第二天隱含狀態為Fever的基礎上，觀察序列為Dizzy，Cold的概率

然後最後的第一個問題的結果就是將這四個結果相加起來就可以了。是不是很簡單，那麼為了還需要前向後向演算法來解決這個事呢？

其實這個演算法在現實中是不可行的。我給的例子由於是為了講解容易，狀態值和觀察值都很小，但是實際中的問題，隱狀態的個數是非常大的。

那麼我們的計算量是不可以忍受的。

我們可以稍微估計一下，加入狀態值是N個，觀察值是K個。總共觀察長度為T。

那麼我們的路徑總個數就是 $N^{T}$ ,我的天，這個複雜度已經接受不了了，到達了每個隱含狀態的時候，還需要算一遍觀察值出現的概率（每個隱狀態計算一遍到觀察值的概率）。又要乘以NT（當然這已經對前面很大複雜度構成不了多少作用了）

所以我們得出結論，暴力法在這裡並不實用，於是就引出了前向後向演算法。它們都是基於動態規劃思想求解。下面分別介紹一下：

前向演算法

我們首先定義一下前向概率

定義：給定隱馬科夫模型 $lambda$ ，定義到時刻t為止的觀測序列為 $o_{1}$ , $o_{2},....,o_{t}$ 且狀態為 $q_{i}$ 的概率為前向概率，記作
可以遞推地求得前向概率 $alpha _{t} (i)$ 及觀測序列概率 $p(o|alpha )$ 。

下面，我們可以整理一下前向演算法的流程：

輸入：隱馬爾可夫模型 $lambda$ ，觀測序列 $O$

輸出：觀測序列概率 $p(o|alpha )$

(1)初值

前向概率的定義中一共限定了兩個條件，一是到當前為止的觀測序列，另一個是當前的狀態。所以初值的計算也有兩項（觀測和狀態），一項是初始狀態概率，另一項是發射到當前觀測的概率。

(2)遞推對t=1,2,3,.....,T-1

每次遞推同樣由兩部分構成，大括弧中是當前狀態為i且觀測序列的前t個符合要求的概率，括弧外的是狀態i發射觀測t+1的概率。

下面稍微解釋一下公式：

(3)終止

由於到了時間T，一共有N種狀態發射了最後那個觀測，所以最終的結果要將這些概率加起來（因為每個隱狀態都可能產生我們需要的觀測值，所以都要加起來）。

公式可以用下面的轉移圖表示，假設我要求第二層某個結點的前向概率，等於前一層所有結點到該結點的轉移，如下圖：

由於每次遞推都是在前一次的基礎上進行的，所以降低了複雜度（計算只存在於相鄰的倆個時間點）。計算如下圖所示：

下方標號表示時間節點，每個時間點都有N種狀態，所以相鄰兩個時間之間的遞推消耗N^2次計算。

而每次遞推都是在前一次的基礎上做的，所以只需累加O(T)次，所以總體複雜度是O(T)個N^2，即0（TN^2）,這比起我們前面說的暴力法的複雜度已經好了太多了。

後向概率

後向概率與前向概率非常類似，也是基於動態規劃的思想，下面介紹一下：

首先給出定義：

定義(後向概率）給定隱馬爾可夫模型 $lambda$ ,定義在時刻t狀態為 $q_{i}$ 的條件下，從t+1到T的部分觀測序列為 $o_{t+1},o_{t+2},...,o_{T}$ 的概率為後向概率，記作

可以用遞推的方法求得後向概率, $eta _{t} (i)$ 及觀測序列概率 $p(o|lambda )$ ,下面給出後向演算法的演算法流程。

輸入：隱馬爾可夫模型 $lambda$ ,觀測序列 $o$

輸出：觀測序列概率 $p(o|lambda )$

(1)初值

根據定義，從T+1到T的部分觀測序列其實不存在，所以硬性規定這個值是1。（這個很重要）

(2)對t = T-1,T-2,...1

$a_{ij}$ 表示狀態i轉移到j的概率， $b_{j}$ 表示發射 $o_{t+1}$ ， $eta _{t+1} (j)$ 表示j後面的序列對應的後向概率。

上面公式的求解可以用下圖轉移表示：假設我現在求第一層的某個結點的後向概率，等於第二層的所有結點轉移到該層結點，如下：

（3）

最後的求和是因為，在第一個時間點上有N種後向概率都能輸出從2到T的觀測序列，所以乘上輸出O1的概率後求和得到最終結果。（這些我在後面的代碼中均有對應的部分）

前向演算法和後向演算法的python編程實現如下：（代碼與上面的偽代碼有著非常清楚的對應關係）

# -*- coding: UTF-8 -*-#代碼是得到HMM的前向與後向演算法,已經驗證成功import numpy as npdef Forward(trainsition_probability,emission_probability,pi,obs_seq): """ :param trainsition_probability:trainsition_probability是狀態轉移矩陣 :param emission_probability: emission_probability是發射矩陣 :param pi: pi是初始狀態概率 :param obs_seq: obs_seq是觀察狀態序列 :return: 返回結果 """ trainsition_probability = np.array(trainsition_probability) emission_probability = np.array(emission_probability) print emission_probability[:,0] pi = np.array(pi) Row = np.array(trainsition_probability).shape[0] F = np.zeros((Row,Col)) #最後要返回的就是F，就是我們公式中的alpha F[:,0] = pi * np.transpose(emission_probability[:,obs_seq[0]]) #這是初始化求第一列,就是初始的概率*各自的發射概率 print F[:,0] for t in range(1,len(obs_seq)): #這裡相當於填矩陣的元素值 for n in range(Row): #n是代表隱藏狀態的 F[n,t] = np.dot(F[:,t-1],trainsition_probability[:,n])*emission_probability[n,obs_seq[t]] #對應於公式,前面是對應相乘 return Fdef Backward(trainsition_probability,emission_probability,pi,obs_seq): """ :param trainsition_probability:trainsition_probability是狀態轉移矩陣 :param emission_probability: emission_probability是發射矩陣 :param pi: pi是初始狀態概率 :param obs_seq: obs_seq是觀察狀態序列 :return: 返回結果 """ trainsition_probability = np.array(trainsition_probability) emission_probability = np.array(emission_probability) pi = np.array(pi) #要進行矩陣運算，先變為array類型 Row = trainsition_probability.shape[0] Col = len(obs_seq) F = np.zeros((Row,Col)) F[:,(Col-1):] = 1 #最後的每一個元素賦值為1 for t in reversed(range(Col-1)): for n in range(Row): F[n,t] = np.sum(F[:,t+1]*trainsition_probability[n,:]*emission_probability[:,obs_seq[t+1]]) return Fif __name__ == "__main__": trainsition_probability = [[0.7,0.3],[0.4,0.6]] emission_probability = [[0.5,0.4,0.1],[0.1,0.3,0.6]] pi = [0.6,0.4] #然後下面先得到前向演算法,在A,B,pi參數已知的前提下，求出特定觀察序列的概率是多少? obs_seq = [0,1] Row = np.array(trainsition_probability).shape[0] Col = len(obs_seq) F = Forward(trainsition_probability,emission_probability,pi,obs_seq) #得到前向演算法的結果 F_backward = Backward(trainsition_probability,emission_probability,pi,obs_seq) #得到後向演算法的結果 res_forward = 0 for i in range(Row): #將最後一列相加就得到了我們最終的結果 res_forward+=F[i][Col-1] #求和於最後一列 emission_probability = np.array(emission_probability) #下面是得到後向演算法的結果 res_backword = 0 res_backward = np.sum(pi*F_backward[:,0]*emission_probability[:,obs_seq[0]]) #一定要乘以發射的那一列 print "res_backward = {}".format(res_backward) print "res_forward = {}".format(res_forward)