#
# Hierachical normal model in Chapter 9, Exercise 6
#
niter <- 25
r <- 4
n <- c(4,4,4,4)
dat <- c(
   98,97,99,96,
   91,90,93,92,
   96,95,97,95,
   95,96,99,98)
x <- matrix(dat,max(n),r)
cat("\n")
cat("Data quoted in P.M. Lee, Bayesian Statistics: An Introduction \n")
cat("(2nd edn), London: Arnold 1997, Chapter 9, Exercise 6.        \n")
cat("\n")
N <- sum(n)
xidot <- rep(0,r)
ssi <- rep(0,r)
for (i in 1:r){
  xidot[i] <- sum(x[1:n[i],i])/n[i]
  ssi[i] <- (n[i]-1)*var(x[1:n[i],i])
}
xdotdot <- sum(x)/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")