#Gibbs sampler for hierarchical normal model

# Data 
y <- list()
y[[1]] = c(62,60,63,59)
y[[2]] = c(63,67,71,64,65,66)
y[[3]] = c(68,66,71,67,68,68)
y[[4]] = c(56,62,60,61,63,64,63,59)
y.mean = unlist(lapply(y,mean))
n.j = unlist(lapply(y,length))
n = sum(n.j)
J = length(y)

#Initial values
tau2 = 1
sigma2 = 1
mu = mean(y.mean)

S = 20000
theta.sim = matrix(nrow=S,ncol=J)
mu.sim = rep(NA,S)
sigma2.sim = rep(NA,S)
tau2.sim = rep(NA,S)
for(s in 1:S)
{
 #Simulate theta_j
 V.theta.inv = (n.j/sigma2)+(1/tau2)
 theta.hat = ((mu/tau2) + (n.j*y.mean/sigma2))/V.theta.inv
 theta = rnorm(J,theta.hat,sqrt(1/V.theta.inv))
 theta.sim[s,] = theta
 #show(theta)

 #Simulate mu
 mu.hat = mean(theta)
 mu = rnorm(1,mu.hat,sqrt(tau2/J))
 mu.sim[s] = mu
 #show(mu)
 
 #Simulate sigma2
 sigma2.hat = 0
 for(j in 1:J)
  sigma2.hat = sigma2.hat + sum((y[[j]]-theta[j])^2)
 sigma2.hat = sigma2.hat/n
 sigma2  = (n-1)*sigma2.hat/rchisq(1,n)
 sigma2.sim[s] = sigma2
 #show(sigma2)
 
 #Simulate tau2
 tau2.hat = sum((theta-mu)^2)/(J-1)
 tau2 = (J-1)*tau2.hat/rchisq(1,J-1)
 tau2.sim[s] = tau2
 #show(tau2)
}
par.sim = cbind(mu.sim,theta.sim,sqrt(sigma2.sim),sqrt(tau2.sim))
par.sum = t(apply(par.sim[-c(1:10000),],2,quantile,probs=c(0.025,0.25,0.5,0.75,0.975)))
par.sum = cbind(apply(par.sim,2,mean),apply(par.sim,2,sd),par.sum)
rownames(par.sum) = c("mu","theta_1","theta_2","theta_3","theta_4","sigma","tau")

colnames(par.sum)[1:2] = c("Mean","sd")
show(round(par.sum,3))
par(mfrow=c(3,2),mar=c(1,2,1,2))
plot.ts(par.sim[,1])
acf(par.sim[,1])
plot.ts(par.sim[,6])
acf(par.sim[,6])
plot.ts(par.sim[,7])
acf(par.sim[,7])