#source("c:/Work07/ShortCourse07/Programs/PoissonCompletionEM.txt",print.eval=TRUE)#
#This does the MCEM algorithm for the poisson completion#
nsim<-50;lam<-array(313/360,dim=c(nsim,1));ybar<-array(4,dim=c(nsim,1))
#Use m values for the mean
m<-25;y<-array(0,dim=c(m,1));
for (j in 2:nsim)  {
for(i in 1:m){while(y[i] < 4) y[i] <- rpois(1,lam[j -1])};
ybar[j]<-mean(y);
lam[j]<-(313+13*ybar[j])/360;
y<-y*0;
}
par(mfrow=c(1,2))
hist(ybar,col="green",breaks=10)
plot(lam,col="red",type="l",ylim=c(.95,1.05))

#To show standard deviation
v1<-var(-360+(313+13*ybar)/lam[nsim])
v2<-mean((313+13*ybar)/lam[nsim]^2)
sd<- sqrt(-1/(v1-v2))
top<-lam[nsim]+sd;bot<-lam[nsim]-sd
par(mfrow=c(1,1))
plot(lam,col="red",type="l",ylim=c(.95,1.10))
par(new=T)
lines(x=c(0,nsim),y=c(top,top),col="green")
par(new=T)
lines(x=c(0,nsim),y=c(bot,bot),col="green")