陳 琦, 胡錫健

(新疆大學(xué) 數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院, 新疆 烏魯木齊 830046)

隱馬爾可夫模型(Hidden Markov Model，HMM)是由Baum L E在20世紀60年代末提出的一個非常重要的統(tǒng)計分析模型，其在語音識別、故障診斷、生物統(tǒng)計等領(lǐng)域中有廣泛的應(yīng)用．評估問題、隱狀態(tài)估計(解碼)問題、參數(shù)估計(學(xué)習(xí))問題是HMM的三個基本問題.其中，隱狀態(tài)估計問題在某種程度上其實就是一個貝葉斯估計問題.當(dāng)HMM的隱狀態(tài)取連續(xù)值時，通常將其稱為狀態(tài)空間模型.對于線性高斯系統(tǒng)，可以利用卡爾曼濾波(Kalman Filter)得到解析解.對于狀態(tài)空間是離散的、有限的馬爾可夫鏈的系統(tǒng)，運用基于網(wǎng)格的濾波方法(Grid-based methods)可以得到后驗概率分布的解析解.但實際上這些條件在很多動態(tài)隨機系統(tǒng)中很難滿足，大部分是非線性、非高斯的，無法求出解析解，因此，只能求助于次優(yōu)或逼近方法.文獻[1]提出擴展卡爾曼濾波(Extended Kalman Filter,EKF),但EKF存在線性化誤差，只適用于弱非線性系統(tǒng).文獻[2]提出UT的變換，并在EKF濾波框架里得到廣泛應(yīng)用，被稱為UKF，這種方法不必對非線性方程進行線性化和計算雅克比矩陣，均值和方差均可精確到二階，比標(biāo)準(zhǔn)的EKF精度更高.但由于這兩種非線性估計方法都是把pk(xk|y1:k)用高斯逼近，如果真實的分布是非高斯的，那么就不能精確表示狀態(tài)的真實分布.然而，對于狀態(tài)空間是連續(xù)的但是可以劃分成有限的Ns個單元，那么仍然可以用基于網(wǎng)格的逼近方法來近似后驗概率密度[3].近年來，將粒子濾波用于處理非線性、非高斯系統(tǒng)的估計，為貝葉斯估計提供了有效的新解法.但當(dāng)隱狀態(tài)空間維數(shù)很高時離子濾波方法有以下兩個不足之處：計算量特別大；不能精確近似目標(biāo)分布.此外，離子濾波還要求似然函數(shù)是逐點可知的，但在實際情況中通常不能滿足，而用近似貝葉斯計算方法(Approximate Bayesian Compute,ABC)可以避免這些問題.本文將HMM隱態(tài)估計問題看作一個貝葉斯濾波問題，首先詳細介紹粒子濾波用于HMM隱狀態(tài)估計，然后將近似貝葉斯計算思想引入HMM隱狀態(tài)估計中.

1 問題描述

HMM是一對離散時間隨機過程{Xn}n≥1、{Yn}n≥1，其中xn∈dx是隱狀態(tài)，yn∈dy是觀測值。HMM隱狀態(tài)估計，就是給定觀測序列y0:n及模型參數(shù)θ，對隱狀態(tài)作出恰當(dāng)?shù)呐袛啵丛谀撤N意義下求與此觀測序列對應(yīng)的最優(yōu)隱狀態(tài)序列x0:n.這一問題的回答取決于對模型狀態(tài)物理意義的理解，根據(jù)不同的準(zhǔn)則，可能會選出不同的最優(yōu)狀態(tài)序列．

對?0≤k≤l及n>0,定義φk:l|n(y0:n，·)為給定Y0:n的條件下Xk:l的條件分布，對于具有轉(zhuǎn)移密度函數(shù)的HMM來說，φk:l|n(y0:n，·)是Yn+1到Xl-k+1的轉(zhuǎn)移核，且有密度函數(shù)φk:l|n(xk:l|y0:n).對于隱狀態(tài)估計問題有兩種原則，單點最優(yōu)原則對應(yīng)于條件分布φk|n(y0:n，·)=φk:k|n(y0:n，·),0≤k≤n；路徑最優(yōu)原則對應(yīng)于條件分布φ0:n|n(y0:n，·).

單點最優(yōu)原則下的隱狀態(tài)估計問題其實是一個最優(yōu)濾波問題[4]，HMM模型的隱狀態(tài)是由φk|n(xk|y0:n)而得到的，即時刻k(0≤k≤n)下最有可能的隱狀態(tài)為

(1)

(2)

2 基于離子濾波的HMM隱狀態(tài)估計

粒子濾波算法是通過非參數(shù)化的蒙特卡洛(MonteCarlo)模擬方法來實現(xiàn)遞推的貝葉斯濾波，適用于任何能用狀態(tài)空間模型描述的非線性系統(tǒng)，其核心思想是得到一批有權(quán)重的離散隨機采樣點(粒子)來近似狀態(tài)變量的后驗概率密度函數(shù)，然后利用他們得到估計值，當(dāng)?shù)玫降牧Ｗ雍芏鄷r，其精度可以逼近最優(yōu)估計.

2.1 離子濾波估計過程描述

HMM隱狀態(tài)估計的核心是后驗分布φ0:n|n(y0:n，·)，假設(shè)其密度函數(shù)存在，則

(3)

(4)

L(y0:n)φ0:n|n(x0:n|y0:n)=L(y0:n-1)φ0:n-1|n-1(x0:n-1|y0:n-1)q(xn-1,xn)g(xn,yn)

(5)

(6)

