EM算法理论及其应用

2009年11月第15卷第4期

安庆师范学院学报(自然科学版)

Jou rn al of A nqing Teachers College(Natural S cience E dition)

Nov.2009

Vol.15No.4 EM算法理论及其应用

杨基栋

(华东师范大学金融与统计学院,上海200062)

摘 要:E M算法是一种迭代算法,主要用来计算后验分布的众数或极大似然估计,广泛地应用于缺损数据、截尾数

据、成群数据、带有讨厌参数的数据等所谓的不完全数据的统计推断问题。在介绍EM算法的基础上,针对E M算法收

敛速度慢的缺陷,具体讨论了加速EM算法:EM B算法和M EM B算法;针对EM算法计算的局限性,给出了EM算法的

推广:GEM和M CEM算法。最后给出了EM的实值实例,结果精确。

关键词:E M算法;极大似然估计;GEM算法;M CE M算法;EM B算法;M EM B算法

中图分类号:O212.8 文献标识码:A 文章编号:1007-4260(2009)04-0030-06

0 引 言

在统计领域里,统计计算技术近年来发展很快,它使许多统计方法,尤其是Bay es统计得到广泛的运用。Bay es计算方法有很多,大体上可分为两大类:一类是直接应用于后验分布以得到后验均值或后验众数的估计,以及这种估计的渐进方差或其近似;另一类算法可以总称为数据添加算法,这是近年发展很快而且应用很广的一种算法,它是在观测数据的基础上加一些 潜在数据,从而简化计算并完成一系列简单的极大化或模拟,该 潜在数据可以是 缺损数据或未知参数。其原理可以表述如下:设我们能观测到的数据是Y, 关于Y的后验分布p( |Y)很复杂,难以直接进行各种统计计算,假如我们能假定一些没有能观测到的潜在数据Z为已知(譬如,Y为某变量的截尾观测值,则Z为该变量的真值),则可能得到一个关于 的简单的添加后验分布p( |Y,Z),利用p( |Y,Z)的简单性我们可以对Z的假定作检查和改进,如此进行,我们就将一个复杂的极大化或抽样问题转变为一系列简单的极大化或抽样。EM算法就是一种常用的数据添加算法。

1 EM算法及其理论

先考虑一个简单情形。设某元件的失效时间Y关于某个变量x有直线回归关系,假设在一次试验中得到一批观测数据,见右图, !表示该种元件的失

效时间对应的值, ?对应元件的(右)截尾时间

(比实际失效时间要小)。如果直线斜率和截距的估

计值已知,则我们可以在真实数据不小于截尾数据

的前提下,将各个被截尾的失效时间估计出来(譬

如,若E(Y|x)>Z,则用E(Y|x)作为真实数据,

否则取Z作为真实数据),从而得到所谓的 完全数

据,由此 完全数据,我们可以对直线斜率和截距

进行估计(如极大似然估计),估计出新的斜率和截

距后,在重新估计各个被截尾的失效时间,得到新的完全数据,如此重复,我们将一个复杂的估计时间替换成一系列简单的估计问题。将之一般化,我们可以给出EM算法。

EM算法是一种迭代方法,最初由Dem pster等于1977年首次提出,主要用来计算后验分布的众数或极大似然估计。近十年来引起了统计学家们的极大兴趣,在统计领域得到广泛应用。这种方法可以广

*收稿日期:2009-05-21

作者简介:杨基栋,女,安徽安庆人,华东师范大学金融与统计学院助理工程师。

泛的应用于缺损数据,截尾数据,成群数据,带有讨厌参数的数据等所谓的不完全数据。

它的每一次迭代有两步组成:E 步(求期望)和M 步(极大化)。一般的,以p ( |Y)表示 的基于观测数据的后验分布密度函数,称为观测后验分布,p ( |Y,Z)表示添加数据Z 后得到的关于 的后验分布密度函数,称为添加后验分布,p (Z | ,Y)表示在给定 和观测数据Y 下潜在数据Z 的条件分布密度函数。我们的目的是计算观测后验分布p ( |Y)的众数,于是,EM 算法如下进行。

