silly.model<-function(y,a=0.01,b=0.01,c=100,n.samples=50000){

   #This function returns n.samples samples from the model:
   #    y|alpha,beta,tau~N(alpha+beta,tau)
   #    alpha,beta~N(0,c^2)
   #    tau~Gamma(a,b)
   #
   # y is the vector of outcomes

   n<-length(y)

   #intial values:
   tau<-rgamma(1,1,1)
   alpha<-rnorm(1,0,100)
   beta<-rnorm(1,0,100)

   #Initialize vectors to store the results:
   keep.tau<-keep.beta<-keep.alpha<-rep(0,n.samples)

   #Start the MCMC sampler!
   for(i in 1:n.samples){
     #update tau:
     tau<-rgamma(1,n/2+a,rate=sum((y-alpha-beta)^2)/2+b)

     #update alpha:
     VVV<-n*tau+1/c^2
     MMM<-tau*sum(y-beta)
     alpha<-rnorm(1,MMM/VVV,1/sqrt(VVV))

     #update beta:
     VVV<-n*tau+1/c^2
     MMM<-tau*sum(y-alpha)
     beta<-rnorm(1,MMM/VVV,1/sqrt(VVV))

     #Use for block updates
     #XXX<-matrix(1,n,2)
     #VVV<-solve(tau*t(XXX)%*%XXX + diag(2)/(c^2))
     #MMM<-tau*t(XXX)%*%y
     #BBB<-VVV%*%MMM + t(chol(VVV))%*%rnorm(2)
     #alpha<-BBB[1]
     #beta<-BBB[2]


     #store results:
     keep.tau[i]<-tau
     keep.alpha[i]<-alpha
     keep.beta[i]<-beta

     #Display the chains every 10000 iterations:
     if(i%%10000==0){
       par(mfrow=c(3,1))
       plot(keep.tau[1:i],type="l",ylab=expression(tau),xlab="Iteration")
       plot(keep.alpha[1:i],type="l",ylab=expression(alpha),xlab="Iteration")
       plot(keep.beta[1:i],type="l",ylab=expression(beta),xlab="Iteration")
     }
   }

#return a list with the posterior samples:
list(tau=keep.tau,alpha=keep.alpha,beta=keep.beta)}

#Generate data:
n<-100                     #samples size of the simulated data
n.samples<-50000           #number of MCMC samples to generate

set.seed(0820)
y<-rnorm(n)

#fit the model
chain1<-silly.model(y,c=100,n.samples=50000)
chain2<-silly.model(y,c=100,n.samples=50000)


#Convergence plots
par(mfrow=c(2,2))
beta<-cbind(chain1$beta,chain2$beta)
plot(beta[,1],type="l",ylim=range(beta))
lines(beta[,2],col=2)
legend("topright",c("Chain 1","Chain 2"),inset=0.05,lty=1,col=1:2)

acf.beta<-acf(beta[,1],lag.max=10000)

alpha<-cbind(chain1$alpha,chain2$alpha)
plot(alpha[,1],type="l",ylim=range(alpha))
lines(alpha[,2],col=2)
legend("topright",c("Chain 1","Chain 2"),inset=0.05,lty=1,col=1:2)

plot(alpha,beta)


#Summarize the posterior of beta
beta.hat<-mean(beta[,1])
beta.N<-length(beta[,1])
beta.ess<-beta.N/(1+2*sum(acf.beta$acf[-1]))
beta.hat.se<-sd(beta[,1])/sqrt(beta.ess)

print("Summary of beta's posterior")
print("Posterior mean")
print(beta.hat)
print("MC standard error")
print(beta.hat.se)
print("Effective sample size")
print(beta.ess)


#Compute Gelman/Rubin statistics
library(Bolstad2)
GelmanRubin(alpha)