####################################################
#
# Bayesian Analysis Using BRugs
#
# A Primer on using BRugs to analyze pumps.odc.
#
# References:
# Title: Bayesian Data Analysis
# Author: Andrew Gelman, John Carlin, Hal Stern, Donald Rubin
# Edition: Second
#
# http://web.maths.unsw.edu.au/~ahayen/Brugs/price.html
# http://www.stat.umn.edu/geyer/5102/examp/reg.html
# http://www.stat.columbia.edu/~gelman/bugsR/
# 
# Change Log:
# Author: Cosenza, Carlo
# Last Changed: 8/17/2007
# 
####################################################


# link to required libraries
library(coda)
library(lattice)
library(BRugs)
#library(R2WinBUGS)

##########
# START  #
##########


# preliminaries. set constants.
bd <- "c:/Program Files/WinBUGS14/"
n.updates=4000
n.burnin=3000
n.chains=2 # note: will need to change init list to match.
modelSetSeed(1237) # set the seed

##########
# Step 1 #
##########
# Define model 
BRugsmodel <- function()
{
    for (i in 1 : N) {
        theta[i] ~ dgamma(alpha, beta)
        lambda[i] <- theta[i] * tt[i]
        x[i] ~ dpois(lambda[i])
    }
    alpha ~ dexp(1)
    beta ~ dgamma(0.1, 1.0)
}

# Define data
# Note: the variable 't' generated errors. change to 'tt'
model.data <- list(tt = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5),
                   x = c(5,    1,    5,    14,  3,    19,   1,    1,    4,   22), 
                   N = 10)
# Define inits
# Note: BRugs requires theta be initialized or generates error.
model.inits <- function() list(alpha = 1, beta = 1,theta=runif(10,0,1))
init.list = list(model.inits(),model.inits()) # 2 inits for 2 chains

# write model, data, and inits to text files.
model.file <- file.path(tempdir(), "pumpmodel.txt")   # temp filename
write.model(BRugsmodel, model.file)                   # write model
file.data <- bugsData(model.data,fileName=tempfile()) # write data
file.inits=bugsInits(init.list,fileName = tempfile()) # write inits

# List Parameters that will be monitored
model.parameters <- c("alpha", "beta", "theta")

##########
# Step 2 #
##########
# Run BUGS
modelCheck(model.file)             # check model file
modelData(file.data)               # read data file
modelCompile(numChains=n.chains)   # compile model with 1 chain
modelInits(file.inits,chainNum=1)  # set inits
modelInits(file.inits,chainNum=2)  # set inits
modelUpdate(n.burnin)              # Run Burn-in
samplesSet(model.parameters)       # set parameters to be monitored
dicSet()                           # Set DIC assuming model has converged
modelUpdate(n.updates)             # more iterations ....

##########
# Step 3 #
##########

# Set to observe all parameters
samplesStats("*") # the summarized results
dicStats() # look at deviance

# Look at graphical info
samplesHistory("*", mfrow = c(4, 2)) # histories,
samplesDensity("*") # densities
samplesBgr("*") # bgr statistics
samplesAutoC("*", 1) # autocorrelations (1st chain)

##########
# Step 4 #
##########
# get Markov Chains
mcmc.chains = buildMCMC("*",firstChain=1,lastChain=n.chains)

plot(mcmc.chains) # view trace and density plot for both chains
gelman.diag(mcmc.chains,0.95) # converges when upper limit is close to 1
gelman.plot(mcmc.chains) # view gelman plots

geweke.diag(mcmc.chains)
geweke.plot(mcmc.chains)

heidel.diag(mcmc.chains) 
raftery.diag(mcmc.chains)


HPDinterval(mcmc.chains[1]) # display Highest Posterior Density intervals

# some other plots
cumuplot(mcmc.chains[1], probs=c(0.025,0.975))
crosscorr(mcmc.chains[1])
autocorr(mcmc.chains[1])
qqmath(mcmc.chains[1])
densityplot(mcmc.chains[1])
acfplot(mcmc.chains[1])
plot(mcmc.chains[2])
HPDinterval(mcmc.chains[2])
densplot(mcmc.chains[2])

##########
# Step 4 #
##########
# Clear monitor and free memory
dicClear()
samplesClear("*") # clear monitor

##########
# FINISH #
##########