记 (i)为第i +1次迭代开始时后验众数的估计值,则第i +1次迭代的两步为

E 步:将p ( |Y,Z)log p ( |Y,Z)后关于Z 的条件分布求期望,从而把Z 积掉,即

Q( | (i),Y) E z [lo g p ( |Y,Z)| (i),Y]=#lo g [p ( |Y,Z)]p (Z | (i),Y)d Z (1)

M 步:将Q( | (i),Y)极大化,即找一个点 (i+1),使Q( (i+1)| (i),Y)=m ax

Q( | (i),Y)(2)如此形成了一次迭代 (i)? (i+1)。将上述E 步和M 步进行迭代直至% (i+1)- (i)%或

%Q( (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+ 4)y 1[14(1- )]y 2[14(1- )]y 3( 4)y 4?(2+ )y 1(1- )y 2+y 3 y 4(3)

假设第一种结果可以分解成两部分,其发生概率分别为12和 4,令Z 和y 1-Z 分别表示试验中结果

落入这两部分的次数(Z 是不能观测的潜在数据),则 的添加后验分布为

p ( |Y ,Z)?&( )p (Y ,Z | )=(12)z ( 4)y 1-z [14(1- )]y 2[14(1- )]y 3( 4)y 4? y 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 4)log +(y 2+y 3)log (1- )| (i),Y]=

[y 1-E z (Z | (i),Y)+y 4]log +(y 2+y 3)log (1- )

因在 (i)和Y 给定下,Z ~b(y 1,2/( (i)+2)),故E z (Z | (i),Y)=2y 1/( (i)+2)

在M 步中,我们将Q( | (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)+2)y 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。

2 EM 算法的收敛性

2.1 EM 算法的收敛性定理

EM 算法的最大优点是简单和稳定。EM 算法的主要目的是提供一个简单的迭代算法来计算极大似然估计,人们自然会问,如此建立的EM 算法是否达到预期的要求,就是说,由EM 算法得到的估计序列是否收敛,如果收敛,其结果是p ( |Y)的最大值或局部最大值。为此我们给出以下两个定理。

记EM 算法得到的估计序列为 (i),i =1,2,(( |Y)=log p ( |Y)。

定理1 EM 算法在每一次迭代后均提高(观测)后验密度函数值,即

p ( (i+1)|Y))p ( (i)|Y)(6)

证明 由全概率公式p ( ,Z |Y)=p (Z | ,Y)p ( |Y)=p ( |Y,Z)p (Z |Y)

将上式后面两项取对数有 log p ( |Y)=log p ( |Y,Z)-log p (Z | ,Y)+lo g p (Z |Y)(7)

设现有估计 (i),将上式对Z 关于p (Z | (i),Y)求期望,有

log p ( |Y)=#[log 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 Z ?31?第4期 杨基栋:EM 算法理论及其应用

分别在(8)中取 为 (i)和 (i+1)并相减,有lo g p ( (i+1)|Y)-lo g p ( (i)|Y)=

[Q( (i+1)| (i),Y)-Q( (i)| (i),Y)]-[H ( (i+1)| (i),Y)-H ( (i)| (i),Y)]

