#source("c:/Work07/ShortCourse07/Programs/censoredGibbs.txt",print.eval=TRUE)#
#This does Gibbs 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)
m<-length(xdata)
n<-30;a<-4.5		#1/3 missing data
nsim<-10000
#----------------------Initial Values----------------------------
xbar<-mean(xdata)
that<-array(xbar,dim=c(nsim,1))
alpha<-a+1;M<-(1/alpha)*exp((alpha^2/2)-alpha*a)
z<-array(0,dim=c(n-m,1))
#----------------------Gibbs sampler-----------------------------
for(i in 1:nsim)
{
  for(j in 1:(n-m))
  {
    while(z[j]==0)		#Accept-Reject for truncated normal
    {
      zcand<-a+rexp(1,rate=alpha)
      test<-exp(zcand^2/2)/(M*exp(alpha*(zcand-a)))
      z[j]<-zcand*(runif(1)<test)
    }
  }
  zbar<-mean(z);z<-z*0
  that[i]<-rnorm(1,mean=(m/n)*xbar+(1-m/n)*zbar,sd=sqrt(1/n))
}
#------------------------Plot------------------------------------
par(mfrow=c(1,2))
plot(cumsum(that)/c(1:nsim),type="l",xlim=c(0,nsim),ylim=c(3.4,4),lwd=2)
hist(that,col="grey",breaks=25)