第19卷第9期計(jì)算機(jī)技術(shù)與發(fā)展Vol 19 No 9CoMPUTER TECHNOLOGY AND DEVELOPMENTSep.2009EM算法研究與應(yīng)用王愛(ài)平,張功營(yíng),劉方(安徽大學(xué)計(jì)算智能與信號(hào)處理教育部重點(diǎn)實(shí)驗(yàn)室,安徽合肥230039)摘要:引入了可處理缺失數(shù)據(jù)的EM算法。EM算法是一種迭代算法,每一次迭代都能保證似然函數(shù)值增加,并且收斂到一個(gè)局部極大值。對(duì)EM算法的基本原理和實(shí)施步驟進(jìn)行了分析。算法的命名,是因?yàn)樗惴ǖ拿恳坏▋刹降诓角笃谕? Expectation Step),稱為E步;第二步求極大值( Maximization Step),稱為M步。EM算法主要用來(lái)計(jì)算基于不完全數(shù)據(jù)的極大似然估計(jì)。在此基礎(chǔ)上把M算法融合到狀態(tài)空間模型的參數(shù)估計(jì)問(wèn)題。給出了基于 Kalman平滑和算法的線性狀態(tài)空間模型參數(shù)估計(jì)方法。關(guān)鍵詞EM算法;狀態(tài)空間模型; Kalman中圖分類號(hào):TP301.6獻(xiàn)標(biāo)識(shí)碼:A文章編號(hào):1673-629X(2009)09-0108-03Research and application of EM algorithmWANG Ai-ping, ZHANG Gong-ying, LIU Fang( Ministry of Education Key Lab. of Intelligent Computing Signal ProcessingAnhui University, Hefei 230039, China)Abstrac: Following the description of traditional maximum likelihood estimation methods and the discussions on their disadvantages. EMalgorithm is an iterative algorithm, every iteration to ensure that the likelihood function can be increased, and the oonvergenoe to a localaximA. Presents an EM algorithm that can be used to deal with missing data problems, where the details of the EM algorithm and its realization procedure have been analyzed. Algorithm named because each iterative algorithm includes two steps: the first step in seeking expectations( Expectation Step), known as the e step; the second step for maxima( Maximization Step), kIown as step-by-step M.EM algorithm used to calculate the principel based on incomplete data, maximum likelihood estimation. This is then followed by applyingthe proposed em algorithm to the parameter estimation of state pace models. The paper also presents the Kalman smoothing basedrameter estimation methods for linear state space modelsKey words: EM algorithm; state- space model; Kalman1EM算法及其應(yīng)用間模型的參數(shù)估計(jì)問(wèn)題EM算法是一種迭代算法,每一次迭代都能保證(1)似然函數(shù)值增加,并且收斂到一個(gè)局部極大值1,2。y, Hrr, Ur(2)算法的命名是因?yàn)樗惴ǖ拿恳坏▋蓚€(gè)步:第一其中初始狀態(tài)x0,隨機(jī)誤差項(xiàng)a和v為獨(dú)立高斯分步求期望( Expectation Step),稱為E步;第二步求極大布,三者之間是相互獨(dú)立:xo~N(x00,P0),at計(jì)算基于不完全數(shù)據(jù)的極大似然估計(jì)要用來(lái)N0,Q),-iN0,R),E【u]=0,E[uxo]值( Maximization Step),稱為M步。EM算法=0,E[ax0]=0。假設(shè)模型噪聲2狀態(tài)的初始條件都服從高斯2EM算法的狀態(tài)空間模型參數(shù)估計(jì)分布,即:下面推導(dǎo)基于 Kalman平滑和算法的狀態(tài)空間模p(ro| 8)=(2x)2 I Poro I-to(x0型參數(shù)估計(jì)方法6-10。給出一般形式的線性狀態(tài)空)P((3)收稿日期:2008-12-26;修回日期:200-03-24基金項(xiàng)目:國(guó)家自然科學(xué)基金項(xiàng)目(60472065)dH中國(guó)煤化工∞p}-號(hào)(CNMHG(4)作者簡(jiǎn)介:王愛(ài)平(1956-),女安徽合肥人,教授,主要從事計(jì)算機(jī)教學(xué)與研究。p(x!1x1,0)=(2)號(hào)1Qs-(x第9期王愛(ài)平等:EM算法研究與應(yīng)用109Fx-1)Q1(x-Fx1-1)(5) HrT)+ HPurh,])則完全數(shù)據(jù)的似然也服從高斯分布,完全數(shù)據(jù)的(10)似然可以表示為:p(x0r,y16)=p(x010)Ⅱp(x1x1,4算法實(shí)施步驟8) p(y, I xp, o)(6)現(xiàn)在把應(yīng)用EM算法進(jìn)行狀態(tài)空間模型參數(shù)因此,完全數(shù)據(jù)的對(duì)數(shù)似然等于:估計(jì)的實(shí)施步驟歸納如下:LnP(Io:T, y1:T I8)=1)對(duì)參數(shù)0={R,Q,F,H,x00,Pool進(jìn)行初始化In I Pono 1-4(xo-IoIo)Polo(xo-xouo)2)E步:根據(jù):-1,執(zhí)行卡爾曼濾波及平滑,計(jì)算吾Q=2x-Em)xB斯值3)M步:計(jì)算參數(shù)={R,Q,F,H,x0,P0新TInIr-1 2( -,x, yR-1( -, 3,的值;n+mN+ min(2r)4)重復(fù)參數(shù)直至收斂14。其中,N是樣本量,m是狀態(tài)變量的維度,n是觀測(cè)變5實(shí)驗(yàn)仿真量的維度?？紤]如下的普通線性高斯?fàn)顟B(tài)1空間模型取3求解對(duì)數(shù)似然的期望的極大化根據(jù)EM算法基本原理,獲得以下式子y= HI t vr而狀態(tài)初值x0、噪聲干擾v2和v均被假設(shè)為高ELInp(Io: T, 1T I 0)]=-5hn I Pc斯分布,即:InIQ1-hnIR|(n+ m)N+mln(2x)xo -N(oo, Poro), w iid N(o, Q),v"iidtr(PoEL(xo-Ioio)(xo-Ioio)'J)N(0,R)其中參數(shù)實(shí)際取值為:F=0.9,H=0.2,Q=0.7-∑(QE[(x1-Fx-1)(x-Fx-1))R=0.2,狀態(tài)初始設(shè)置為x0~N(0,10)。1S(RE[(y-Hx)(y-Hx)1)(8)依據(jù) Kalman濾波和平滑公式,然后聯(lián)合EM算法對(duì)模型參數(shù)日=|F,H,Rn,Qu}進(jìn)行估計(jì)EM算法迭代次數(shù)設(shè)置為400(當(dāng)然選代次數(shù)設(shè)置越大,能得到更A=2(ur a P)好的估計(jì)結(jié)果)。程序運(yùn)行結(jié)果θ=10.9835,0.2049,0.2107,0B=∑(x1xt1r+P-1r)683},仿真結(jié)果如圖1~3所示C=∑(x1rgx11r+P:-1r)這樣,完全數(shù)據(jù)的對(duì)數(shù)似然期望可以表示為:ELInt(T, yi:TIni QI-lI R I(2x)-ftr(POOEIGCoT -Io)(Ic中國(guó)煤化工Por)-2∑(QCNMHGFBFJ)tr([( -Hrr)(y圖1似然函教的期望的Q(日,-1)收斂性110計(jì)算機(jī)技術(shù)與發(fā)展第19卷用D]廣州暨南大學(xué),2007[2]張杰,陽(yáng)憲惠.多變量統(tǒng)計(jì)過(guò)程控制[M].北京:化工工業(yè)出版社,20[3]陳志國(guó),劉文舉基于EM算法的極大似然參數(shù)估計(jì)探討J河南大學(xué)學(xué)報(bào):自然科學(xué)版,2002,32(4):35-41[4]孫大飛, Dempster A P, Laird M,etl.MaxiSeries b,1997,39(1):1-38.[5] Meng X L, Rubin D B. Recent Extension to theEM algorithm[M]. Bayesian Statistics 4.Oord圖2參數(shù)F真實(shí)值(粗線)和估計(jì)值[6]胡壽松,王執(zhí)銓,胡維禮.最優(yōu)控制理論與系統(tǒng)[M].北京科學(xué)出版社,2005[7]郭雷,程代展,馮德興控制理論導(dǎo)論[M]北京:科學(xué)出版社,2005:1-1078 Andrieu C, Doucet A Online Expection-Maximization Type algorithms for Parameter Estimation inGeneral State Space Models[C]//in Proc. IEEE[s.L.]:[s.n.],2003:69-72.[9]賈沛璋朱征桃最優(yōu)估計(jì)及其應(yīng)用[M]北京科學(xué)出版社,1994[10] Parzen E On the estimation of a probability densityfunction andmode [J]. annals of Mathematical圖3參數(shù)H的收效過(guò)程Statistics,1962,33:1065-1076從上面的仿真圖形可以看出,EM算法用于狀態(tài)[1]waAP, Wang H. Minimising entropy and mean tracking空間模型參數(shù)估計(jì),其估計(jì)結(jié)果還是不錯(cuò)的,可以為時(shí)control for affine nonlinear and non Gaussian dynamic變不同的參數(shù)模型估計(jì)提供很好的幫助。stochastic system[J]. IEE Proceedings Control Theory & Ap-plication,2005,151(4):405-520[12] Wang A P, Wang H, Tan J. Optimal Filtering for Multivari6結(jié)束語(yǔ)ble Stochastic System via residual Probability Density Func-針對(duì)傳統(tǒng)的極大似然參數(shù)估計(jì)方法在解決實(shí)際問(wèn)ion Shaping[c]//Proceedings of SICE 2005 Annual Confer-題時(shí)所存在的不足,引入了EM算法。討論了算法原ence.[sL]:[sn.],2005:215-219理和實(shí)施步驟,并把EM算法推廣到狀態(tài)空間模型的13]任光張均東時(shí)間序列狀態(tài)空間模型建模及應(yīng)用(M參數(shù)估計(jì)問(wèn)題,然后推導(dǎo)了基于Kamn平滑和EM算大連:大連海事大學(xué)出版社,2000102-112法的狀態(tài)空間模型參數(shù)估計(jì)的過(guò)程,最后給出實(shí)驗(yàn)仿14cuL. Wang H. Mininum entropy filtering for multivariatestochastic systems with non-Gaussian noises [J].IEEETransactions on Automatic Control, 2006, 51(4): 670-695[15]梁應(yīng)敞張賢達(dá)李衍達(dá),等非高斯相關(guān)噪聲中高斯信號(hào)參考文獻(xiàn):的時(shí)延估計(jì)[門]電子科學(xué)學(xué)刊,1997,19(5):607-612[1]陳學(xué)華狀態(tài)空間模型理論與算法及其在金融計(jì)量中的應(yīng)多史國(guó)計(jì)算機(jī)學(xué)會(huì)余我史《計(jì)算機(jī)技術(shù)與發(fā)展歡迎訂閱VLl中國(guó)煤化工CNMHG