EM算法研究与应用

第19卷第9期计算机技术与发展Vol 19 No 9CoMPUTER TECHNOLOGY AND DEVELOPMENTSep.2009EM算法研究与应用王爱平,张功营,刘方(安徽大学计算智能与信号处理教育部重点实验室,安徽合肥230039)摘要:引入了可处理缺失数据的EM算法。EM算法是一种迭代算法,每一次迭代都能保证似然函数值增加,并且收敛到一个局部极大值。对EM算法的基本原理和实施步骤进行了分析。算法的命名,是因为算法的每一迭代包括两步第步求期望( Expectation Step),称为E步;第二步求极大值( Maximization Step),称为M步。EM算法主要用来计算基于不完全数据的极大似然估计。在此基础上把M算法融合到状态空间模型的参数估计问题。给出了基于 Kalman平滑和算法的线性状态空间模型参数估计方法。关键词EM算法;状态空间模型; Kalman中图分类号:TP301.6献标识码:A文章编号: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算法及其应用间模型的参数估计问题EM算法是一种迭代算法,每一次迭代都能保证(1)似然函数值增加,并且收敛到一个局部极大值1,2。y, Hrr, Ur(2)算法的命名是因为算法的每一迭代包括两个步:第一其中初始状态x0,随机误差项a和v为独立高斯分步求期望( Expectation Step),称为E步;第二步求极大布,三者之间是相互独立:xo~N(x00,P0),at计算基于不完全数据的极大似然估计要用来N0,Q),-iN0,R),E【u]=0,E[uxo]值( Maximization Step),称为M步。EM算法=0,E[ax0]=0。假设模型噪声2状态的初始条件都服从高斯2EM算法的状态空间模型参数估计分布,即:下面推导基于 Kalman平滑和算法的状态空间模p(ro| 8)=(2x)2 I Poro I-to(x0型参数估计方法6-10。给出一般形式的线性状态空)P((3)收稿日期:2008-12-26;修回日期:200-03-24基金项目:国家自然科学基金项目(60472065)dH中国煤化工∞p}-号(CNMHG(4)作者简介:王爱平(1956-),女安徽合肥人,教授,主要从事计算机教学与研究。p(x!1x1,0)=(2)号1Qs-(x第9期王爱平等:EM算法研究与应用109Fx-1)Q1(x-Fx1-1)(5) HrT)+ HPurh,])则完全数据的似然也服从高斯分布,完全数据的(10)似然可以表示为:p(x0r,y16)=p(x010)Ⅱp(x1x1,4算法实施步骤8) p(y, I xp, o)(6)现在把应用EM算法进行状态空间模型参数因此,完全数据的对数似然等于:估计的实施步骤归纳如下:LnP(Io:T, y1:T I8)=1)对参数0={R,Q,F,H,x00,Pool进行初始化In I Pono 1-4(xo-IoIo)Polo(xo-xouo)2)E步:根据:-1,执行卡尔曼滤波及平滑,计算吾Q=2x-Em)xB斯值3)M步:计算参数={R,Q,F,H,x0,P0新TInIr-1 2( -,x, yR-1( -, 3,的值;n+mN+ min(2r)4)重复参数直至收敛14。其中,N是样本量,m是状态变量的维度,n是观测变5实验仿真量的维度。考虑如下的普通线性高斯状态1空间模型取3求解对数似然的期望的极大化根据EM算法基本原理,获得以下式子y= HI t vr而状态初值x0、噪声干扰v2和v均被假设为高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)其中参数实际取值为:F=0.9,H=0.2,Q=0.7-∑(QE[(x1-Fx-1)(x-Fx-1))R=0.2,状态初始设置为x0~N(0,10)。1S(RE[(y-Hx)(y-Hx)1)(8)依据 Kalman滤波和平滑公式,然后联合EM算法对模型参数日=|F,H,Rn,Qu}进行估计EM算法迭代次数设置为400(当然选代次数设置越大,能得到更A=2(ur a P)好的估计结果)。程序运行结果θ=10.9835,0.2049,0.2107,0B=∑(x1xt1r+P-1r)683},仿真结果如图1~3所示C=∑(x1rgx11r+P:-1r)这样,完全数据的对数似然期望可以表示为:ELInt(T, yi:TIni QI-lI R I(2x)-ftr(POOEIGCoT -Io)(Ic中国煤化工Por)-2∑(QCNMHGFBFJ)tr([( -Hrr)(y图1似然函教的期望的Q(日,-1)收敛性110计算机技术与发展第19卷用D]广州暨南大学,2007[2]张杰,阳宪惠.多变量统计过程控制[M].北京:化工工业出版社,20[3]陈志国,刘文举基于EM算法的极大似然参数估计探讨J河南大学学报:自然科学版,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参数F真实值(粗线)和估计值[6]胡寿松,王执铨,胡维礼.最优控制理论与系统[M].北京科学出版社,2005[7]郭雷,程代展,冯德兴控制理论导论[M]北京:科学出版社,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]贾沛璋朱征桃最优估计及其应用[M]北京科学出版社,1994[10] Parzen E On the estimation of a probability densityfunction andmode [J]. annals of Mathematical图3参数H的收效过程Statistics,1962,33:1065-1076从上面的仿真图形可以看出,EM算法用于状态[1]waAP, Wang H. Minimising entropy and mean tracking空间模型参数估计,其估计结果还是不错的,可以为时control for affine nonlinear and non Gaussian dynamic变不同的参数模型估计提供很好的帮助。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结束语ble Stochastic System via residual Probability Density Func-针对传统的极大似然参数估计方法在解决实际问ion Shaping[c]//Proceedings of SICE 2005 Annual Confer-题时所存在的不足,引入了EM算法。讨论了算法原ence.[sL]:[sn.],2005:215-219理和实施步骤,并把EM算法推广到状态空间模型的13]任光张均东时间序列状态空间模型建模及应用(M参数估计问题,然后推导了基于Kamn平滑和EM算大连:大连海事大学出版社,2000102-112法的状态空间模型参数估计的过程,最后给出实验仿14cuL. Wang H. Mininum entropy filtering for multivariatestochastic systems with non-Gaussian noises [J].IEEETransactions on Automatic Control, 2006, 51(4): 670-695[15]梁应敞张贤达李衍达,等非高斯相关噪声中高斯信号参考文献:的时延估计[门]电子科学学刊,1997,19(5):607-612[1]陈学华状态空间模型理论与算法及其在金融计量中的应多史国计算机学会余我史《计算机技术与发展欢迎订阅VLl中国煤化工CNMHG