r <- 10
phi <- rep(NA,r)
mu <- rep(NA,r)
S <- rep(NA,r)
n <- 12; xbar <- 119; SS <- 13045
phi0 <- 20; mu0 <- 110; S0 <- 2700; nu0 <- 11
S[1] <- S0
nustar <- nu0 + n
for (i in 2:r) {
  phi[i] <- (1/phi0 + n*nustar/S[i-1])^{-1} 
  mu[i] <- phi[i]*(mu0/phi0 + n*xbar*nustar/S[i-1])
  S[i] <- S0 + (SS+n*xbar^2) - 2*n*xbar*mu[i] + 
            n*(mu[i]^2 + phi[i])
  cat("i",i,"phi",phi[i],"mu",mu[i],"S",S[i],"\n")
}
mustar <- mu[r]; phistar <- phi[r]; Sstar <- S[r]
cat("mu has mean",mustar,"and s.d.",sqrt(phistar),"\n")
cat("phi has mean",Sstar/(nustar-2),"and s.d.",
    (Sstar/(nustar-2))*sqrt(2/(nustar-4)),"\n")