由Jensen 不等式,E Z| (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 算法求出的 (!)为 的极大似然估计。这是从纯数学的,实际应用有一定的困难。

定理2 (1)如果p ( |Y)有上界,则L ( (i)|Y)收敛到某个L !。

(2)如果Q( | )关于 和 都连续,则在关于L 的很一般的条件下,由EM 算法得到的估计序列 (i)的收敛值 !是L 的稳定点。

证明 (1)的证明由单调收敛定理立得;(2)的证明见文献[4]。

关于定理2,我们指出,定理的条件在大多数场合是满足的,定理的收敛性结论是针对(对数)后验密度函数数值给出的,而后验密度函数值序列的收敛性比估计序列本身的收敛性更具意义。另外,在定理2的条件下,EM 算法的结果只能保证收敛到后验密度函数的稳定点,并不能保证收敛到极大值点,事实上,任何一种算法都很难保证其结果为极大值点。较可行的办法是选取几个不同的初值进行迭代,然后在诸估计间加以选择,这可以减轻初值选取对结果的影响。

2.2 EM 算法收敛速度的探讨

EM 算法是一种求参数极大似然估计的迭代算法,在处理不完全数据中有重要应用。EM 算法实现简单,数值计算稳定,存储量小,并具有良好的全局收敛性。但是,EM 算法收敛速度相当慢,只是次线性的收敛速度,这个缺点防碍了EM 算法的应用。现已提出了多种加速EM 算法收敛的方法,其中使用非线性规划中的Broyden 对称秩1校正公式,提出了一种加速EM 算法收敛的方法。与其他加速收敛方法相比,本方法简便易行,不必对不完全数据的似然函数一维搜索,在收敛性方面,与EM 算法一样具有全局收敛性。数值试验结果表明,本方法的收敛速度比EM 算法快的多。

一般地,在EM 算法的M 步求 (i+1)时,可求解方程组 Q( | (i))

=0,得到参数的迭代公式 (i+1)=G ( (i))。这种迭代公式通常使EM 算法的收敛速度很慢,为加速收敛,可考虑使用其他方法求解。

根据著名的Fisher 公式 Q( | (i))

| = (i)=g( (i)),这里g( (i))是不完全数据对数似然函数L ( )在 (i)处的梯度g ( (i))=▽L ( )| = (i),于是问题转化为求 (!),使g( (!))=0,从而可以使用非线性规划中的有效方法求解,达到加速收敛的目的。

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

算法1(EMB 算法)

1)取初始值 (0),常数 ?(0,1),H (0)=I ,r (0)=%g( (0))%,k =0。

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)使用Broy den 对称秩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)。

关于EM B 算法的收敛性有如下结论:

定理3 在EM 算法收敛的条件下,EMB 算法是收敛的。证明 分三种情形讨论

?32?安庆师范学院学报(自然科学版) 2009年

1)若在EM B 算法中,每次迭代都取EM 步长. (k )计算 (k+1),即EMB 算法即为EM 算法。

2)若在EM B 算法中,每次迭代都按公式 (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,因此EM B 算法收敛。

3)一般情形时,若设EM B 算法局部收敛,即存在某个?>0及k 1,当k >k 1时,%g( (k))%>?,而

由非负单调递减数列{r (k )}的构造可知,存在子列{r (k )},r (k )?0。因此存在k 2,当k >max {k 1,k 2}时,

r (k)0

从数值试验结果来看,EM B 算法能显著减少迭代次数,但有些时候效果不好,为此考虑在EM B 算法中加进一个延迟参数#,目的是将校正尽量靠后,当接近最优解时,再用Broy den 对称秩1校正公式。

算法2(MEM B 算法)

1)取初始值 (0),常数 ?(0,1),非负整数#,r (0)=%g( (0))%,k =0。

2)若%g( (k ))%=0,停止计算, (k)即为所求;否则转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)使用Broy den 对称秩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.3 EM 算法的推广算法

2.3.1 GEM 算法

EM 算法得到广泛应用的一个重要原因是在M 步中,求极大化的方法与完全数据下求极大化的方

法完全一样。在许多场合,这样的极大化有显式表示,然而并不总是这样,有时要找Q( | (i),Y)达到最

大的 是困难的,一个比较简单的方法是找一个 (i+1),使

Q( (i+1)| (i),Y)>Q( (i)| (i),Y)(9)

由(1)和(9)组成的EM 算法称为广义EM 算法(GEM 算法)。

