#################################################################################
#
#  Unit 6 - Example i.
#
#  Hierarchical model in R2WinBugs
#
#
#################################################################################
#
# Author : Hedibert Freitas Lopes
#          Graduate School of Business
#          University of Chicago
#          5807 South Woodlawn Avenue
#          Chicago, Illinois, 60637
#          Email : hlopes@ChicagoGSB.edu
#
#################################################################################
packages("R2WinBUGS")
library("R2WinBUGS")

#################### BEGIN: save in "hierarchicalmodel.bug"
model;
{
   for( j in 1 : J ) {
      theta[j] ~ dnorm(mu.theta,tau.theta)
   }
   for( j in 1 : J ) {
      y[j] ~ dnorm(theta[j],tau.y[j])
   }
   mu.theta ~ dnorm( 0.0,1.0E-6)
   tau.theta <- pow(sigma.theta, -2)
   for( j in 1 : J ) {
      tau.y[j] <- pow(sigma.y[j], -2)
   }
   sigma.theta ~ dunif(0,100)
}
#################### END: save in "hierarchicalmodel.bug"

data(schools)
schools

J       <- nrow(schools)
y       <- schools$estimate
sigma.y <- schools$sd
data    <- list("J","y","sigma.y")
inits   <- function(){
  list(theta=rnorm(J,0,100),mu.theta=rnorm(1,0,100),sigma.theta=runif(1,0,100))
}

schools.sim = bugs(data,inits,model.file="hierarchicalmodel.bug",
                   parameters=c("theta","mu.theta","sigma.theta"),
                   n.chains=3,n.iter=2000,n.burnin=1000,n.thin=1,
                   bugs.directory="c:/Program Files/WinBUGS14/",codaPkg=FALSE)

print(schools.sim)

plot(schools.sim)