#`
#This does the gibbs sampler for the poisson completion#
#Adds Rao-Blackwellized Estimator
#Eliminates WHILE step - direct simulation of truncated Poisson
nsim<-5000
lam<-array(313/360,dim=c(nsim,1))
z<-array(0,dim=c(13,1));
RB<-array(313/360,dim=c(nsim,1))
#--------------------------------------------------------
for (j in 2:nsim)  
{
  top<-round(lam[j -1]+6*sqrt(lam[j -1]))
  prob<-dpois(c(4:top),lam[j -1])
  cprob<-cumsum(prob/sum(prob))
  for(i in 1:13){z[i] <- 4+sum(cprob<runif(1))}
  a<-313+sum(z)+1
  lam[j]<-rgamma(1,a,scale=1/360);
  RB[j]<-(313+sum(z))/360
  z<-z*0;
}
#----------------------Plot-------------------------------
den<-1:(nsim)
meanlam<-cumsum(lam)/den;
par(mfrow=c(1,2))
plot(meanlam,type="l",ylim=c(.9,1.1),xlab="iteration",ylab="estimate",col="red")
par(new=T)
plot(cumsum(RB)/den,type="l",ylim=c(.9,1.1),xlab="iteration",ylab="estimate",col="blue")
hist(lam,main="Mean",freq=F,col="green")