我们可以利用无约束最优化方法中的一些方法,像最速下降法,New ton 法,共轭方向法和共轭梯度法,拟N ew ton 法,Po w ell 方向加速法,来构造广义EM 算法(GEM 算法),既可避免EM 算法M 步中表达式的局限,又可加速EM 算法,使EM 算法的收敛速度得到很大的提高。

在EM 算法的M 步求 (i+1)时,根据著名的Fisher 公式 Q( | (i))

| = (i)=g( (i)),这里g( (i))是不完全数据对数似然函数L ( )在 (i)处的梯度,于是问题转化为求 (*),使g ( (*))=0,从而可以使用非线

性规划中的有效方法求解,达到加速收敛的目的。

像上一章介绍的EMB 算法和M EMB 算法,就是两种广义的EM 算法,及GEM 算法。还有Laird 等利用数值分析中的外推Aitken 加速来加速EM 算法,Louis 在N ew ton 法的基础上关于EM 算法提出了一种加速的方法,Am shidian 等提出了EM 算法的广义共轭梯度加速方法等。

2.3.2 M CEM 算法

而对于EM 算法的E 步,有时要获得期望的显式表示是不可能的,即使近似计算也很困难,这时用Mo nte Carlo 方法来完成,就是所谓的M onte Car lo EM (M CEM )方法,它将E 步改为:

