# edps 590BAY
# Fall 2019
# C.J.Anderson
#
# Estimate the mean and variance of normaly distributed variable
#
# Illustration of jags input and output for
#          rjags     
#          run.jags 
#          jagsUI    (most verbose output)
#
library(coda)

setwd("C:\\Users\\cja\\Dropbox\\edps 590BAY\\Lectures\\6 Gibbs and jags\\6 Gibbs and jags")

source("C:\\Users\\cja\\Dropbox\\edps 590BAY\\Lectures\\6 Gibbs and jags\\DBDA2E-utilities.txt")

ano <- read.table("anorexia_data.txt",header=T)

ano$change <- ano$weight2 - ano$weight1

############################################################################################
# rjags
############################################################################################
set.seed(234590)

library(rjags)

dataList <- list(y=ano$change, 
                 Ntotal=length(ano$change),
                 meanY = mean(ano$change),
                 sdY = sd(ano$change)
                 )

Model1 = "model {                                 
                 for (i in 1:Ntotal){
                  y[i] ~ dnorm( mu, precision )                
                 }
                 mu ~ dnorm( meanY , 1/(100*sdY^2) )    
                 precision <- pow(sigma, -2)  # = 1/var(y)
                 sigma ~ dunif( sdY/1000, sdY*1000 )                   }
               "

writeLines(Model1, con="Model1.txt")

thetaInit = mean(ano$change)
sigmaInit = sd(ano$change)
initsList = list(mu=thetaInit, sigma=sigmaInit )

jagsModel1 <- jags.model(file="Model1.txt",     # compiles and intializes model
                        data=dataList,
                        inits=initsList,
                        n.chains=4,
                        n.adapt=500,
						mc.cores=4)            # adaptive phase for 500 iterations

update (jagsModel1, n.iter=500)                 # gets samples from posterior with a 
                                                #  burn in of 500 iterations
												
# run samples from all chains with 500 iterations
#Samples <- coda.samples(jagsModel1, variable.names=c("mu","sigma"), n.iter=4000)

# A bit Faster if turn progress bar off
Samples <- coda.samples(jagsModel1, variable.names=c("mu","sigma"), n.iter=4000,
           progress.bar="none")

# output summary information
summary(Samples)
plot(Samples)

### These can be used to compare models
#
dic   <- dic.samples(jagsModel1,n.iter=1000,type="pD",progress.bar="none")

s <- jags.samples(jagsModel1, c("deviance", "WAIC"), type="mean", n.iter=10000, thin=1)
s <- lapply(s, unclass)
sapply(s, sum)


autocorr.plot(Samples,auto.layout=TRUE)
geweke.diag(Samples,frac1=0.1,frac2=0.5)
effectiveSize(Samples)
cumuplot(Samples,probs=c(.25,.50,.75),lwd=c(1,2),lty=c(2,1),col=c("blue","red"))
gelman.plot(Samples)

#
# These can be used to compare models
#
dic   <- dic.samples(jagsModel1,n.iter=1000,type="pD",progress.bar="none")

s <- jags.samples(jagsModel1, c("deviance", "WAIC"), type="mean", n.iter=10000, thin=1)
s <- lapply(s, unclass)
sapply(s, sum)

###########################################################
##  runjags: Alternative way to run model and get output 
###########################################################
library(runjags)

inits1 <- list("mu"=mean(ano$change), "sigma"=sd(ano$change),
                .RNG.name="base::Super-Duper", .RNG.seed=1)
inits2 <- list("mu"=rnorm(1,2,4), "sigma"=1 ,
               .RNG.name="base::Wichmann-Hill", .RNG.seed=2)				 
inits3 <- list("mu"=rnorm(1,4,1), "sigma"=8 ,
               .RNG.name="base::Wichmann-Hill", .RNG.seed=3)
inits4 <- list("mu"=rnorm(1,-4,2),"sigma"=.5,
               .RNG.name="base::Wichmann-Hill", .RNG.seed=4)
			   
initsList <- list(inits1,inits2,inits3,inits4)	
		

out.runjags <- run.jags(model=Model1, 
                        monitor=c("mu","sigma"),
                        data=dataList, 
						n.chains=4, 
						inits=initsList,
						method="parallel"
)

print(out.runjags)

#################################################
#  jagsUI 
#################################################
library(jagsUI)
initsList = list(list("mu"=mean(ano$change), "sigma"=sd(ano$change)),
                 list("mu"=rnorm(1,2,4), "sigma"=1 ),
                 list("mu"=rnorm(1,4,1), "sigma"=8 ),
                 list("mu"=rnorm(1,-4,2),"sigma"=.5 )
                 )

out.jagsUI <- jags(model.file="Model1.txt",
              data=dataList,inits=initsList, 
              parameters.to.save=c("mu","sigma","precision"), 
              n.iter=2000,
              n.burnin=500,
              n.chains=4)
print(out.jagsUI)