# Edps 590BAY
# Fall 2019
# C.J.Anderson
#
# To check to make sure that you have everything in R set up to
#  run Gibbs sampling via rjags, below is a test program that runs
#  fine.   You just need to change your working directory to where
#  the data set lives and change director to where DBDS2E-utilities.txt lives.
#
# Getting what paid for JAGS.....mean and variance
#
#

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

sat <- read.table("C:\\Users\\cja\\Dropbox\\edps 590BAY\\Lectures\\6 Gibbs and jags\\getting_what_paid_for_data.txt",header=T)

dataList <- list(y=sat$ave_tot,           
				 N=length(sat$ave_tot),
                 sdY = sd(sat$ave_tot),
                 meanY = mean(sat$ave_tot)
                )
mod1 = "model {  for (i in 1:N){
# Likelihood          
                    y[i] ~ dnorm(mu, precision)
                      }
                    mu ~ dnorm(meanY,1/(100*sdY^2))       
                    sigma ~ dunif( 1E-3, 1E+30 )
                    precision <- 1/sigma^2           
              } "

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

Ymean.int = mean(sat$ave_tot)
sigma.int = sd(sat$ave_tot)
initsList = list(mu=Ymean.int, sigma=sigma.int )

jagsMod1 <- jags.model(file="mod1.txt",  
        data=dataList,
        inits=initsList,
        n.chains=4,
        n.adapt=1000)           

update (jagsMod1, n.iter=1000)   

# gets samples from all chains for 2000 iterations
Samp1 <- coda.samples(jagsMod1, variable.names=c("mu","sigma","precision"), n.iter=2000)

# output summary information
summary(Samp1)
plot(Samp1)
#
# these functions are from the DBSDA2E-utilities.txt 
#
par(ask=T)

diagMCMC(codaObject=Samp1, parName="mu")

diagMCMC(codaObject=Samp1, parName="sigma")

diagMCMC(codaObject=Samp1, parName="precision")

plotPost( Samp1[,'mu'], main="mu", xlab=bquote(mu), cenTend="median", compVal=945, credMass=0.95)

plotPost( Samp1[,'sigma'], main="Sigma", xlab=bquote(Sigma), cenTend="median", compVal=76, credMass=0.95)