(E1)由p (Z | (i),Y)随机的抽取m 个随机数Z 1,(,Z m ;

(E2)计算Q ^( | (i),Y)=1m 1m j =1

log p ( |Z j ,Y)。由大数定理,只要m 足够大,Q ^( | (i),Y)与( | (i),Y)很接近,从而我们可以在M 步中对Q ^( | (i),Y)

求极大化。

在MCEM 算法中有两点是主要考虑的:一是m 的大小的确定,从精确角度来讲,m 自然越大越好,但过大的m 使得计算的效率太低,一般在开始时m 不需要很大;另一点是收敛性的判断,因在E 步中是采用的Mo nte Carlo 方法,若要求这样得到的 (i)收敛到一点显然是不现实的。在M CEM 中,收敛性的判别往往可借助图形来进行。若经过若干次迭代后,迭代值围绕直线 = *小幅波动,则可以认为算法收?

33?第4期 杨基栋:EM 算法理论及其应用

敛了,此时,为增加估计精度,可增加m 的值再运行一段时间,就可停止。

MCEM 算法比较灵活,但是需要仔细选择模拟容量和确保正确的收敛性准则。我们可以通过增加迭代次数来提高模拟容量,例如,在M cCo lloch(1994)方法中,模拟容量随着迭代次数线性增加。

除此以外,由于蒙特卡罗误差,该EM 算法不具有单调性,难以估计其收敛性。在McColloch 方法中,M CEM 算法通过规定迭代的次数来结束算法。Shi 和Lee(2000)用临时抽样的方法直观的观察MCEM 算法的收敛性。

3 EM 算法的应用实例

EM 算法可以应用于医学研究中,尤其是临床医学中十分常见的一种数据观测形式为重复观测,其特点是在同一实验单位上进行多次重复观测,这个过程由于各种原因经常导致实验观测数据缺失,如动物的意外死亡,记录仪器发生故障,被调查者拒绝回答相关调查项目等,然而,目前绝大多数重复观测数据统计分析方法都有一个基本假定,即每个实验单位具有完全的重复观测数据(通常假定服从多元正态分布)。因此有必要对缺失数据进行适当的处理,否则数据分析时,通常的做法是删去具有缺失的部分观察记录而不考虑记录数据所蕴涵的信息,从而造成信息的损失以及分析结果的偏性。

缺失数据的处理方式在很大程度上依赖医学数据观测的实际背景,例如,常规的实验设计得到的独立结构数据方差分析前的缺项估计的原理,都是解误差离均差平方和对缺项变量的导数等于0的方程,从而获得缺项的估计值。而医学重复观测设计得到的数据有别于独立结构数据,其数据间存在一定的相关性,因此需要进行专门的缺项估计方法的研究。为此。我们在软件编程的基础上,应用EM 算法进行了这方面的数据处理。

目前,可用于进行缺项估计的方法很多,但EM 算法较为优秀,它是利用所有的资料信息来进行缺项估计,以 补缺 的方式将含有缺失值的 不完全资料 转化为 完全资料 ,因此,对EM 算法得到的 完全资料 可以应用所有的重复观测资料分析方法进行处理分析。EM 方法补缺使估计系数的方差明显变小,从而提高了生长曲线系数的估计精度。

当医院科研数据量较大时,为方便缺失值的估计,就需用MAT LAB 进行编程,另外,在实际工作中为了了解多个指标间的关系及变化规律,常常需要同时观测个体的多个反应指标,此时,虽然数据不是重复观测数据,但如果缺失数据随机发生,即缺失数据发生的可能性不被该数据值的大小所影响,则同样可用EM 算法进行估计,因此EM 算法不但适用于重复观测缺失数据的处理,而且适用于一般情况下多变量缺失数据的处理分析,其适用范围更广泛。

以下的具体实例中,我们将介绍EM 算法在医学研究中的应用。例2 有8名受试者用庆大霉素80mg 后的血药浓度(m g ?L -1)动态变化观察结果如表1所示,其中有2处数据存在缺失:第4例120

min 时观察缺失和第8例50min 时

观察缺失。

设Y 1,Y 2,(,Y n 为p 维正态分

布总体N p (?,1)的随机样本,完

表1 受试者8名注射药物后血药浓度%/(mg ?L -1)subject 10min 30m in 50min 90min 120min 150min 180min

1 1.5158 3.6541 3.4466 2.9308 2.7925 1.93080.8962

2 2.6006 4.4290 4.6864 3.0858 3.2574 2.0858 1.1287

3 2.8778 4.6270 5.0032 3.125

4 3.3762 2.1254 1.1881

40.8844 3.2031 2.7250 2.8406N AN 1.84060.760950.9913 3.2795 2.8472 2.8559 2.5677 1.85590.78396 2.3998 4.2856 4.4569 3.0571 3.1713 2.0571 1.08577 1.7201 3.8001 3.6801 2.9600 2.8800 1.96000.94008 2.4830 4.3450NA N 3.06903.20702.06901.1035全数据的两个充分统计量为T 1=

1n k =1Y k 和T 2=1n k =1Y k Y 0k ,假定?和1均为未知,则缺失的估计需要进行

一下迭代计算:

1)给出?和1的初始估计值?~和1~可以用未缺失的数据计算得到的算术均数作为?的初始值,用算

术均数代替缺失值后,可求得1~的极大似然估计初始值。

2)计算条件数学期望对于具有缺失值的每一个向量Y k ,令Y (1)k 表示缺失的分量,Y (2)k 表示已知的分

量,如果对?~和1~进行相应的分块变化,则Y (1)k 的条件正态分布的均值为Y ~(1)k =E(Y (1)k ∣Y (2)k ,?~,1~)=?~(1)+1~121~-122(Y (2)k -?~

(2)),Y (1)k Y (1)0k 和Y (1)k Y (2)0k 的条件正态分布的均值分别为Y (1)k Y (1)0k =E(Y (1)k Y (1)0k |Y (2)k ;?~,1~)=1~11-1~121~-1221~21+Y ~(1)k Y ~(1)0k

Y (1)k Y (2)0k =E(Y (1)k Y (2)0k |Y (2)k ;?~,1~)=Y (1)k Y (2)0k 将具有缺失分量Y k 的结果与样本实际观察数据的计算结果结合起来,就可得到充分的估计值T ~1和T ~2。?34?安庆师范学院学报(自然科学版) 2009年

3)计算?和1的最大似然估计校正值?~和1~:?~=T ~1n ,1~=T ~2n -?~?~0。

