#
# Hierachical normal model at end of Section 9.2
#
niter <- 25
r <-  4
n <- c(4,6,6,8)
dat <- c(62,60,63,59,NA,NA,NA,NA,
         63,67,71,64,65,66,NA,NA,
         68,66,71,67,68,68,NA,NA,
         56,62,60,61,63,64,63,59)
x <- matrix(dat,max(n),r)
cat("\n")
cat("Based on A. Gelman, J.B. Carlin, H.S. Stern and D.B. Rubin    \n")
cat("Bayesian Data Analysis, London: Chapman & Hall 1995, Sec. 9.8 \n")
cat("\n")
N <- sum(n)
xidot <- rep(0,r)
ssi <- rep(0,r)
for (i in 1:r){
  xidot[i] <- sum(x[,i],na.rm=TRUE)/n[i]
  ssi[i] <- (n[i]-1)*var(x[,i],na.rm=TRUE)
}
xdotdot <- sum(x,na.rm=TRUE)/N
ssw <- sum(ssi)
ssb <- (r-1)*var(xidot)
mu <- xdotdot
phi <- ssw/(N-1)
psi <- ssb/(r-1)
muold <- mu
phiold <- phi
psiold <- psi
for (t in 1:niter){
  muold <- mu
  phiold <- phi
  psiold <- psi
  mu <- 0
  phi <- 0
  psi <- 0
  v <- 1/(1/psiold + n/phiold)
  theta <- v*(muold/psiold + n*xidot/phiold)
  mu <- mean(theta)
  for (i in 1:r)
    for (j in 1:n[i])
      phi <- phi + (v[i]+(x[j,i]-theta[i])^2)/(N+2)
  psi <- sum(v + (mu-theta)^2)/r
}
for (i in 1:r)
  cat("Theta[",i,"] =",theta[i],"\n")
cat("\n")
cat("mu =",mu,"\n")
cat("phi =",phi,"\n")
cat("psi =",psi,"\n")
cat("\n")