鮑姆-韋爾奇算法求解HMM參數

1. HMM模型參數求解概述

　　　　HMM模型參數求解根據已知的條件可以分為兩種情況。

　　　　第一種情況較為簡單，就是我們已知DD個長度為TT的觀測序列和對應的隱藏狀態序列，即{(O1,I1),(O2,I2),...(OD,ID)}{(O1,I1),(O2,I2),...(OD,ID)}是已知的，此時我們可以很容易的用最大似然來求解模型參數。

　　　　假設樣本從隱藏狀態qiqi轉移到qjqj的頻率計數是AijAij,那么狀態轉移矩陣求得為：

A=[aij],其中aij=Aij∑s=1NAisA=[aij],其中aij=Aij∑s=1NAis

 

　　　　假設樣本隱藏狀態為qjqj且觀測狀態為vkvk的頻率計數是BjkBjk,那么觀測狀態概率矩陣為：

B=[bj(k)],其中bj(k)=Bjk∑s=1MBjsB=[bj(k)],其中bj(k)=Bjk∑s=1MBjs

 

　　　　假設所有樣本中初始隱藏狀態為qiqi的頻率計數為C(i)C(i),那么初始概率分布為：

Π=π(i)=C(i)∑s=1NC(s)Π=π(i)=C(i)∑s=1NC(s)

 

　　　　可見第一種情況下求解模型還是很簡單的。但是在很多時候，我們無法得到HMM樣本觀察序列對應的隱藏序列，只有DD個長度為TT的觀測序列，即{(O1),(O2),...(OD)}{(O1),(O2),...(OD)}是已知的，此時我們能不能求出合適的HMM模型參數呢？這就是我們的第二種情況，也是我們本文要討論的重點。它的解法最常用的是鮑姆-韋爾奇算法，其實就是基於EM算法的求解，只不過鮑姆-韋爾奇算法出現的時代，EM算法還沒有被抽象出來，所以我們本文還是說鮑姆-韋爾奇算法。

2. 鮑姆-韋爾奇算法原理

　　　　鮑姆-韋爾奇算法原理既然使用的就是EM算法的原理，那么我們需要在E步求出聯合分布P(O,I|λ)P(O,I|λ)基於條件概率P(I|O,λ¯¯¯)P(I|O,λ¯)的期望，其中λ¯¯¯λ¯為當前的模型參數，然后再M步最大化這個期望，得到更新的模型參數λλ。接着不停的進行EM迭代，直到模型參數的值收斂為止。

　　　　首先來看看E步，當前模型參數為λ¯¯¯λ¯, 聯合分布P(O,I|λ)P(O,I|λ)基於條件概率P(I|O,λ¯¯¯)P(I|O,λ¯)的期望表達式為：

L(λ,λ¯¯¯)=∑d=1D∑IP(I|O,λ¯¯¯)logP(O,I|λ)L(λ,λ¯)=∑d=1D∑IP(I|O,λ¯)logP(O,I|λ)

 

　　　　在M步，我們極大化上式，然后得到更新后的模型參數如下：　

λ¯¯¯=argmaxλ∑d=1D∑IP(I|O,λ¯¯¯)logP(O,I|λ)λ¯=argmaxλ∑d=1D∑IP(I|O,λ¯)logP(O,I|λ)

 

　　　　通過不斷的E步和M步的迭代，直到λ¯¯¯λ¯收斂。下面我們來看看鮑姆-韋爾奇算法的推導過程。

3. 鮑姆-韋爾奇算法的推導

　　　　我們的訓練數據為{(O1,I1),(O2,I2),...(OD,ID)}{(O1,I1),(O2,I2),...(OD,ID)}，其中任意一個觀測序列Od={o(d)1,o(d)2,...o(d)T}Od={o1(d),o2(d),...oT(d)},其對應的未知的隱藏狀態序列表示為：Od={i(d)1,i(d)2,...i(d)T}Od={i1(d),i2(d),...iT(d)}

　　　　首先看鮑姆-韋爾奇算法的E步，我們需要先計算聯合分布P(O,I|λ)P(O,I|λ)的表達式如下：

P(O,I|λ)=πi1bi1(o1)ai1i2bi2(o2)...aiT−1iTbiT(oT)P(O,I|λ)=πi1bi1(o1)ai1i2bi2(o2)...aiT−1iTbiT(oT)

 

　　　　我們的E步得到的期望表達式為：

L(λ,λ¯¯¯)=∑d=1D∑IP(I|O,λ¯¯¯)logP(O,I|λ)L(λ,λ¯)=∑d=1D∑IP(I|O,λ¯)logP(O,I|λ)

 