4)重复以上的2-3步,直至?~和1~收敛为止。

经过EM 迭代算法可得:

?~=[1.9341 3.9529 3.9247 2.9906 2.9717 1.9906 0.9859]

1~=0.5081 0.3629 0.5807 0.0726 0.2187 0.0726 0.1089

0.3629 0.2592 0.4148 0.0518 0.1566 0.0518 0.0778

0.5807 0.4148 0.6636 0.0830 0.2489 0.0830 0.1244

0.0726 0.0518 0.0830 0.0104 0.0311 0.0104 0.0156

0.2187 0.1556 0.2489 0.0311 0.0933 0.0311 0.0467

0.0726 0.0518 0.0830 0.0104 0.0311 0.0104 0.0156

0.1089 0.0778 0.1244 0.0156 0.0469 0.0156 0.0233

迭代算法补入的两个数据:第4行第5个为2.5216,第8行第3个为4.5522。

5)M onte Carlo 模拟,随机取x ~N 4(?,1),样本容量n =4,其中?~=[1234]

1~=0.2500 0.3000 0.3500 0.4000

0.3000 0.3600 0.4200 0.48000.3500 0.4200 0.4900 0.5600

0.4000 0.4800 0.5600 0.6400

去掉其中第2行第3列和第3行第2列的两个数。(1)如果删去,不考虑这部分数据蕴含的信息,得到?1=[0.7894 1.2284 2.2466 3.6630]

11=0.1469 0.0344 0.4909 0.2350

0.0344 0.6786 0.1527 0.05510.4909 0.1527 1.7173 0.7854

0.2350 0.0551 0.7854 0.3760

(2)如果用EM 算法计算,得到?2=[0.7894 1.5158 2.9236 3.6630]

11=0.1469 0.1031 0.0697 0.2350

0.1031 0.0887 0.0318 0.16490.0697 0.0318 0.0512 0.1115

0.2350 0.1649 0.1115 0.3760

通过计算比较可得:no rm (?1-?~)=1.1493,no rm (11-1~)=1.3936,norm(11-1~)=0.6311,

no rm (?2-?~)=1.2333

我们发现通过EM 算法计算得到的结果比直接删除得到的结果更接近精确值。

参考文献:

[1]茆诗松,王静龙,等.高等数理统计[M ].北京:高等教育出版社,1998:432-438.

[2]W u C F J.On th e convergence properties of the EM algorithm[J ].Th e Annals of Statistics,1983(11):95-103.

[3]高旅瑞,陈志,王家润.一种加速EM 算法收敛的方法[J].数理统计与应用概率,1998,13(4):342-348.

[4]陈长生,王彤,徐勇勇,尚磊.医学科研中缺失数据的E M 估计[J].第四军医大学学报,2002,23(1):59-61.

EM Algorithm and Its Applications

YANG Ji !dong

(Sch ool of Finance and Statis tics,East China Normal University,Shanghai 200062,Ch ina)

Abstract:EM algorithm,a method of iteration,is mainly used to calculate th e m ode of a p os terior distribution or the max imum lik elih ood es tim ate.EM algorithm has been wid ely applied to s tatistical in feren ces involving in complete data such as mis sing data,censorin g data,group data and data bearing disgusting parameters.T his thes is firstly introduces EM alg orith m.T o deal with the defects of EM algorithm 2s slow convergence sp eed,the accelerating EM algorithms,nam ely EM B algorithm and M EM B algorithm ar e in trodu ced.We also briefly introduce the tw o generaliz ed m ethods,GE M algorithm and M CEM alg orith m,to avoid its lim itations.T he th esis gives th e exam ples and M on te Carlo simula !tion s in the end.By designing M AT LAB programs w e obtain and analyze th e res ults.

K ey words:EM algorithm ,max imum lik elih ood es tim ate,GE M algorithm ,M CEM algorithm,EM B alg orith m,M EM B algorith m ?

35?第4期 杨基栋:EM 算法理论及其应用