EM算法及应用实例

统计学模拟及其R实现论文报告论文（设计）题目：EM算法及应用实例学院：数学与统计学院专业：2014级统计学任课老师：黄恒振姓名：钟思敏学号：201410700258EM算法及应用实例2EM算法及应用实例【内容摘要】EM算法是一种迭代算法,主要用来计算后验分布的众数或极大似然估计,广泛地应用于缺损数据、截尾数据、成群数据、带有讨厌参数的数据等所谓的不完全数据的统计推断问题。在介绍EM算法的基础上,针对EM算法收敛性定理,给出了EM的实值实例,最后给出EM算法的医学应用。【关键词】Bayes；后验分布；极大似然估计；医学实例【Abstract】TheEMalgorithmisaniterativealgorithm,isusedtocalculatetheposteriordistributionofthenumberormaximumlikelihoodestimation,widelyusedinthedefectdata,ofaccordingto,groupsofdata,withhateparameterssuchastheso-calledincompletedatastatisticalinferenceproblem.BasedontheintroductionoftheEMalgorithm,inviewoftheEMalgorithmconvergencetheorem,givesarealvalueofEMinstance,themedicalapplicationofEMalgorithmisgiven.【Keywords】Bayestheorem；Posteriordistribution；Maximumlikelihoodestimation；MedicalinstanceEM算法及应用实例3目录【内容摘要】....................................................2【关键词】......................................................20引言..........................................................41EM算法及其理论..............................................42EM算法的收敛性..............................................62.1EM算法的收敛性定理.................................62.2EM算法收敛速度的探讨...............................83EM算法的应用实例............................................9参考文献:......................................................11EM算法及应用实例40引言在统计领域里,统计计算技术近年来发展很快,它使许多统计方法,尤其是Bayes统计得到广泛的运用。Bayes计算方法有很多,大体上可分为两大类:一类是直接应用于后验分布以得到后验均值或后验众数的估计,以及这种估计的渐进方差或其近似;另一类算法可以总称为数据添加算法,这是近年发展很快而且应用很广的一种算法,它是在观测数据的基础上加一些“潜在数据”,从而简化计算并完成一系列简单的极大化或模拟,该“潜在数据”可以是“缺损数据”或未知参数。其原理可以表述如下:设我们能观测到的数据是Y,θ关于Y的后验分布p(θ|Y)很复杂,难以直接进行各种统计计算,假如我们能假定一些没有能观测到的潜在数据Z为已知(譬如,Y为某变量的截尾观测值,则Z为该变量的真值),则可能得到一个关于θ的简单的添加后验分布p(θ|Y,Z),利用p(θ|Y,Z)的简单性我们可以对Z的假定作检查和改进,如此进行,我们就将一个复杂的极大化或抽样问题转变为一系列简单的极大化或抽样。EM算法就是一种常用的数据添加算法。1EM算法及其理论EM（Expectation-Maximization）算法是一种迭代方法,最初由Dempster等于1977年首次提出,主要用来计算后验分布的众数(即极大似然估计)。近十年来引起了统计学家们的极大兴趣,在统计领域得到广泛应用。这种方法可以广泛的应用于缺损数据,截尾数据,成群数据,带有讨厌参数的数据等所谓的不完全数据。它的每一次迭代有两步组成:E步(求期望)和M步(极大化)。一般地,以p(θ|Y)表示θ的基于观测数据的后验分布密度函数,称为观测后验分布,p(θ|Y,Z)表示添加数据Z后得到的关于θ的后验分布密度函数,称为添加后验分布,p(Z|θ,Y)表示在给定θ和观测数据Y下潜在数据Z的条件分布密度函数。我们的目的是计算观测后验分布p(θ|Y)的众数,于是,EM算法按下述步骤进行。记θ(i)为第i+1次迭代开始时后验众数的估计值,则第i+1次迭代的两步为EM算法及应用实例5E步:将p(θ|Y,Z)或logp(θ|Y,Z)后关于Z的条件分布求期望,从而把Z积掉,即Q(θ|θ(i),Y)Ez[logp(θ|Y,Z)|θ(i),Y]=∫log[p(θ|Y,Z)]p(Z|θ(i),Y)dZ(1)M步:将Q(θ|θ(i),Y)极大化,即找一个点θ(i+1),使Q(θ(i+1)|θ(i),Y)=maxθQ(θ|θ(i),Y)(2)如此形成了一次迭代θ(i)→θ(i+1)。将上述E步和M步进行迭代直至‖θ(i+1)-θ(i)‖或‖Q(θ(i+1)|θ(i),Y)-Q(θ(i)|θ(i),Y)‖充分小时停止。例1（Rao（1965））设有197种动物服从多项分布，将其分成四类，观测到的数据为y=(y1,y2,y3,y4)=(125,18,20,34),再设属于各类的概率分布为（+，(1-θ)，(1-θ)，），其中θ（0，1）试估计θ。取θ的先验分布Π(θ)为(0,1)上均匀分布,则θ的观测后验分布为p(θ|Y)∝Π(θ)p(Y|θ)=(+)y1[(1-θ)]y2[(1-θ)]y3()y4∝(2+θ)y1(1-θ)y2+y3θy4(3)假设第一种结果可以分解成两部分,其发生概率分别为和θ,令Z和y1-Z分别表示试验中结果落入这两部分的次数(Z是不能观测的潜在数据),则θ的添加后验分布为p(θ|Y,Z)∝Π(θ)p(Y,Z|θ)=()z()y1-z[(1-θ)]y2[(1-θ)]y3()y4∝θy1-z+y4(1-θ)y2+y3(4)EM算法及应用实例6用(3)求θ的后验众数比较麻烦,而用(4)求后验众数则十分简单。在第i+1次迭代中,假设有估计值θ(i),则可通过E步和M步得到θ一个新的估计,在E步中有:Q(θ|θ(i),Y)=Ez[(y1-Z+y4)logθ+(y2+y3)log(1-θ)|θ(i),Y]=[y1-Ez(Z|θ(i),Y)+y4]logθ+(y2+y3)log(1-θ)因在θ(i)和Y给定下,Z～b(y1,2/(θ(i)+2)),故Ez(Z|θ(i),Y)=2y1/(θ(i)+2)在M步中,我们将Q(θ|θ(i),Y)对θ求导并令其为0,有θ(i+1)=θ=θ=(5)(5)式给出了有EM算法得到的迭代公式。由此公式进行迭代,可得到观测后验分布的众数。由(5)式不用迭代也可求出θ的估计。事实上,它是一个迭代公式,若收敛到θ=(159θ(i)+68)/(197θ(i)+144),解之有θ=0.626821497。2EM算法的收敛性2.1EM算法的收敛性定理EM算法的最大优点是简单和稳定。EM算法的最大目的是提供一个简单的迭代算法来计算后验众数（或MLE极大似然估计）,我们不禁问,由EM算法得到的估计序列是否收敛？如果收敛,其结果是否为p(θ|Y)的最大值或局部最大值。为此我们给出以下两个定理。记EM算法得到的估计序列为θ(i),i=1,2,…，L(θ|Y)=logp(θ|Y)。定理1EM算法在每一次迭代后均提高(观测)后验密度函数值,即EM算法及应用实例7p(θ(i+1)|Y)≥p(θ(i)|Y)(6)证明由全概率公式（乘数公式）有p(θ,Z|Y)=p(Z|θ,Y)p(θ|Y)=p(θ|Y,Z)p(Z|Y)将上式后面两项取对数有logp(θ|Y)=logp(θ|Y,Z)-logp(Z|θ,Y)+logp(Z|Y)(7)设现有估计θ(i),将上式对Z关于p(Z|θ(i),Y)求期望,有logp(θ|Y)=∫[logp(θ|Y,Z)-logp(Z|θ,Y)+logp(Z|Y)]p(Z|θ(i),Y)dZ==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)dZ,K(θ(i),Y)=∫log[p(Z|Y)]p(Z|θ(i),Y)dZ分别在(8)中取θ为θ(i)和θ(i+1)并相减,有logp(θ(i+1)|Y)-logp(θ(i)|Y)=[Q(θ(i+1)|θ(i),Y)-Q(θ(i)|θ(i),Y)]-[H(θ(i+1)|θ(i),Y)-H(θ(i)|θ(i),Y)]由Jensen不等式,知Ez[log(θθ)|θ(i),Y]log{Ez(θθ)|θ(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.所以（6）成立。这说明用EM算法求出的θ(i)的确使其对数似然L(θ(i)|Y)关于EM算法及应用实例8i单增,但只有在严格单增的情况下,才能保证EM算法求出的θ(×)为θ的极大似然估计。这是从纯数学的,实际应用有一定的困难。定理2(1)如果p(θ|Y)有上界,则L(θ(i)|Y)收敛到某个L*。(2)如果Q(θ|φ)关于θ和φ都连续,则在关于L的很一般的条件下,由EM得到的估计序列θ(i)的收敛值θ*是L的稳定点。证明由定理1的证明及单调收敛定理易知（1）成立;(2)的证明见文献[4]。定理2表明EM算法的结果只能保证收敛到后验密度函数的稳定点,并不能保证收敛到极大值点,事实上,任何一种算法都很难保证其结果为极大值点。通常选取几个不同的初值进行迭代,然后在诸估计间加以选择,以减轻初值选取对结果的影响。2.2EM算法收敛速度的探讨EM算法是一种求参数极大似然估计的迭代算法,在处理不完全数据中有重要应用。EM算法实现简单,数值计算稳定,存储量小,并具有良好的全局收敛性。但是,EM算法收敛速度相当慢,只是次线性的收敛速度,这个缺点防碍了EM算法的应用。现已提出了多种加速EM算法收敛的方法,其中使用非线性规划中的Broyden对称秩1校正公式,提出了一种加速EM算法收敛的方法。与其他加速收敛方法相比,本方法简便易行,不必对不完全数据的似然函数一维搜索,在收敛性方面,与EM算法一样具有全局收敛性。数值试验结果表明,本方法的收敛速度比EM算法快的多。一般地,在EM算法的M步求θ(i+1)时,可求解方程组θθθ=0,得到参数的迭代公式θ(i+1)=G(θ(i))。这种迭代公式通常使EM算法的收敛速度很慢,为加速收敛,可考虑使用其他方法求解。根据著名的Fisher公式θθθθ=θ(i)=g(θ(i)),这里g(θ(i))是不完全数据对数似然函数L(θ)在θ(i)处的梯度g(θ(i))=▽L(θ)|θ=θ(i),于是问题转化为求θ(*),使g(θ(*))=0,从而可以使用非线性规划中的有效方法求解,达到加速收敛的目的。在非线性规划的求解方法中,Broyden对称秩1校正公式具有二次终止性EM算法及应用实例9和超线性的收敛速度,使用它可望改善EM步长的缓慢变化,以加速EM算法的收敛。但是,如果不对L(θ)进行一维搜索,Broyden对称秩1校正