####################################################
#
# 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(lattice)
library(coda)
library(BRugs)
library(R2WinBUGS)

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


# preliminaries. set constants.
n.update2=20000
n.burnin=5000
n.chains=2 # note: will need to change init list to match.
modelSetSeed(12) # set the seed
temp.filename="oring2.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))
    P[1]~ dbeta(1.6,1)
    P[2] ~ dbeta(1,1.6)
    beta[1] <- (75/20)*logit(P[1])+(-55/20)*logit(P[2])
    beta[2] <- (-1/20)*logit(P[1])+(1/20)*logit(P[2])
}



# 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(P=c(0.5,0.5))  
model.inits2 <- function() list(P=c(0.1,0.9))  

init.list = list(model.inits1(),model.inits2())

# 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","P")
file.inits1=bugsInits(model.inits1,fileName = tempfile()) # write inits
file.inits2=bugsInits(model.inits2,fileName = tempfile()) # write inits
file.vector=c(file.inits1,file.inits2)
##########
# 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.vector)  # 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 ....
samplesStats("*",firstChain=1,lastChain=1) # the summarized results
samplesStats("*",firstChain=2,lastChain=2) # the summarized results
# Look at graphical info
samplesHistory("*", mfrow = c(3, 2),firstChain=1,lastChain=1) # histories,
samplesHistory("*", mfrow = c(3, 2),firstChain=2,lastChain=2) # histories,
samplesDensity("*",mfrow = c(3, 2),firstChain=1,lastChain=2) # densities
samplesBgr("*",mfrow = c(3, 2)) # bgr statistics
samplesAutoC("*", chain=1,mfrow = c(3, 2)) # autocorrelations (1st chain)
samplesAutoC("*", chain=2,mfrow = c(3, 2)) # autocorrelations (2nd chain)
dicStats() # look at deviance
mcmc.chains = buildMCMC("*",firstChain=1,lastChain=n.chains)
source("summariseMCMC")
summariseMCMC("Prob1","Probability 0-Ring Failure 55")
summariseMCMC("Prob2","Probability 0-Ring Failure 75")

##########
# 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
# Note: gelman plot seems to generate error with oring2.odc
gelman.diag(mcmc.chains,0.95) 
gelman.plot(mcmc.chains) 

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


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


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