中国科学技术大学：《实用统计软件》课程课件讲义（统计计算与软件）第十讲 Expectation-Maximization（EM算法）方法

Lecture 10:Expectation-Maximization (EM) 张伟平 Monday 16th November,2009

Lecture 10: Expectation-Maximization (EM) ê{ ‹ï² Monday 16th November, 2009

Contents 1 EM optimization method 1.1 EM algorithm........···.·..···.· 2 l.2 Convergence..:·。。··.·····。······ 15 l.3 Usage in exponential families·.:.......·.· 19 l.4 Usage in finite normal mixtures...·..····.···· 20 1.5 Variance estimation 23 1.5.1 Louis method..... 24 1.5.2 SEM algorithm . 28 1.5.3 Bootstrap method.. 36 1.5.4 Empirical Information 37 1.6 EM Variants . 444444 38 l.6.1 Improving the E step:...·。.。·。.··· 38 1.6.2 Improving the M step 4+。·4。·。4.4。·。。。· 39 1.7 Pros and Cons.····.············· 。4 40 Previous Next First Last Back Forward 1

Contents 1 EM optimization method 1 1.1 EM algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 1.3 Usage in exponential families . . . . . . . . . . . . . . . . . 19 1.4 Usage in finite normal mixtures . . . . . . . . . . . . . . . . 20 1.5 Variance estimation . . . . . . . . . . . . . . . . . . . . . . 23 1.5.1 Louis method . . . . . . . . . . . . . . . . . . . . . . 24 1.5.2 SEM algorithm . . . . . . . . . . . . . . . . . . . . . 28 1.5.3 Bootstrap method . . . . . . . . . . . . . . . . . . . 36 1.5.4 Empirical Information . . . . . . . . . . . . . . . . . 37 1.6 EM Variants . . . . . . . . . . . . . . . . . . . . . . . . . . 38 1.6.1 Improving the E step . . . . . . . . . . . . . . . . . . 38 1.6.2 Improving the M step . . . . . . . . . . . . . . . . . 39 1.7 Pros and Cons . . . . . . . . . . . . . . . . . . . . . . . . . 40 Previous Next First Last Back Forward 1

Chapter 1 EM optimization method 期望-最大化(EM)算法是一种在观测到数据后，估计未知参数的迭代优化方 法.其能非常简单地执行并且能够通过稳定，上升的步骤非常可靠地找到全局 最优值.对EM方法详细介绍请参考非常好的教材McLachlan and Krishnan, 1997,EM algorithm and extensions,Wiley. 在频率概率框架下，我们可以认为：除了观测到的样本X外，还有一些与 之伴随的未知的缺失量(missing)或者没有观察到的量Z.那么，完整的样本 应该是Y=(X,Z).问题的目的是在得到观测的样本x后，最大化似然函 数L(x)=∫f(x,z旧)dz.此似然函数由于涉及到边际概率函数一般难以处 理，而采用f(x,旧)和f(zx,)则可能容易处理，EM算法就是通过采用这些较 容易的密度而避免直接考虑L(x): Previous Next First Last Back Forward

Chapter 1 EM optimization method œ"-Ååz(EM)é{¥ò´3*ˇÍ‚￾, OôÎÍSì`zê {. Ÿ Uö~{¸/â1øÖU œL­½,˛,⁄½ö~åÇ/Ȥ Å`ä. ÈEMê{ç[0ûÎö~–· McLachlan and Krishnan, 1997, EM algorithm and extensions, Wiley. 3™«V«µee, ·Çå±@è:ÿ *ˇX , Ñkò Ü Éäëô "î˛(missing)½ˆvk* ˛Z. @o,  AT¥Y = (X, Z). ØK8 ¥3*ˇx￾, Ååzq,º ÍL(θ|x) = R f(x, z|θ)dz. dq,ºÍ du9>SV«ºÍòÑJ±? n, Ê^f(x, z|θ)⁄f(z|x, θ)KåUN¥?n, EMé{“¥œLÊ^˘  N¥ó› ;ù܃L(θ|x). Previous Next First Last Back Forward 1

另外，在Byes理论框架下，感兴趣的通常是对后验概率函数f(0x)最大化， 以后验众数（似然思想）来估计参数.此时，最优化问题有时可以通过考虑引入 一个未观测的参数而得到简化 有时候，缺失数据Z并非真正的缺失了，它们可能仅仅是为了简化我们的问 题而引入的变量.此时，Z常常被称为隐变量(latent variable) 1.1 EM algorithm 记可观测的量为X,缺失量为Z,完全数据为Y=(X,Z),待估参数为8.则 我们有的信息是观测数据的似然L(x),最大化此似然或者说是求的极大 似然估计是我们的目标.EM算法通过迭代方式来寻求最大化L(x)的解.假 设0（表示在第t次迭代后的最大值点，t=0,1,2,··.定义Q(8)为在观测 到X=x,以及在参数0=)的条件下完全数据的联合对数似然函数的期望： 此期望是对f21x(zx,8()计算.即 Q(0l0())=EflogL(0Y)I,0()} Previous Next First Last Back Forward 2

, , 3Bayesnÿµee, a,œ~¥È￾V«ºÍf(θ|x)Ååz, ±￾ØÍ(q,gé)5OÎÍ. dû, Å`zØKkû屜Lƒ⁄\ òáô*ˇÎÍψ {z. kûˇ, "îÍ‚Zøö˝"î , ßÇåU==¥è {z·ÇØ K ⁄\C˛. dû, Z~~°è ¤C˛(latent variable). 1.1 EM algorithm På*ˇ˛èX, "î˛èZ, Í‚èY = (X, Z), ñÎÍèθ. K ·Çk&E¥*ˇÍ‚ q,L(θ|x), Ååzdq,½ˆ`¥¶θ4å q,O¥·Ç8I. EMé{œLSìê™5œ¶ÅåzL(θ|x)). b θ (t)L´31tgSì￾Ååä:, t = 0, 1, 2, · · · . ½¬Q(θ|θ (t) )è3*ˇ X = x,±93ÎÍθ = θ (t)^áe Í‚È‹ÈÍq,ºÍœ", dœ"¥ÈfZ|X(z|x, θ(t) )Oé. = Q(θ|θ (t) ) = E{logL(θ|Y )|x, θ(t) } Previous Next First Last Back Forward 2

Eflogf(,Z10)r,0() logf(z,210)f(0())dz 最后一式强调一旦我们给定了X=x,Z就是Y唯一的随机部分 EM算法从8(O)开始，然后在两步之间交替：E表示期望，M表示最大化.该 算法概括如下： (1)E步：计算Q(ee)=EL(0Y)z.9)1 (2)M步：关于0最大化Q((t).并记9(+1)表示此时的最大值点 (3)返回到E步，直至收敛准则达到 例1椒花蛾(peppered moth)一个进化和工业污染的例子. 白桦尺蛾(peppered moth),麟翅目(Lepidoptera)尺镀科(Geometridae)昆 虫，学名Biston betularia.翅黑或白色，上有斑点，分布欧洲.其生命周 期有四个阶段：卵，毛毛虫.前.成虫.在由蛹变为成虫的时刻，它们会停留 在树干上完成此过程，此时就会被鸟类捕食.1848年，在英国曼彻斯特首 次注意到黑色型，到1898年黑色型反而以99：1的比例超过了淡色型.对这 种现象的解释是：工业区的树干被煤烟染黑，黑色的蛾栖于树上，目标不显 著，不易为鸟类捕食.这是通过工业黑化现象进行自然选择的一个例子. 一百度百科 Previous Next First Last Back Forward 3

= E{logf(x, Z|θ)|x, θ(t) } = Z logf(x, z|θ)f(z|x, θ(t) )dz Å￾ò™rNò￾·Çâ½ X = x, Z“¥Y çòëÅ‹©. EMé{lθ (0)m©, ,￾3¸⁄ÉmO: EL´œ", ML´Ååz. T é{V)Xe: (1) E⁄: OéQ(θ|θ (t) ) = E[L(θ|Y )|x, θ(t) ]. (2) M⁄: 'uθÅåzQ(θ|θ (t) ), øPθ (t+1)L´dûÅåä:. (3) à£E⁄, Üñ¬ÒOKà. ~1 ˛sˇ(peppered moth) òá?z⁄Ûí¿/~f. xºˇ(peppered moth),º8(Lepidoptera)ºˇâ(Geometridae)& £, ƶBiston betularia. ºÁ½x⁄, ˛kÄ:, ©ŸÓ³. Ÿ)·± œkoá„: Œ, ff£, W, §£. 3dWC解ûè, ßÇ ¨ 3 3‰Z˛§dLß, dû“¨ja”†. 1848c, 3=I˘îdA ƒ g5øÁ⁄., 1898cÁ⁄.á ±99:1'~áL ⁄.. È˘ ´yñ)º¥: Ûí´‰ZuÎ/Á, Á⁄ˇ—u‰˛, 8Iÿw Õ, ÿ¥èja”†. ˘¥œLÛíÁzyñ?1g,¿Jòá~f. —–z›zâ Previous Next First Last Back Forward 3

这种娥子的色彩己经被确认是通过某单个基因决定的，该基因有三个可能 的等位基因：C,I和T.其中C对1是显性的，T对I是隐性的.因此基因型CC CI和CT导致（深色）黑色的表型.基因型TT导致浅色的表型.Ⅱ和T导致一种 中间的表型，外观上变化很广泛，但通常以中间色彩杂色而成.因此有6种基因 型，3种表现型 在英因和北美，在燃煤的工业区，浅色的蛾子几乎以及被深色的蛾子所替 代.这种等位基因频率在种群内的变化被引为是在人类时间刻度下可以观察 到微进化的一个例子.（被试验证实的）理论结果是”鸟类在不同反射背景下对 蛾子的捕食程度明显不同”.在燃煤的工业区，污染减弱了黑色表型蛾子栖息 Previous Next First Last Back Forward 4

˘´ˇf⁄ÁƲ(@¥œL,¸áƒœ˚½, TƒœknáåU †ƒœ: C, I⁄T. Ÿ• CÈI¥w5, TÈI¥¤5. œdƒœ.CC, CI⁄CTó(⁄)Á⁄L.. ƒœ.TTóf⁄ L.. II⁄ITóò´ •mL., *˛CzÈ2ç, œ~±•m⁄Á,⁄ §. œdk6´ ƒœ ., 3´Ly.. 3=I⁄{, 3-uÛí´, f⁄ˇfA±9⁄ˇf§O ì. ˘´†ƒœ™«3´+ SCz⁄è¥3<aûmè›eå±* á?zòá~f. (£y¢)nÿ(J¥”ja 3ÿ”áµeÈ ˇf”†ß›²wÿ””. 3-uÛí´, ¿/~f Á⁄L.ˇf—E Previous Next First Last Back Forward 4

在树皮表面的反射程度，因而对其有利.当环境改善后，浅色表型增加，黑色表 型减少，这并不令人奇怪 因此，监视这种蛾子的等位基因C,I,T的变化，除了可以研究微进化过 程外，也可以用于环境污染监控.假设Hardy-Veinberg平衡律成立，记 各个等位基因在种群中的频率为Pc,P1,Pr(pC+PI十Pr=1),则基因 型CC,CI,CT,IL,IT,TT的频率分别为P吃，2pcpI,2pcpT,p听，2ppT,P咛 假设我们随机捕了m只蛾子，其中黑色表型的有nc只，中间表型的有n1只， 浅色表型的有T只.但是各个基因型的频数不可观测： CC CI CT Π IT TT 不可观测 ncc nCI nCT nIl nIT nTT nc n nT 观测 黑色 中间 浅色 因此，此时观测数据为x=(nc,nI,nr),而完全数据为y=(ncc,ncI,ncT, nII,nIT,nTT),这里ncc+ncI+ncT=nc,nl+nIT=nI,以及 nrT=nr.我们希望由此数据来估计各个等位基因的概率pc,p1,pT. Previous Next First Last Back Forward 5

3‰ôL°áß›, œ ÈŸk|. ǸUı￾, f⁄L.O\, Á⁄L .~, ˘øÿ-<¤%. œd, i¿˘´ˇf†ƒœC,I,TCz, ÿ å±Ôƒá?zL ß , èå±^u Ǹ¿/iõ. bHardy-Weinberg²ÔƧ·, P àᆃœ3´+•™«èpC , pI , pT (pC + pI + pT = 1), Kƒœ .CC,CI,CT,II,IT,TT ™«©Oèp 2 C , 2pC pI , 2pC pT , p2 I , 2pI pT , p2 T . b·ÇëÅ” nêˇf, Ÿ•Á⁄L.knCê, •mL.knIê, f⁄L.knT ê. ¥àნ.™Íÿå*ˇ: CC CI CT II IT TT ÿå*ˇ nCC nCI nCT nII nIT nT T nC nI nT *ˇ Á⁄ •m f⁄ œd, dû*ˇÍ‚èx = (nC , nI , nT ), Í‚èy = (nCC , nCI , nCT , nII , nIT , nT T ), ˘pnCC + nCI + nCT = nC , nII + nIT = nI , ±9 nT T = nT . ·ÇF"ddÍ‚ 5OàᆃœV«pC , pI , pT . Previous Next First Last Back Forward 5

从而完全数据的似然函数为 logf(yp) ncclogpe +ncilog(2PcPI)+ncrlog(2PCPT) nirlog(pi)+nrrlog(2pIpr)+nrrlog(pir) n! log nccmncrncn 由于完全数据不可观测，若记NcC,Nc,NcT,N,Nr,Nrr分别为各 个基因型的个数，p=(PCC,pPCI,PCT,pII,pIT,PTT,则有 1.NCC:NCI,NcTlnc,nI,nT,p~ MN (pc)2 2PCPI 2PCPT nc:1-(pI+Pr)1-(PI+PT)2'1-(PI+PT)2 2.NII,NITlnc,nI,nT:P~ (pr)2 2PCPI MN I2I(2PPI+() 3.NTrlnc:nI,nT:P~B(nr:pi) Previous Next First Last Back Forward 6

l Í‚q,ºÍè logf(y|p) = nCC logp2 C + nCI log(2pC pI ) + nCT log(2pC pT ) + nII log(p 2 I ) + nIT log(2pI pT ) + nT T log(p 2 T T ) + log n! nCC !nCI !nCT !nII !nIT !nT T ! duÍ‚ÿå*ˇ, ePNCC , NCI , NCT , NII , NIT , NT T ©Oèà ნ .áÍ, p = (pCC , pCI , pCT , pII , pIT , pT T ), Kk 1. NCC , NCI , NCT |nC , nI , nT , p ∼ MN  nC , (pC ) 2 1 − (pI + pT ) 2 , 2pC pI 1 − (pI + pT ) 2 , 2pC pT 1 − (pI + pT ) 2  2. NII , NIT |nC , nI , nT , p ∼ MN  nI , (pI ) 2 2pC pI + (pT ) 2 , 2pC pI 2pC pI + (pT ) 2  3. NT T |nC , nI , nT , p ∼ B(nT , p2 T ) Previous Next First Last Back Forward 6

从而在EM算法的第步中，对完全似然函数取期望得到 Q(plp())=ncclog(pe)+nlog(log(2pcPr)+nerlog(2pcpr) +nlog(pi)+nlog(2pipr)+nerlog(pir) k(nc,nI,nr,p(t)), 其中 nce E(Ncclnc,ni,nr.p(}=nc- 892 1-(p9+p92 ne E(Ncclnc,ni,nT,p()=nc. 2p99 1-p9+p92 ner E(Ncclnc,ni,nr,p())=nc- 2p99 1-p9+p92 p92 n E(Nuncr Previous Next First Last Back Forward 7

l 3EMé{1t⁄•, Èq,ºÍœ" Q(p|p (t) ) = n (t) CC log(p 2 C ) + n (t) CI log(log(2pC pI ) + n (t) CT log(2pC pT ) + n (t) II log(p 2 I ) + n (t) IT log(2pI pT ) + n (t) T T log(p 2 T T ) + k(nC , nI , nT , p(t) ), Ÿ• n (t) CC = E{NCC |nC , nI , nT , p(t) } = nC (p (t) C ) 2 1 − (p (t) I + p (t) T ) 2 , n (t) CI = E{NCC |nC , nI , nT , p(t) } = nC 2p (t) C p (t) I 1 − (p (t) I + p (t) T ) 2 , n (t) CT = E{NCC |nC , nI , nT , p(t) } = nC 2p (t) C p (t) T 1 − (p (t) I + p (t) T ) 2 , n (t) II = E{NII |nC , nI , nT , p(t) } = nI (p (t) I ) 2 2p (t) C p (t) I + (p (t) T ) 2 Previous Next First Last Back Forward 7

E{NITlnc,nI:nr,p())=nI 2p9 2p9p9+(p92 最后一项 k(nc,nI:nr:p()) n! Eflog- cIncrncrIniinITInTr! c,nI,nr,p(t)y 与p无关. 下面对Q(plpO)进行最大化，注意pc+pr十pT=1,于是关于pC,p1求导 aQ(plp() 2m思+n8+n_2n鼎+n+n唱 opc PC 1-PC-PI aQ(plp() 2m9+n鼎+n8 () 2n鼎+n唱+n出 OpI PI 1-PC -PI 令导数为零，得到 p+= 2n+n思+n温 2n Previous Next First Last Back Forward 8

n (t) IT = E{NIT |nC , nI , nT , p(t) } = nI 2p (t) I p (t) T 2p (t) C p (t) I + (p (t) T ) 2 . Å￾òë k(nC , nI , nT , p(t) ) = E{log n! nCC !nCI !nCT !nII !nIT !nT T ! |nC , nI , nT , p(t) } ÜpÃ'. e°ÈQ(p|p (t) )?1Ååz, 5øpC + pI + pT = 1, u¥'upC , pI¶ ∂Q(p|p (t) ) ∂pC = 2n (t) CC + n (t) CI + n (t) CT pC − 2n (t) T T + n (t) CT + n (t) IT 1 − pC − pI , ∂Q(p|p (t) ) ∂pI = 2n (t) II + n (t) II + n (t) CI pI − 2n (t) T T + n (t) CT + n (t) IT 1 − pC − pI -Íè", p (t+1) C = 2n (t) CC + n (t) CT + n (t) CI 2n , Previous Next First Last Back Forward 8