#source("c:/Work07/ShortCourse07/Programs/censoredEM-MCEM.txt",print.eval=TRUE)#
#This does the em and mcem  algorithm for censored data#
#n obs normal(t0,1), m uncensored values, n-m censored at a#
xdata<-c(3.64, 2.78, 2.91,2.85,2.54,2.62,3.16,2.21,4.05,2.19,2.97,4.32,
3.56,3.39,3.59,4.13,4.21,1.68,3.88,4.33)
n<-25;m<-20;t0<-4;a<-4.5;nt<-50
xbar<-mean(xdata);that<-array(xbar,dim=c(nt,1));
for (j in 2:nt)  {
that[j] <-(m/n)*xbar+(1-m/n)*(that[j-1]+dnorm(a-that[j-1])/(1-pnorm(a-that[j-1])))}
#now do MCEM, z=missing data, nz=size of MC sample
tmc<-array(xbar,dim=c(nt,1));nz<-25;
for (j in 2:nt)  {
z<-array(a-1,dim=c(nz,1));
for (k in 1:nz)  {while(z[k] <a) z[k] <- rnorm(1,mean=tmc[j-1],sd=1)}
zbar<-mean(z)
tmc[j] <-(m/n)*xbar+(1-m/n)*zbar}
plot(that,type="l",xlim=c(0,nt),ylim=c(3.3,3.7),lwd=2)
par(new=T)
plot(tmc,type="l",xlim=c(0,nt),ylim=c(3.3,3.7),xlab="",ylab="",xaxt="n",yaxt="n",lwd=2)