　　　　在M步我們要極大化上式。由於P(I|O,λ¯¯¯)=P(I,O|λ¯¯¯)/P(O|λ¯¯¯)P(I|O,λ¯)=P(I,O|λ¯)/P(O|λ¯),而P(O|λ¯¯¯)P(O|λ¯)是常數，因此我們要極大化的式子等價於：

λ¯¯¯=argmaxλ∑d=1D∑IP(O,I|λ¯¯¯)logP(O,I|λ)λ¯=argmaxλ∑d=1D∑IP(O,I|λ¯)logP(O,I|λ)

 

　　　　我們將上面P(O,I|λ)P(O,I|λ)的表達式帶入我們的極大化式子，得到的表達式如下：

λ¯¯¯=argmaxλ∑d=1D∑IP(O,I|λ¯¯¯)(logπi1+∑t=1T−1logaitait+1+∑t=1Tbit(ot))λ¯=argmaxλ∑d=1D∑IP(O,I|λ¯)(logπi1+∑t=1T−1logaitait+1+∑t=1Tbit(ot))

 

　　　　我們的隱藏模型參數λ=(A,B,Π)λ=(A,B,Π),因此下面我們只需要對上式分別對A,B,ΠA,B,Π求導即可得到我們更新的模型參數λ¯¯¯λ¯　

 

　　　　首先我們看看對模型參數ΠΠ的求導。由於ΠΠ只在上式中括號里的第一部分出現，因此我們對於ΠΠ的極大化式子為：

πi¯¯¯¯¯=argmaxπi1∑d=1D∑IP(O,I|λ¯¯¯)logπi1=argmaxπi∑d=1D∑i=1NP(O,i(d)1=i|λ¯¯¯)logπiπi¯=argmaxπi1∑d=1D∑IP(O,I|λ¯)logπi1=argmaxπi∑d=1D∑i=1NP(O,i1(d)=i|λ¯)logπi

 

　　　　由於πiπi還滿足∑i=1Nπi=1∑i=1Nπi=1，因此根據拉格朗日子乘法，我們得到πiπi要極大化的拉格朗日函數為：

argmaxπi∑d=1D∑i=1NP(O,i(d)1=i|λ¯¯¯)logπi+γ(∑i=1Nπi−1)argmaxπi∑d=1D∑i=1NP(O,i1(d)=i|λ¯)logπi+γ(∑i=1Nπi−1)

 

　　　　其中，γγ為拉格朗日系數。上式對πiπi求偏導數並令結果為0， 我們得到：

∑d=1DP(O,i(d)1=i|λ¯¯¯)+γπi=0∑d=1DP(O,i1(d)=i|λ¯)+γπi=0

 

　　　　令ii分別等於從1到NN，從上式可以得到NN個式子，對這NN個式子求和可得：

∑d=1DP(O|λ¯¯¯)+γ=0∑d=1DP(O|λ¯)+γ=0

 

　　　　從上兩式消去γγ,得到πiπi的表達式為：

πi=∑d=1DP(O,i(d)1=i|λ¯¯¯)∑d=1DP(O|λ¯¯¯)=∑d=1DP(O,i(d)1=i|λ¯¯¯)DP(O|λ¯¯¯)=∑d=1DP(i(d)1=i|O,λ¯¯¯)D=∑d=1DP(i(d)1=i|O(d),λ¯¯¯)Dπi=∑d=1DP(O,i1(d)=i|λ¯)∑d=1DP(O|λ¯)=∑d=1DP(O,i1(d)=i|λ¯)DP(O|λ¯)=∑d=1DP(i1(d)=i|O,λ¯)D=∑d=1DP(i1(d)=i|O(d),λ¯)D

 

　　　　利用我們在隱馬爾科夫模型HMM（二）前向后向算法評估觀察序列概率里第二節中前向概率的定義可得：

P(i(d)1=i|O(d),λ¯¯¯)=γ(d)1(i)P(i1(d)=i|O(d),λ¯)=γ1(d)(i)

 

　　　　因此最終我們在M步πiπi的迭代公式為：

πi=∑d=1Dγ(d)1(i)Dπi=∑d=1Dγ1(d)(i)D

 

 

　　　　現在我們來看看AA的迭代公式求法。方法和ΠΠ的類似。由於AA只在最大化函數式中括號里的第二部分出現，而這部分式子可以整理為：

∑d=1D∑I∑t=1T−1P(O,I|λ¯¯¯)logaitait+1=∑d=1D∑i=1N∑j=1N∑t=1T−1P(O,i(d)t=i,i(d)t+1=j|λ¯¯¯)logaij∑d=1D∑I∑t=1T−1P(O,I|λ¯)logaitait+1=∑d=1D∑i=1N∑j=1N∑t=1T−1P(O,it(d)=i,it+1(d)=j|λ¯)logaij

 

　　　　由於aijaij還滿足∑j=1Naij=1∑j=1Naij=1。和求解πiπi類似，我們可以用拉格朗日子乘法並對aijaij求導，並令結果為0，可以得到aijaij的迭代表達式為：