得φk|n(xk|y0:n)的估計為

(7)

2.2 建議分布選取

故(4)式的最優(yōu)建議分布為

(8)

2.3 重抽樣

在進行序貫重要性采樣時，隨著 k的增大，將會出現(xiàn)權(quán)值退化問題，當(dāng)n較大時，就將大量的時間浪費在一些對估計結(jié)果不起任何作用的粒子更新上，針對權(quán)值退化問題，Gordon等提出采用重抽樣(Resampling)方法，即評估粒子權(quán)值后重抽樣粒子以減少小權(quán)值的粒子而復(fù)制大權(quán)值的粒子.

2.4 粒子濾波算法

(a)當(dāng)k=0時

(b)當(dāng)k=1,2,…，n時

第k步抽樣：

第k步重抽樣：

由于

(9)

由(7)式知

(10)

(11)

設(shè)

(12)

(13)

3 基于ABC(Approximate Bayesian Compute)的HMM隱狀態(tài)估計

3.1 ABC的引入

粒子濾波算法是基于重點抽樣來近似目標(biāo)分布的序列，通常目標(biāo)分布的維數(shù)隨著時間參數(shù)的增加而增加.在維數(shù)很高時，離子濾波有計算量大、不能準(zhǔn)確地近似目標(biāo)分布等缺點.第一個缺點歸因于進行模擬和估計權(quán)值時的計算量都很大，第二個缺點是由于目標(biāo)分布本身的復(fù)雜性.此外，離子濾波通常要求g(xk,yk)是逐點可得的，但在一些應(yīng)用中不能滿足該要求.采用基于近似貝葉斯計算的方法可以解決此類問題.

近似貝葉斯計算已經(jīng)成為貝葉斯范例的一個重要方面，在人們感興趣的許多實際問題中，貝葉斯模型不可得，或者由于計算量太大而不能被擬合，或者所需的模擬方法不能處理問題的復(fù)雜性.文獻[5]介紹了后驗密度π(x|y)∝g(y|x)μ(x)的如下近似：

(14)

其中：Aε,y={u:ρ(s(u),s(y))<ε}；u是一個輔助數(shù)據(jù)點，∏A是集合A的示性函數(shù)，g是似然函數(shù)，μ是先驗密度，s:dy→ds是一個數(shù)據(jù)的匯總統(tǒng)計量，ρ是一個距離測度.通常ds比dy小很多.如果s是充分的，可以看到當(dāng)ε趨于0時，πε(x|y)趨于π(x|y).可以采用拒絕抽樣[5]、馬爾可夫蒙特卡洛[6]或序貫蒙特卡洛[7-8])等方法對(14)式進行統(tǒng)計推斷，其中的關(guān)鍵點是，若不能對g(u|x)進行數(shù)值估計，但對任一x，如果可以從g(·|x)進行抽樣，仍就可以從πε(x|y)生成樣本.

3.2 ABC目標(biāo)的構(gòu)建

πn(x1:n|y1:n)的邊際ABC目標(biāo)為

(15)

3.3 ABC目標(biāo)的SMC算法

第1步：

4 結(jié)束語

文中將HMM的隱狀態(tài)估計問題看成一個貝葉斯估計問題，在單點原則下就是一個最優(yōu)濾波問題，而解決濾波問題的方法很多，但適用的條件及范圍各不相同，其得到的估計結(jié)果精度也不相同.離子濾波可以很好地對大部分HMM的隱狀態(tài)的后驗密度進行估計，但其也有自身的一些缺點，并且當(dāng)遇到g(xk,yk)逐點不可得時，無法對后驗密度進行估計.而基于ABC的方法可以對此類情況進行估計，從而使HMM隱狀態(tài)估計問題得到了全面解決，而且也可以考慮將ABC方法用到其它更廣泛的領(lǐng)域中.

[1] Anderson B D,Moore J B.Optimal filtering[M].New Jersey:Prentice-Hall,1979.

[2] Julier S J, Uhlmann J K.A new extension of the Kalman filter to nonlinear systems[C]// The 11th Int Symp on Aerospace:Simulation and Controls. London: 1997:10-12.

[3] 張詩桂，朱立新，趙義正.離子濾波算法研究進展與展望[J].自動化技術(shù)與應(yīng)用 ,2010,29(6):1-16.

[4] 朱成文，李兵，胡奎，等.HMM隱狀態(tài)的粒子濾波估計[J].計算機工程與應(yīng)用,2012,48(8):161-167.

[5] Pritchard J K, Seielstad M T, Perez-Lezaun A,etal. Population growth of human Y chromosome microsatellites[J]. Mol Biol Evol.,1999,16, 1 791-1 798.

[6] Majoram P, Molitor J , Plagnol V,etal. Markov chain Monte Carlo without likelihoods[J]. Proc Natl Acad Sci,2003,100, 15 324-15 328.

[7] Beaumont M A, Cornuet J M, Marin J M, et al.Adaptivity for ABC algorithms: the ABC- PMC scheme[J]. Biometrika, 2009,96, 983-990.

[8] Doucet A, Johansen A M. A tutorial on particle filtering and smoothing: fifteen years later[J]. Oxford Handbook of Nonlinear Filtering,2009,7：54-89.

[9] Jasra A ,Singh S S,Martin J S,etal. Filtering via approximate Bayesian computation[J]. Statistics and Computing, 2012,22(6):1 223-1 237.

[10] Moral P D, Doucet A, Jasra A. An adaptive sequential Monte Carlo method for approximate Bayesian computation[J]. Statistics and Computing, 2012,22(5):1 009-1 020.