#####################################################################################
# Code for Gibbs Sampler for the Thermonuclear Test Example in the article:         #
# "Critical Differences in Bayesian and Non-Bayesian Inference."                    #
# In \emph{Current Methodological Developments of Statistics in the Social Sciences #
# Stanislav Kolenikov, Lori Thombs, and Douglas Steinley (eds.).  Forthcoming 2009. #
# John Wiley & Sons.                                      Copyright 2008, Jeff Gill #
#####################################################################################

read.table("http://jgill.wustl.edu/data/nukes.vec.data",header=FALSE)

# R CODE FOR PRIOR -> POSTERIOR CONJUGATE ANALYSIS WITH GRAPHS

par(mfrow=c(1,2),lwd=1.5,mar=c(5,5,2,2),oma=c(1,1,1,1),cex.lab=1.5)
ruler <- seq(0,100,length=500)

alpha <- 50; beta=2; n <- nrow(nukes)
prior.sims <- rgamma(10000,alpha,rate=beta)
post.sims <- rgamma(10000,shape=alpha+n*mean(nukes[,2]),rate=beta+n)
plot(ruler, dgamma(ruler,shape=alpha+n*mean(nukes[,2]),rate=beta+n), type="l", lwd=3,
        xlim=c(0,28), xlab=expression(mu),ylab="Density" )
lines(ruler, dgamma(ruler,alpha,rate=beta), lwd=4, col="grey")

alpha <- 5; beta=1; n <- nrow(nukes)
prior.sims <- rgamma(10000,alpha,rate=beta)
post.sims <- rgamma(10000,shape=alpha+n*mean(nukes[,2]),rate=beta+n)
plot(ruler, dgamma(ruler,shape=alpha+n*mean(nukes[,2]),rate=beta+n), type="l", lwd=3,
        xlim=c(0,28), xlab=expression(mu),ylab="" )
lines(ruler, dgamma(ruler,alpha,rate=beta), lwd=4, col="grey")

# R CODE FOR THE CHANGEPOINT ANALYSIS GIBBS SAMPLER 

bcp <- function(theta.matrix,y,a,b,g,d)  {
    n <- length(y)
    k.prob <- rep(0,length=n)
    for (i in 2:nrow(theta.matrix))  {
        if (i %% 1000 == 0) print(i)
        lambda <- rgamma(1,a+sum(y[1:theta.matrix[(i-1),3]]), b+theta.matrix[(i-1),3])
        phi    <- rgamma(1,g+sum(y[theta.matrix[(i-1),3]:n]), d+length(y)-theta.matrix[(i-1),3])
        for (j in 1:n)  {
                prod1 <- j*(phi-lambda)
                prod2 <- sum(y[1:j])*log(lambda)
                prod3 <- sum(y[1:j])*log(phi)
                k.prob[j] <- exp(prod1 + prod2 - prod3)
        }
        k.prob <- (k.prob/sum(k.prob))
        k      <- sample(1:n,size=1,prob=k.prob)
        theta.matrix[i,] <- c(lambda,phi,k)
    }
    return(theta.matrix)
}

num.reps <- 1000000
nukes.mat <- matrix(NA,ncol=3,nrow=num.reps)
nukes.mat[1,] <- c(20,10,20)
alpha <- 1; beta <- 1; gamma <- 1; delta <- 1
nukes.mat <- bcp(nukes.mat,nukes[,2],alpha,beta,gamma,delta)
#nukes.mat <- bcp(nukes.mat,nukes[1:36,2],alpha,beta,gamma,delta)

summary(nukes.mat[500000:1000000,])

# BUGS CODE FOR JAGS MODEL WITH COVARIATES
# R SETUP:
nukes2 <- read.table("http://jgill.wustl.edu/data/nukes.data",header=TRUE)
nukes2.df <- data.frame("US.TEST"=nukes[,2],nukes2)

# GET THIS FUNCTION FROM THE BUGS SITE
source("writeDatafileR.R")	
# MODIFIED BY HAND LATER FOR JAGS, THERE ARE NOW DIRECT FUNCTIONS FOR JAGS FORMATTED DATA
 writeDatafileR(nukes2.df,towhere="Chapter.Winemiller/nukes2.data") 

# BUGS MODEL:

model {
    for (i in 1:n)  {
        log(mu[i]) <- b0 + b2*US.GDP[i] + b3*US.DEBT[i] + b4*(SOV.TEST[i]-mean(SOV.TEST))
        US.TEST[i] ~ dpois(mu[i])
    }

    b0 ~ dnorm(0,1.0E-4)
    b2 ~ dnorm(0,1.0E-4)
    b3 ~ dnorm(0,1.0E-4)
    b4 ~ dunif(-100,100)
}

# DATA:

"n" <- 40
"US.TEST" <- c(11, 6, 18, 18, 32, 77, 0, 0, 10, 98, 47, 47, 39, 48, 42, 56, 46, 39, 24, 27, 24, 23, 22, 21, 20, 21, 16, 17, 17, 19, 19, 20, 18, 15, 15, 15, 12, 9, 8, 6) 
"YEAR" <- c(1953, 1954, 1955, 1956, 1957, 1958, 1959, 1960, 1961, 1962, 1963, 1964, 1965, 1966, 1967, 1968, 1969, 1970, 1971, 1972, 1973, 1974, 1975, 1976, 1977, 1978, 1979, 1980, 1981, 1982, 1983, 1984, 1985, 1986, 1987, 1988, 1989, 1990, 1991, 1992) 
"US.GDP" <- c(0.37, 0.38, 0.4, 0.43, 0.45, 0.46, 0.49, 0.52, 0.53, 0.57, 0.6, 0.64, 0.69, 0.75, 0.81, 0.87, 0.95, 1.01, 1.08, 1.18, 1.31, 1.44, 1.55, 1.73, 1.97, 2.21, 2.5, 2.72, 3.05, 3.21, 3.42, 3.81, 4.1, 4.37, 4.61, 4.95, 5.35, 5.68, 5.86, 6.14) 
"US.DEBT" <- c(0.49, 0.51, 0.55, 0.58, 0.6, 0.64, 0.69, 0.72, 0.77, 0.82, 0.88, 0.94, 1.01, 1.07, 1.15, 1.24, 1.33, 1.42, 1.56, 1.71, 1.9, 2.07, 2.26, 2.51, 2.83, 3.21, 3.6, 3.96, 4.35, 4.77, 5.34, 6.12, 7.09, 7.95, 8.67, 9.44, 10.18, 10.87, 11.35, 11.9) 
"SOV.TEST" <- c(0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 9, 14, 18, 17, 17, 19, 16, 23, 24, 17, 21, 19, 21, 24, 31, 31, 24, 21, 19, 25, 27, 10, 0, 23, 16, 7, 1, 0, 0)

# INITIAL VALUES:

"b0" <- 0
"b2" <- 0
"b3" <- 0
"b4" <- 0

# JAGS COMMANDS:
model in "nukes2.bug"
data in "nukes2.data"
compile
inits in "nukes2.init"
initialize
update 500000
monitor set b0
monitor set b2
monitor set b3
monitor set b4
update 1000000
coda *
exit