###############################################
#WinBUGS Demo
#Dye Yields Data
#2/14/06
###############################################

#Set working directory
setwd("U:\\Teaching\\Stat694\\Computing\\")

#Load BRugs library (use install.packages("BRugs") the first time you use it
library(BRugs)

#Read in the dye data
dye.data <- read.table("dye-yields.dat",header=T,sep='\t')

N.batches <- dim(dye.data)[2]    #Number of batches
N.samples <- dim(dye.data)[1]    #Number of samples of each batch

#Directory with WinBUGS files
model.dir <- "BUGS-Dye-Yields/"

#Write the data to the BUGS directory
data.bugs <- list("y"=as.matrix(dye.data),"N.batches"=N.batches,"N.samples"=N.samples)
bugsData(data.bugs, fileName = file.path(getwd(), paste(model.dir,"data.txt",sep="")))

#Write the initial values to the BUGS directory
inits.bugs <- list("theta"=1500,"tau.with"=1,"tau.btw"=1,"mu"=rep(1500,N.batches))
bugsInits(list(inits.bugs), fileName = file.path(getwd(), paste(model.dir,"init.txt",sep="")))

#WinBUGS setup
modelCheck(fileName = file.path(getwd(), paste(model.dir,"model.txt",sep="")))          # check model file
modelData(file.path(getwd(), paste(model.dir,"data.txt",sep="")))                       # read data file
modelCompile(numChains=1)                                                               # compile model with 1 chain
modelInits(rep(file.path(getwd(), paste(model.dir,"init.txt",sep="")),1))               # read init data file

#Set nodes (parameters to monitor)
samplesSet(c("theta","mu","tau.with","tau.btw","sigma2.with","sigma2.btw"))

#Run the chain for 50 iteration and look at trace plots
ITER <- 50
modelUpdate(ITER)
samplesHistory('*')

#Run the chain for 25000 more iteration and look at trace plots
ITER <- 25000
modelUpdate(ITER)
samplesHistory('*')

#Set BURNIN and thin the chain
samplesSetBeg(5000)
samplesSetThin(10)
samplesHistory('*')

#Look at density plots
samplesDensity('*')

#Look at sample statistics
samplesStats('*')

#Get samples from the posterior distribution of the total variance
samples.sigma2.with <- samplesSample('sigma2.with')
samples.sigma2.btw <- samplesSample('sigma2.btw')
samples.sigma2.total <- samples.sigma2.with+samples.sigma2.btw
par(mfrow=c(1,1))
hist(samples.sigma2.total)

#Clear samples
samplesClear('*')

#Write another set of initial values
inits2.bugs <- list("theta"=0,"tau.with"=1,"tau.btw"=1,"mu"=rep(0,N.batches))
bugsInits(list(inits2.bugs), fileName = file.path(getwd(), paste(model.dir,"init2.txt",sep="")))

#WinBUGS setup -- 2 chains
modelCheck(fileName = file.path(getwd(), paste(model.dir,"model.txt",sep="")))
modelData(file.path(getwd(), paste(model.dir,"data.txt",sep="")))
modelCompile(numChains=2)
modelInits(c(file.path(getwd(), paste(model.dir,"init.txt",sep="")),file.path(getwd(), paste(model.dir,"init2.txt",sep=""))))

#Set nodes (parameters to monitor)
samplesSet(c("theta","mu","tau.with","tau.btw","sigma2.with","sigma2.btw"))

#Run the chain for 25000 iteration and look at trace plots
ITER <- 25000
modelUpdate(ITER)
samplesSetBeg(5000)
samplesSetThin(10)
samplesHistory('*')