EM算法理论及其应用

1、2009年11月第15卷第4期安庆师范学院学报(自然科学版Journal of Anqing Teachers College(Natural Science EditionNov.2009杨基栋(华东师范大学金融与统计学院,上海200062摘要:EM算法是一种迭代算法,主要用来计算后验分布的众数或极大似然估计,广泛地应用于缺损数据、截尾数据、成群数据、带有讨厌参数的数据等所谓的不完全数据的统计推断问题。在介绍EM算法的基础上,针对EM算法收敛速度慢的缺陷,具体讨论了加速EM算法:EMB算法和MEMB算法;针对EM算法计算的局限性,给出了EM算法的推广:GEM和MCEM算法。最后给出了EM的实

2、值实例,结果精确。关键词:EM算法;极大似然估计;GEM算法;MCEM算法;EMB算法;MEMB算法中图分类号:O212.8文献标识码:A文章编号:1007-4260(200904-0030-060引言在统计领域里,统计计算技术近年来发展很快,它使许多统计方法,尤其是Bayes统计得到广泛的运用。Bayes计算方法有很多,大体上可分为两大类:一类是直接应用于后验分布以得到后验均值或后验众数的估计,以及这种估计的渐进方差或其近似;另一类算法可以总称为数据添加算法,这是近年发展很快而且应用很广的一种算法,它是在观测数据的基础上加一些“潜在数据”,从而简化计算并完成一系列简单的极大化或模拟,该“潜在

3、数据”可以是“缺损数据”或未知参数。其原理可以表述如下:设我们能观测到的数据是Y,关于Y的后验分布p(|Y很复杂,难以直接进行各种统计计算,假如我们能假定一些没有能观测到的潜在数据Z为已知(譬如,Y为某变量的截尾观测值,则Z为该变量的真值,则可能得到一个关于的简单的添加后验分布p(|Y,Z,利用p(|Y,Z的简单性我们可以对Z的假定作检查和改进,如此进行,我们就将一个复杂的极大化或抽样问题转变为一系列简单的极大化或抽样。EM算法就是一种常用的数据添加算法。1EM算法及其理论先考虑一个简单情形。设某元件的失效时间Y关于某个变量x有直线回归关系,假设在一次试验中得到一批观测数据,见右图,“

4、5;”表示该种元件的失 效时间对应的值,“”对应元件的(右截尾时间(比实际失效时间要小。如果直线斜率和截距的估计值已知,则我们可以在真实数据不小于截尾数据的前提下,将各个被截尾的失效时间估计出来(譬如,若E(Y|x>Z,则用E(Y|x作为真实数据,否则取Z作为真实数据,从而得到所谓的“完全数据”,由此“完全数据”,我们可以对直线斜率和截距进行估计(如极大似然估计,估计出新的斜率和截距后,在重新估计各个被截尾的失效时间,得到新的完全数据,如此重复,我们将一个复杂的估计时间替换成一系列简单的估计问题。将之一般化,我们可以给出EM算法。EM算法是一种迭代方法,最初由Demp ster等于197

5、7年首次提出,主要用来计算后验分布的众数或极大似然估计。近十年来引起了统计学家们的极大兴趣,在统计领域得到广泛应用。这种方法可以广3收稿日期:2009-05-21作者简介:杨基栋,女,安徽安庆人,华东师范大学金融与统计学院助理工程师。泛的应用于缺损数据,截尾数据,成群数据,带有讨厌参数的数据等所谓的不完全数据。它的每一次迭代有两步组成:E 步(求期望和M 步(极大化。一般的,以p (|Y 表示的基于观测数据的后验分布密度函数,称为观测后验分布,p (|Y ,Z 表示添加数据Z 后得到的关于的后验分布密度函数,称为添加后验分布,p (Z |,Y 表示在给定和观测数据Y 下潜在数据Z 的条件分布密

6、度函数。我们的目的是计算观测后验分布p (|Y 的众数,于是,EM 算法如下进行。记(i 为第i +1次迭代开始时后验众数的估计值,则第i +1次迭代的两步为E 步:将p (|Y ,Z lo g p (|Y ,Z 后关于Z 的条件分布求期望,从而把Z 积掉,即Q (|(i ,Y E z lo g p (|Y ,Z |(i ,Y =log p (|Y ,Z p (Z |(i ,Y d Z (1M 步:将Q (|(i ,Y 极大化,即找一个点(i+1,使Q (i+1|(i ,Y =max Q (|(i ,Y (2如此形成了一次迭代(i (i+1。将上述E 步和M 步进行迭代直至(i+1-(i 或Q

7、 (i+1|(i ,Y -Q (i |(i ,Y 充分小时停止。例1假设一次试验可能有四个结果,其发生的概率分别为125,18,20,34。此处观测数据为Y =(y 1,y 2,y 3,y 4=(125,18,20,34,取的先验分布(为(0,1上均匀分布,则的观测后验分布为p (|Y (p (Y |=(12+4y 114(1-y 214(1-y 3(4y 4(2+y 1(1-y 2+y 3y 4(3假设第一种结果可以分解成两部分,其发生概率分别为12和4,令Z 和y 1-Z 分别表示试验中结果落入这两部分的次数(Z 是不能观测的潜在数据,则的添加后验分布为p (|Y ,Z (p (Y ,Z

8、|=(12z (4y 1-z 14(1-y 214(1-y 3(4y 4y 1-z +y 4(1-y 2+y 3(4用(3求的后验众数比较麻烦,而用(4求后验众数则十分简单。在第i +1次迭代中,假设有估计值(i ,则可通过E 步和M 步得到一个新的估计,在E 步中有:Q (|(i ,Y =E z (y 1-Z +y 4log +(y 2+y 3lo g (1-|(i ,Y =y 1-E z (Z |(i ,Y +y 4lo g+(y 2+y 3log (1-因在(i 和Y 给定下,Z b (y 1,2/(i +2,故E z (Z |(i ,Y =2y 1/(i +2在M 步中,我们将Q (|

9、(i ,Y 对求导并令其为0,有(i+1=y 1+y 4-E z (Z |(i ,Y y 1+y 2+y 3+y 4-E z (Z |(i ,Y =(i y 1+(i +2y 4(i y 1+(i +2(y 2+y 3+y 4=159(i +68197(i +144(5(5式给出了有EM 算法得到的迭代公式。由此公式进行迭代,可得到观测后验分布的众数。由(5式不用迭代也可求出的估计。事实上,它是一个迭代公式,若收敛到=(159+68/(197+144,解之有=0.626821497。2EM 算法的收敛性2.1EM 算法的收敛性定理EM 算法的最大优点是简单和稳定。EM 算法的主要目的是提供一个

10、简单的迭代算法来计算极大似然估计,人们自然会问,如此建立的EM 算法是否达到预期的要求,就是说,由EM 算法得到的估计序列是否收敛,如果收敛,其结果是p (|Y 的最大值或局部最大值。为此我们给出以下两个定理。记EM 算法得到的估计序列为(i ,i =1,2,(|Y =lo g p (|Y 。定理1EM 算法在每一次迭代后均提高(观测后验密度函数值,即p (i+1|Y p (i |Y (6证明由全概率公式p (,Z |Y =p (Z |,Y p (|Y =p (|Y ,Z p (Z |Y 将上式后面两项取对数有lo g p (|Y =log p (|Y ,Z -log p (Z |,Y +lo

11、g p (Z |Y (7设现有估计(i ,将上式对Z 关于p (Z |(i ,Y 求期望,有log p (|Y =lo g p (|Y ,Z -log p (Z |,Y +log p (Z |Y p (Z |(i ,Y d Z Q (|(i ,Y -H (|(i ,Y +K (i ,Y (8其中Q (|(i ,Y 已在(1中定义,H (|(i ,Y =log p (Z |,Y p (Z |(i ,Y d Z ,K (i ,Y =log p (Z |Y p (Z |(i ,Y d Z13第4期杨基栋:EM 算法理论及其应用分别在(8中取为(i 和(i+1并相减,有log p (i+1|Y -lo

12、 g p (i |Y =Q (i+1|(i ,Y -Q (i |(i ,Y -H (i+1|(i ,Y -H (i |(i ,Y 由J ensen 不等式,EZ|(i ,y log (p (Z |(i ,Y p (Z |(i ,Y log E Z|(i ,y (p (Z |(i+1,Y p (Z |(i ,Y =0,故H (i+1|(i ,Y -H (i |(i ,Y 0,而(i+1是使Q (|(i ,Y 达到最大的,显然Q (i+1|(i ,Y -Q (i |(i ,Y 0。这就说明,用EM 算法求出的(i 的确使其对数似然L (i |Y 关于i 单增,但只有在严格单增的情况下,才能保证EM

13、 算法求出的(×为的极大似然估计。这是从纯数学的,实际应用有一定的困难。定理2(1如果p (|Y 有上界,则L (i |Y 收敛到某个L ×。(2如果Q (|关于和都连续,则在关于L 的很一般的条件下,由EM 算法得到的估计序列(i 的收敛值×是L 的稳定点。证明(1的证明由单调收敛定理立得;(2的证明见文献4。关于定理2,我们指出,定理的条件在大多数场合是满足的,定理的收敛性结论是针对(对数后验密度函数数值给出的,而后验密度函数值序列的收敛性比估计序列本身的收敛性更具意义。另外,在定理2的条件下,EM 算法的结果只能保证收敛到后验密度函数的稳定点,并不能保证收敛

14、到极大值点,事实上,任何一种算法都很难保证其结果为极大值点。较可行的办法是选取几个不同的初值进行迭代,然后在诸估计间加以选择,这可以减轻初值选取对结果的影响。2.2EM 算法收敛速度的探讨EM 算法是一种求参数极大似然估计的迭代算法,在处理不完全数据中有重要应用。EM 算法实现简单,数值计算稳定,存储量小,并具有良好的全局收敛性。但是,EM 算法收敛速度相当慢,只是次线性的收敛速度,这个缺点防碍了EM 算法的应用。现已提出了多种加速EM 算法收敛的方法,其中使用非线性规划中的Broyden 对称秩1校正公式,提出了一种加速EM 算法收敛的方法。与其他加速收敛方法相比,本方法简便易行,不必对不完

15、全数据的似然函数一维搜索,在收敛性方面,与EM 算法一样具有全局收敛性。数值试验结果表明,本方法的收敛速度比EM 算法快的多。一般地,在EM 算法的M 步求(i+1时,可求解方程组5Q (|(i 5=0,得到参数的迭代公式(i+1=G (i 。这种迭代公式通常使EM 算法的收敛速度很慢,为加速收敛,可考虑使用其他方法求解。根据著名的Fisher 公式5Q (|(i 5|=(i =g (i ,这里g (i 是不完全数据对数似然函数L (在(i 处的梯度g (i =L (|=(i ,于是问题转化为求(×,使g (×=0,从而可以使用非线性规划中的有效方法求解,达到加速收敛的目的

16、。在非线性规划的求解方法中,Broyden 对称秩1校正公式具有二次终止性和超线性的收敛速度,使用它可望改善EM 步长的缓慢变化,以加速EM 算法的收敛。但是,如果不对L (进行一维搜索,Broyden 对称秩1校正公式仅具有局部收敛性,因此考虑将Broyden 对称秩1校正公式的二次终止性和EM 算法良好的全局收敛性相结合,形成一种能够加速EM 算法收敛的新方法EM B 算法。在迭代初期,使用EM 算法,而当接近最优解时,EM 步长的变化极为缓慢,这时使用Broyden 对称秩1校正公式进行校正。算法1(EMB 算法1取初始值(0,常数(0,1,H (0=I ,r (0=g (0,k =0。

17、2若g (k =0,停止计算,(k 即为所求;否则转3。3若g (k <r (k ,则取r (k+1满足g (k ÷r (k+1<r (k ,置(k+1=(k -H (k g (k ,转5,否则转4。4置r (k+1=r (k ,取EM 步长(k ,置(k+1=(k +(k 。5对H (k 使用Broyden 对称秩1校正公式,置H (k+1=H (k +(k -H (k (k (k -H (k (k T (k -H (k (k T (k ,其中k =g (k+1-g (k ,(k =(k+1-(k 。6置k =k +1,转2。关于EMB 算法的收敛性有如下结论:定理3在

18、EM 算法收敛的条件下,EMB 算法是收敛的。证明分三种情形讨论23安庆师范学院学报(自然科学版2009年1若在EMB 算法中,每次迭代都取EM 步长(k 计算(k+1,即EMB 算法即为EM 算法。2若在EMB 算法中,每次迭代都按公式(k+1=(k -H (k g (k 计算(k+1,则由g (k ÷r (k+1<r (k 可得r (k+1<r (k <2r (k-1<<k+1r (0,因<1,故k 时,记(k (×,则g (×=0,因此EMB 算法收敛。3一般情形时,若设EMB 算法局部收敛,即存在某个>0及k 1,

19、当k >k 1时,g (k >,而由非负单调递减数列r (k 的构造可知,存在子列r (k ,r (k 0。因此存在k 2,当k >max k 1,k 2时,r (k <,于是g (k <r (k <(t 1,矛盾。因此EMB 算法收敛。从数值试验结果来看,EMB 算法能显著减少迭代次数,但有些时候效果不好,为此考虑在EMB 算法中加进一个延迟参数,目的是将校正尽量靠后,当接近最优解时,再用Broyden 对称秩1校正公式。算法2(M EMB 算法1取初始值(0,常数(0,1,非负整数,r (0=g (0,k =0。2若g (k =0,停止计算,(k 即为所

20、求;否则转3。3若k >,转4,否则转5。4若g (k <r (k ,则取r (k+1满足g (k ÷r (k+1<r (k ,置(k+1=(k -H (k g (k ,转6,否则转5。5置r (k+1=r (k ,取EM 步长(k ,置(k+1=(k +(k 。6对H (k 使用Broyden 对称秩1校正公式,置H (k+1=H (k +(k -H (k (k (k -H (k (k T (k -H (k (k T (k ,其中k =g (k+1-g (k ,(k =(k+1-(k 。7置k =k +1,转2。2.3EM 算法的推广算法EM 算法得到广泛应用的一

21、个重要原因是在M 步中,求极大化的方法与完全数据下求极大化的方法完全一样。在许多场合,这样的极大化有显式表示,然而并不总是这样,有时要找Q (|(i ,Y 达到最大的是困难的,一个比较简单的方法是找一个(i+1,使Q (i+1|(i ,Y >Q (i |(i ,Y (9由(1和(9组成的EM 算法称为广义EM 算法(GEM 算法。我们可以利用无约束最优化方法中的一些方法,像最速下降法,Newto n 法,共轭方向法和共轭梯度法,拟Newto n 法,Powell 方向加速法,来构造广义EM 算法(GEM 算法,既可避免EM 算法M 步中表达式的局限,又可加速EM 算法,使EM 算法的收敛

22、速度得到很大的提高。在EM 算法的M 步求(i+1时,根据著名的Fisher 公式5Q (|(i 5|=(i =g (i ,这里g (i 是不完全数据对数似然函数L (在(i 处的梯度,于是问题转化为求(3,使g (3=0,从而可以使用非线性规划中的有效方法求解,达到加速收敛的目的。像上一章介绍的EMB 算法和M EMB 算法,就是两种广义的EM 算法,及GEM 算法。还有Laird 等利用数值分析中的外推Ait ken 加速来加速EM 算法,Lo uis 在Newto n 法的基础上关于EM 算法提出了一种加速的方法,Amshidian 等提出了EM 算法的广义共轭梯度加速方法等。而对于EM

23、 算法的E 步,有时要获得期望的显式表示是不可能的,即使近似计算也很困难,这时用Mo nte Carlo 方法来完成,就是所谓的Mo nte Carlo EM (MCEM 方法,它将E 步改为:(E1由p (Z |(i ,Y 随机的抽取m 个随机数Z 1,Z m ;(E2计算Q (|(i ,Y =1m m j =1log p (|Z j ,Y 。由大数定理,只要m 足够大,Q (|(i ,Y 与(|(i ,Y 很接近,从而我们可以在M 步中对Q (|(i ,Y 求极大化。在MCEM 算法中有两点是主要考虑的:一是m 的大小的确定,从精确角度来讲,m 自然越大越好,但过大的m 使得计算的效率太低

24、,一般在开始时m 不需要很大;另一点是收敛性的判断,因在E 步中是采用的Mo nte Carlo 方法,若要求这样得到的(i 收敛到一点显然是不现实的。在MCEM 中,收敛性的判别往往可借助图形来进行。若经过若干次迭代后,迭代值围绕直线=3小幅波动,则可以认为算法收33第4期杨基栋:EM 算法理论及其应用敛了,此时,为增加估计精度,可增加m 的值再运行一段时间,就可停止。MCEM 算法比较灵活,但是需要仔细选择模拟容量和确保正确的收敛性准则。我们可以通过增加迭代次数来提高模拟容量,例如,在McColloch (1994方法中,模拟容量随着迭代次数线性增加。除此以外,由于蒙特卡罗误差,该EM 算

25、法不具有单调性,难以估计其收敛性。在McColloch 方法中,MCEM 算法通过规定迭代的次数来结束算法。Shi 和Lee (2000用临时抽样的方法直观的观察MCEM 算法的收敛性。3EM 算法的应用实例EM 算法可以应用于医学研究中,尤其是临床医学中十分常见的一种数据观测形式为重复观测,其特点是在同一实验单位上进行多次重复观测,这个过程由于各种原因经常导致实验观测数据缺失,如动物的意外死亡,记录仪器发生故障,被调查者拒绝回答相关调查项目等,然而,目前绝大多数重复观测数据统计分析方法都有一个基本假定,即每个实验单位具有完全的重复观测数据(通常假定服从多元正态分布。因此有必要对缺失数据进行适

26、当的处理,否则数据分析时,通常的做法是删去具有缺失的部分观察记录而不考虑记录数据所蕴涵的信息,从而造成信息的损失以及分析结果的偏性。缺失数据的处理方式在很大程度上依赖医学数据观测的实际背景,例如,常规的实验设计得到的独立结构数据方差分析前的缺项估计的原理,都是解误差离均差平方和对缺项变量的导数等于0的方程,从而获得缺项的估计值。而医学重复观测设计得到的数据有别于独立结构数据,其数据间存在一定的相关性,因此需要进行专门的缺项估计方法的研究。为此。我们在软件编程的基础上,应用EM 算法进行了这方面的数据处理。目前,可用于进行缺项估计的方法很多,但EM 算法较为优秀,它是利用所有的资料信息来进行缺项

27、估计,以“补缺”的方式将含有缺失值的“不完全资料”转化为“完全资料”,因此,对EM 算法得到的“完全资料”可以应用所有的重复观测资料分析方法进行处理分析。EM 方法补缺使估计系数的方差明显变小,从而提高了生长曲线系数的估计精度。当医院科研数据量较大时,为方便缺失值的估计,就需用MA TL AB 进行编程,另外,在实际工作中为了了解多个指标间的关系及变化规律,常常需要同时观测个体的多个反应指标,此时,虽然数据不是重复观测数据,但如果缺失数据随机发生,即缺失数据发生的可能性不被该数据值的大小所影响,则同样可用EM 算法进行估计,因此EM 算法不但适用于重复观测缺失数据的处理,而且适用于一般情况下多变