#
#  Rat data in Chapter 9, Exercise 11
#
#  Remember to load the file wishart.r first
#
p <- P <- 2
m <-  500
k <- 30
ni <- 5
epsilon <- 0.001
x <- c(8, 15, 22, 29, 36)
dat <- read.table("rats.dat")
y <- dat[,2:(ni+1)]
alpha0 <- 0.0
beta0 <- 0.0
a <- aalpha <- abeta <- epsilon    #  B P Carlin and T A Louis p. 169
b <- balpha <- bbeta <- 1/epsilon  #  B P Carlin and T A Louis p. 170
sigma2 <- 1.0                      #  Initially sigma2 is IG(a,b)
sigmaa2 <- 100                     #  A E Gelfand et al. p. 979 col.1
sigmab2 <- 0.1                     #  A E Gelfand et al. p. 979 col.1
#  Thus R = (100   0 )
#           ( 0   0.1)
#
#  Take values for alpha[i] and beta[i] given
#  alpha0, beta0, sigmaa2, sigmab2 and sigma2
alphabar <- 0.0
betabar <- 0.0
alpha <- rep(0,k)
beta <- rep(0,k)
for (i in 1:k)
  {
     vara <- ni/sigma2 + 1/sigmaa2
     suma <- sum(y[i,])
     meana <- (suma/sigma2 + alpha0/sigmaa2)/vara
     #  alpha[i] ~ N(meana,vara)
     alpha[i] <- meana+sqrt(vara)*rnorm(1)
     alphabar <- alpha[i]/k
     varb <- var(x)/sigma2 + 1/sigmab2
     sumb <- sum((x-mean(x))*unlist(y[i,]))
     meanb <- (sumb/sigma2 + 1/sigmab2)/varb
     #  beta[i] ~ N(meanb,varb)
     beta[i] <- meanb+sqrt(varb)*rnorm(1)
     betabar <- betabar + beta[i]/k
  }
#  Initialize var (capital sigma)
v <- matrix(c(sigmaa2,0,0,sigmab2),P,P)
#  Take values for alpha0 and beta0 given
#  alpha[i], beta[i], sigmaa2, sigmab2 and sigma2
wish(p,k,var)
#  alpha0 ~ N(alphabar,sigmaa2/k)
alpha0 <- alphabar+sqrt(sigmaa2/k)*rnorm(1)
#  beta0  ~ N(betabar,sigmab2/k)
beta0 <- betabar+sqrt(sigmab2/k)*rnorm(1)
#  See last displayed formula on p. 168 of
#  B P Carlin and T A Louis
#  simplified by taking C^{-1} = 0
#
#  Take values for sigmaa2 and sigmab2 given
#  alpha[i], beta[i], alpha0, beta0 and sigma2
#  sigmaa2 ~ IG(alpha0,beta0)
sigmaa2 <- 1/(beta0*rgamma(1,alpha0))
#  sigmab2 ~ IG(alpha)
sigmab2 <- 1/(beta0*rgamma(1,alpha0))
#
#  Take value for sigma2 given
#  alpha[i], beta[i], alpha0, beta0, sigmaa2 and sigmab2
#  sigma2 ~ IG(alpha0,beta0)
sigma2 <- 1/(beta0*rgamma(1,alpha0))