* * EMEXAMPLE.RPF * Example of the EM optimization algorithm * From User's Guide, Section 4.14 * all 6237 open data g10xrate.xls data(format=xls,org=columns) / usxjpn usxfra usxsui * * This estimates a "mixture" model, with a zero mean and two branches * for the variance. * set xjpn = 100.0*log(usxjpn/usxjpn{1}) compute rstart=2,rend=6237 * * Get guess values based on the sample variance. * stats xjpn compute lowv=.5*%variance,highv=2.0*%variance ****************************************************************** * * Direct maximization * nonlin p sigsq1 sigsq2 frml logl = log(p*exp(%logdensity(sigsq1,xjpn))+\$ (1-p)*exp(%logdensity(sigsq2,xjpn))) * compute t1=%cputime() compute p=.5 compute sigsq1=lowv,sigsq2=highv maximize logl rstart rend ****************************************************************** * * Classic "EM" calculation * compute t2=%cputime() compute p=.5 compute sigsq1=lowv,sigsq2=highv compute oldparms=||p,sigsq1,sigsq2|| do iters=1,500 set pstar = f1=exp(%logdensity(sigsq1,xjpn)),\$ f2=exp(%logdensity(sigsq2,xjpn)),\$ p*f1/(p*f1+(1-p)*f2) sstats(mean) rstart rend pstar>>p pstar*xjpn^2>>sumsq1 \$ (1-pstar)*xjpn^2>>sumsq2 compute sigsq1=sumsq1/p,sigsq2=sumsq2/(1-p) compute newparms=||p,sigsq1,sigsq2|| if %testdiff(newparms,oldparms)<1.e-6 break compute oldparms=newparms end disp "Estimates from classical EM" p sigsq1 sigsq2 disp "Iterations" iters ****************************************************************** * * EM using the ESTEP option * function EStep set pstar = f1=exp(%logdensity(sigsq1,xjpn)),\$ f2=exp(%logdensity(sigsq2,xjpn)),\$ p*f1/(p*f1+(1-p)*f2) end EStep * compute t3=%cputime() frml emlogl = pstar*%logdensity(sigsq1,xjpn)+\$ (1-pstar)*%logdensity(sigsq2,xjpn)+\$ pstar*log(p)+(1-pstar)*log(1-p) * compute p=.5 compute sigsq1=lowv,sigsq2=highv nonlin p sigsq1 sigsq2 maximize(estep=estep(),method=bhhh,cvcrit=1.e-6,iters=500) emlogl rstart rend compute t4=%cputime() * disp "Timing" disp "ML" @15 t2-t1 disp "Classical EM" @15 t3-t2 disp "Max(ESTEP)" @15 t4-t3