#This approach was suggested in Bayesian Data Analysis by Gelman et al.
#I like the simplicity of the update functions and the updating
#algorithm
mu.update=function(){
rnorm(1,(st*mu0+s0*sum(theta))/(st+K* s0), sqrt((st * s0)/(st + K * s0)))
}
st.update=function(){
(b1 + 0.5 * sum((theta - mu) * (theta - mu)))/ rgamma(1, a1 + K/2)
}
se.update=function(){
(b2 + 0.5 * sum((y - theta %o% rep(1,J))^2))/rgamma(1, a2 + (K * J)/2)
}
theta.update=function(){
theta1 <- rnorm(K, (se * mu)/(J * st + se), sqrt((st * se)/(J * st + se)))
theta1 + ((J * st)/(J * st + se)) * ybar
}
J=dim(y)[2]
K=dim(y)[1]
ybar=apply(y,1,mean)
#These should be reasonable start values.  If you choose larger values of
#(a1,b1) and (a2,b2) with the same expectation, your prior has smaller
#variance
mu0=mean(y)
s0=var(y)
a2=2
b2=1
a1=2
b1=170
n.chains=5
n.iter=100
sims=array(NA,c(n.iter,n.chains,K+3))
dimnames(sims)=list(NULL,NULL,c(paste("Theta[",1:6,"]",sep=""),"Mu","VarTheta","VarError"))
for(m in 1:n.chains){ 
  mu=rnorm(1,mean(y),sd(y))
  theta=rnorm(K,apply(y,1,mean),diag(apply(y,1,var)))
  for(t in 1:n.iter){
    st=st.update()    
    mu=mu.update()    
    se=se.update()
    theta=theta.update()
    sims[t,m,]=c(theta,mu,st,se)
   }
}
apply(sims,3,quantile,probs=c(.05,.50,.95))
#Change plot.parm (values from 1 to 9) to look at plots for other parameters
plot.parm=3
plot.ts(sims[,,plot.parm],main=paste(as.character(n.chains)," MCMC Chains for ",as.character(dimnames(sims)[[3]][plot.parm])),plot.type="single",
lty=rep(1,n.chains),col=seq(1,n.chains),xlab="Index",ylab=as.character(dimnames(sims)[[3]][plot.parm]))