#source("c:/Work07/ShortCourse07/Programs/PoissonCompletion2.txt",print.eval=TRUE)#
#This does the gibbs sampler for the poisson completion#
#Adds Rao-Blackwellized Estimator
nsim<-500
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)  {
for(i in 1:13)
{
	while(z[i] < 4) z[i] <- rpois(1,lam[j -1])}
	a<-313+sum(z)+1
	lam[j]<-rgamma(1,a,scale=1/360);
	RB[j]<-(313+sum(z))/360
	z<-z*0;
}
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")