p+)= 2m鼎+n鼎+n8 2n 2m唱+n愠+n温 (t) (t) +1)= PT 2n 因此，使用以上结论便可进行EM算法： 千Code moth<-function(p,n.obs){ n<-sum(n.obs) nc<-n.obs[1] ni<-n.obs[2] nt<-n.obs [3] ntt<-nt cat(p,"\n") pct<-pit<-ptt<-rep(0,20) pct[1]<-p[1] pit[1]<-p[2] ptt[1]<-1-p[1]-p[2] for(i in 2:20){ Previous Next First Last Back Forward 9
p (t+1) I = 2n (t) II + n (t) IT + n (t) CI 2n , p (t+1) T = 2n (t) T T + n (t) CT + n (t) IT 2n . œd, ¶^±˛(ÿBå?1EMé{: ↑Code moth<-function(p,n.obs){ n<-sum(n.obs) nc<-n.obs[1] ni<-n.obs[2] nt<-n.obs[3] ntt<-nt cat(p,"\n") pct<-pit<-ptt<-rep(0,20) pct[1]<-p[1] pit[1]<-p[2] ptt[1]<-1-p[1]-p[2] for(i in 2:20){ Previous Next First Last Back Forward 9
pc.old<-pct [i-1] pi.old<-pit [i-1] pt.old<-ptt [i-1] den<-pc.old~2+2*pc.old*pi.old+2*pc.old*pt.old ncc<-nc*pc.old~2/den nci<-2*nc*pc.old*pi.old/den nct<-2*nc*pc.old*pt.old/den nii<-ni*pi.old2/(pi.old-2+2*pi.old*pt.old) nit<-2*ni*pi.old*pt.old/(pi.old-2+2*pi.old*pt.old) pct [i]<-(2*ncc+nci+nct)/(2*n) pit[i]<-(2*nii+nit+nci)/(2*n) ptt[i]<-(2*ntt+nct+nit)/(2*n) return(list(pct=pct,pit=pit,ptt=ptt)) n.obs<-c(85,196,341)#observed data,n_c,n_I,n_T p<-c(1/3,1/3) a<-moth(p,n.obs) Previous Next First Last Back Forward 10
pc.old<-pct[i-1] pi.old<-pit[i-1] pt.old<-ptt[i-1] den<-pc.old^2+2*pc.old*pi.old+2*pc.old*pt.old ncc<-nc*pc.old^2/den nci<-2*nc*pc.old*pi.old/den nct<-2*nc*pc.old*pt.old/den nii<-ni*pi.old^2/(pi.old^2+2*pi.old*pt.old) nit<-2*ni*pi.old*pt.old/(pi.old^2+2*pi.old*pt.old) pct[i]<-(2*ncc+nci+nct)/(2*n) pit[i]<-(2*nii+nit+nci)/(2*n) ptt[i]<-(2*ntt+nct+nit)/(2*n) } return(list(pct=pct,pit=pit,ptt=ptt)) } n.obs<-c(85,196,341) # observed data,n_c,n_I,n_T p<-c(1/3,1/3) a<-moth(p,n.obs) Previous Next First Last Back Forward 10
pct<-aSpct pit<-aSpit ptt<-aSptt #convergence diagnostics statistic R rcc=sqrt((diff(pct)-2+diff(pit)2)/(pct [-20]-2+pit[-20]-2)) rcc=c(0,rcc)#adjusts the length to make the table below d1=(pct[-1]-pct[20])/(pct[-20]-pct[20]) d1=c(d1,0) d2=(pit[-1]-pit[20])/(pit[-20]-pit[20]) d2=c(d2,0) #Table output print(cbind(pct,pit,rcc,d1,d2)[1:9,],digits=5) Code 其中，收敛的标准为)=p山为由一次迭代到下一次送代 p-1)‖ 在p:-1)上相对改变的总量.（本例中我们以迭代10次为例，没有使用此标准控 制收敛).最后两列是验证EM算法的收敛速度为线性的. Previous Next First Last Back Forward 11
pct<-a$pct pit<-a$pit ptt<-a$ptt #convergence diagnostics # statistic R rcc=sqrt( (diff(pct)^2+diff(pit)^2)/(pct[-20]^2+pit[-20]^2) ) rcc=c(0,rcc) #adjusts the length to make the table below d1=(pct[-1]-pct[20])/(pct[-20]-pct[20]) d1=c(d1,0) d2=(pit[-1]-pit[20])/(pit[-20]-pit[20]) d2=c(d2,0) #Table output print(cbind(pct,pit,rcc,d1,d2)[1:9,],digits=5) ↓Code Ÿ•,¬ÒIOèR(t) = kp (t)−p (t−1)k kp(t−1)k èdògSì eògSì 3p (t−1)˛ÉÈUCo˛. ( ~•·Ç±Sì10gè~,vk¶^dIOõ õ¬Ò). Ÿ¥yEMé{¬ÒÑ›èÇ5. Previous Next First Last Back Forward 11
例2 Bayes后验众数考虑一个具有似然L(y),先验π(8)以及缺失数据 或者参数Z(即y=(x,z)的Bayes问题.为找到后验众数，E步需要 Q(00())E{log(L(Oly)(0)k(y)),0()} EflogL(0ly)l,0()}+logn(0)+Eflogk(y)r,0()} 其中最后一项可以在最大化Q中略去，因其与无关.显然，对此Bayes.后验众 数问题，和经典统计方法下的差异在于多了一项先验的对数.因此此时EM算 法为 1.E步：计算Q(00(t)=E{logL(y)z,a(t)} 2.M步：计算0(t+1)=argmazQ(0lee) 考虑X=(X1,X2,X3)~MN(n,(2+θ)/4，(1-)/2,0/4)，为应用EM算 法估计8，我们视完全数据为Y=(Z11,Z12,X2,X3)~MN(n,1/2,0/4,(1- 0)/2,8/4),其中Z11+Z12=X1 Previous Next First Last Back Forward 12
~2 Bayes ØÍ ƒòá‰kq,L(θ|y), kπ(θ)±9"îÍ‚ ½ˆÎÍZ(=y = (x, z)) BayesØK. èÈØÍ, E⁄Iá Q(θ|θ (t) ) = E{log(L(θ|y)π(θ)k(y))|x, θ(t) } = E{logL(θ|y)|x, θ(t) } + logπ(θ) + E{logk(y)|x, θ(t) } Ÿ•Åòëå±3ÅåzQ•—, œŸÜθÃ'. w,, ÈdBayes Ø ÍØK, ⁄²;⁄Oê{e …3uı òëkÈÍ. œddûEMé {è 1. E⁄: OéQ(θ|θ (t) ) = E{logL(θ|y)|x, θ(t)}, 2. M⁄: Oéθ (t+1) = argmaxQ(θ|θ (t) ) ƒX = (X1, X2, X3) ∼ MN(n, (2+θ)/4, (1−θ)/2, θ/4), èA^EMé {Oθ, ·Ç¿ Í‚èY = (Z11, Z12, X2, X3) ∼ MN(n, 1/2, θ/4, (1− θ)/2, θ/4), Ÿ•Z11 + Z12 = X1. Previous Next First Last Back Forward 12
因此，有 I(0Y)=(Z12+X3)loge X2log(1-0)+constants E[l(Y)1X,6]=(E[Z12lZ11+Z12=X1,]+X3)log9 +X2log(1-0)+constants X19e) =(+X)og+Xalog(1-0)+constants 考虑0的先验为Beta(a,b), x(0)= r(a+b）ga-11-b-1. (a)+T(b) 因此 Q(0j0())= X10() (+a-1)og0+(2+-1)og(1-0)+constants 所以得到 8(t+1)= X10() 2+9+X3+a- ）)/(+++a+-2） Previous Next First Last Back Forward 13
œd, k l(θ|Y ) = (Z12 + X3)logθ + X2log(1 − θ) + constants E[l(θ|Y )|X, θ(t) ] = (E[Z12|Z11 + Z12 = X1, θ(t) ] + X3)logθ +X2log(1 − θ) + constants = ( X1θ (t) 2 + θ (t) + X3)logθ + X2log(1 − θ) + constants ƒθkèBeta(a, b), π(θ) = Γ(a + b) Γ(a) + Γ(b) θ a−1 (1 − θ) b−1 . œd Q(θ|θ (t) ) = ( X1θ (t) 2 + θ (t) + X3 + a − 1)logθ + (X2 + b − 1)log(1 − θ) + constants §± θ (t+1) = X1θ (t) 2 + θ (t) + X3 + a − 1 ! . X1θ (t) 2 + θ (t) + X3 + X2 + a + b − 2 ! Previous Next First Last Back Forward 13