####################################################
#
# Bayesian Analysis Using BRugs
#
# A Primer on using BRugs to analyze oring1.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/18/2007
# 
####################################################


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

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


# preliminaries. set constants.
n.updates=20000
n.burnin=5000
n.chains=2 # note: will need to change init list to match.
modelSetSeed(1237) # set the seed
temp.filename="oring1.txt"
##########
# Step 1 #
##########
# Define model 
BRugsmodel <- function()
{ 
    for(i in 1:n){
        y[i] ~ dbern(p[i])
        logit(p[i]) <- beta[1] + beta[2]*T[i]
    }
    Prob1 <- exp(beta[1] + beta[2]*55)/(1+exp(beta[1]+beta[2]*55))
    Prob2 <- exp(beta[1] + beta[2]*75)/(1+exp(beta[1]+beta[2]*75))
    for(j in 1:2){
        beta[j] ~ dnorm(0,0.25)
    }
}

# Define data
model.data <- list(T= c(53,57,58,63,66,67,
                        67,67,68,69,70,70,
                        70,70,72,73,75,75,
                        76,76,78,79,81), 
                   y= c(1,1,1,1,0,0,
                        0,0,0,0,0,0,
                        1,1,0,0,0,1,
                        0,0,0,0,0),
                   n=23)

# Define inits
model.inits1 <- function() list(beta=c(0,0))  
model.inits2 <- function() list(beta=c(-1,1)) 
file1.inits=bugsInits(model.inits1,fileName = tempfile()) # write inits
file2.inits=bugsInits(model.inits2,fileName = tempfile()) # write inits
file.inits=c(file1.inits,file2.inits)

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


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

##########
# Step 2 #
##########
# Run BUGS
modelCheck(model.file)             # check model file
modelData(file.data)               # read data file
modelCompile(numChains=n.chains)   # compile model n chains
modelInits(file.inits)  # 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("*",firstChain=1,lastChain=1) # the summarized results
samplesStats("*",firstChain=2,lastChain=2) # the summarized results

dicStats() # look at deviance
# Look at graphical info
samplesHistory("*", mfrow = c(3, 2)) # histories,
samplesDensity("*",mfrow = c(2, 2)) # densities
samplesBgr("*",mfrow = c(2, 2)) # bgr statistics
samplesAutoC("*", chain=1,mfrow = c(2, 2)) # autocorrelations (1st chain)
samplesAutoC("*", chain=2,mfrow = c(2, 2)) # autocorrelations (2nd chain)
mcmc.chains = buildMCMC("*",firstChain=1,lastChain=n.chains)

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

# analyze Markov 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)


HPDinterval(mcmc.chains) # display Highest Posterior Density intervals
summary(mcmc.chains) # general display


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