531f10EM2 - x=c(58,12,9,13) n=sum(x) start=EM=cur=diff=0.1...

Info iconThis preview shows pages 1–4. Sign up to view the full content.

View Full Document Right Arrow Icon
Sheet1 Page 1 ############################################################# # STAT 531 EM2 ############################################################# da=rnorm(400)+2.5*sample(0:1, 400, rep=T, prob=c(1,3) ) like=function(mu){ -sum(log((.25*dnorm(da-mu[1])+.75*dnorm(da-mu[2])))) } x0 <-c(1,1) nlm(like, x0, hessian =T)$estimate x <- da EM=cur =c(-2,2) iter=20 for (i in 1: iter){ probs = (1/(1+3*dnorm(x-cur[2])/(dnorm(x-cur[1] )) )) mu1hat =sum(probs*x)/(sum(probs)) mu2hat=sum((1-probs)*x)/(sum(1-probs)) EM =rbind(EM, c(mu1hat, mu2hat)) cur=EM[dim(EM)[1],] } EM ############################################################# gen <- function (m1,m2,p,N) { w <- runif(N) < p x1 <- rexp(N,1/m1) x2 <- rexp(N,1/m2) w*x1 + (1-w)*x2 } xa<- gen(4, 20, 0.3, 1000) ############################################################# ### MCEM #############################################################
Background image of page 1

Info iconThis preview has intentionally blurred sections. Sign up to view the full version.

View Full DocumentRight Arrow Icon
Background image of page 2
Background image of page 3

Info iconThis preview has intentionally blurred sections. Sign up to view the full version.

View Full DocumentRight Arrow Icon
Background image of page 4
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: x=c(58,12,9,13) n=sum(x) start=EM=cur=diff=0.1 while (diff&gt;.001){ EM=c(EM,((cur*x[1]/(2+cur))+x[4])/((cur*x[1]/(2+cur))+x[2]+x[3]+x[4])) diff=abs(cur-EM[length(EM)]) cur=EM[length(EM)] } M=10^2 MCEM=start for (i in 2:length(EM)){ MCEM[i]=1/(1+(x[2]+x[3])/(x[4]+ mean( rbinom(100, x[1], prob=1/(1+2/MCEM[i-1])) ))) } Sheet1 Page 2 ############################################################# M=10^2 MCEM=matrix(start,ncol=length(EM),nrow=500) for (i in 2:length(EM)){ MCEM[,i]=1/(1+(x[2]+x[3])/(x[4]+rbinom(500,M*x[1], prob=1/(1+2/MCEM[,i-1]))/M)) } plot(EM,type=&quot;l&quot;,xlab=&quot;iterations&quot;,ylab=&quot;MCEM sequences&quot;) upp=apply(MCEM,2,max) polygon(c(1:length(EM),length(EM):1),c(upp,rev(dow)),col=&quot;grey78&quot;) lines(EM,col=&quot;gold&quot;,lty=2,lwd=2) ############################################################# Sheet1 Page 3 Sheet1 Page 4 dow=apply(MCEM,2,min)...
View Full Document

This note was uploaded on 02/04/2011 for the course STAT 531 taught by Professor Gaborlukacs during the Spring '11 term at Manitoba.

Page1 / 4

531f10EM2 - x=c(58,12,9,13) n=sum(x) start=EM=cur=diff=0.1...

This preview shows document pages 1 - 4. Sign up to view the full document.

View Full Document Right Arrow Icon
Ask a homework question - tutors are online