aij=∑d=1D∑t=1T−1P(O(d),i(d)t=i,i(d)t+1=j|λ¯¯¯)∑d=1D∑t=1T−1P(O(d),i(d)t=i|λ¯¯¯)aij=∑d=1D∑t=1T−1P(O(d),it(d)=i,it+1(d)=j|λ¯)∑d=1D∑t=1T−1P(O(d),it(d)=i|λ¯)

 

　　　　利用隱馬爾科夫模型HMM（二）前向后向算法評估觀察序列概率里第二節中前向概率的定義和第五節ξt(i,j)ξt(i,j)的定義可得們在M步aijaij的迭代公式為：

aij=∑d=1D∑t=1T−1ξ(d)t(i,j)∑d=1D∑t=1T−1γ(d)t(i)aij=∑d=1D∑t=1T−1ξt(d)(i,j)∑d=1D∑t=1T−1γt(d)(i)

 

 

　　　　現在我們來看看BB的迭代公式求法。方法和ΠΠ的類似。由於BB只在最大化函數式中括號里的第三部分出現，而這部分式子可以整理為：

∑d=1D∑I∑t=1TP(O,I|λ¯¯¯)logbit(ot)=∑d=1D∑j=1N∑t=1TP(O,i(d)t=j|λ¯¯¯)logbj(ot)∑d=1D∑I∑t=1TP(O,I|λ¯)logbit(ot)=∑d=1D∑j=1N∑t=1TP(O,it(d)=j|λ¯)logbj(ot)

 

　　　　由於bj(ot)bj(ot)還滿足∑k=1Mbj(ot=vk)=1∑k=1Mbj(ot=vk)=1。和求解πiπi類似，我們可以用拉格朗日子乘法並對bj(k)bj(k)求導，並令結果為0，得到bj(k)bj(k)的迭代表達式為：

bj(k)=∑d=1D∑t=1TP(O,i(d)t=j|λ¯¯¯)I(o(d)t=vk)∑d=1D∑t=1TP(O,i(d)t=j|λ¯¯¯)bj(k)=∑d=1D∑t=1TP(O,it(d)=j|λ¯)I(ot(d)=vk)∑d=1D∑t=1TP(O,it(d)=j|λ¯)

 

　　　　其中I(o(d)t=vk)I(ot(d)=vk)當且僅當o(d)t=vkot(d)=vk時為1，否則為0. 利用隱馬爾科夫模型HMM（二）前向后向算法評估觀察序列概率里第二節中前向概率的定義可得bj(ot)bj(ot)的最終表達式為：

bj(k)=∑d=1D∑t=1,o(d)t=vkTγ(d)t(i)∑d=1D∑t=1Tγ(d)t(i)bj(k)=∑d=1D∑t=1,ot(d)=vkTγt(d)(i)∑d=1D∑t=1Tγt(d)(i)

 

　　　　有了πi,aij,bj(k)πi,aij,bj(k)的迭代公式，我們就可以迭代求解HMM模型參數了。

4. 鮑姆-韋爾奇算法流程總結

　　　　這里我們概括總結下鮑姆-韋爾奇算法的流程。

　　　　輸入： DD個觀測序列樣本{(O1),(O2),...(OD)}{(O1),(O2),...(OD)}

　　　　輸出：HMM模型參數

　　　　1)隨機初始化所有的πi,aij,bj(k)πi,aij,bj(k)

　　　　2) 對於每個樣本d=1,2,...Dd=1,2,...D，用前向后向算法計算γ(d)t(i)，ξ(d)t(i,j),t=1,2...Tγt(d)(i)，ξt(d)(i,j),t=1,2...T

　　　　3)  更新模型參數：

 

πi=∑d=1Dγ(d)1(i)Dπi=∑d=1Dγ1(d)(i)D

 

 

aij=∑d=1D∑t=1T−1ξ(d)t(i,j)∑d=1D∑t=1T−1γ(d)t(i)aij=∑d=1D∑t=1T−1ξt(d)(i,j)∑d=1D∑t=1T−1γt(d)(i)

 

 

bj(k)=∑d=1D∑t=1,o(d)t=vkTγ(d)t(i)∑d=1D∑t=1Tγ(d)t(i)bj(k)=∑d=1D∑t=1,ot(d)=vkTγt(d)(i)∑d=1D∑t=1Tγt(d)(i)

 

　　　　4) 如果πi,aij,bj(k)πi,aij,bj(k)的值已經收斂，則算法結束，否則回到第2）步繼續迭代。

　　　　以上就是鮑姆-韋爾奇算法的整個過程。

 轉載：http://www.cnblogs.com/pinard/p/6972299.html

http://www.itdadao.com/articles/c15a132036p0.html