## This file shows some later post-processing, namely making some
## JAGS/Stan comparision figures that appear in the plot.

library(coda)
## Read back in the data for the two runs I'm re-analyzing.
## Note R will decompress on the fly (these loads take a while though).
source("K2.Simulation.JAGS.2.R.gz")
source("K2.Simulation.Stan.unordered.2.R.gz")

## Compare Deviance/Log Parameter (this is slightly different in Stan
## and JAGS).
png("DeviancePlots.png")
par(mfrow=c(2,2),oma=c(0,0,3,0))
plot(K2.Simulation.JAGS.2$scalars$deviance,auto.layout=FALSE,
     main="JAGS")
plot(K2.Simulation.Stan.unordered.2$scalars$lp__,auto.layout=FALSE,
     main="Stan (unconstrained model)")
mtext("Deviance (JAGS)/Log Posterior (Stan)",outer=TRUE,line=1)
mtext(paste("Rhat (JAGS) =",
         round(gelman.diag(K2.Simulation.JAGS.2$scalars$deviance)$psrf[1,1],2),
         "      effective sample size (JAGS) = ",
         round(effectiveSize(K2.Simulation.JAGS.2$scalars$deviance))),
      outer=TRUE,line=0)
mtext(paste("Rhat (Stan) =",
            round(gelman.diag(K2.Simulation.Stan.unordered.2$scalars$lp__)$psrf[1,1],2),
            "     effective sample size (Stan) = ",
            round(effectiveSize(K2.Simulation.Stan.unordered.2$scalars$lp__))),
      outer=TRUE,line=-1)
dev.off()


png("alpha0Plots.png")
par(mfrow=c(2,2),oma=c(0,0,3,0))
plot(K2.Simulation.JAGS.2$hparams$alpha0[,1],auto.layout=FALSE,
     main="JAGS")
plot(K2.Simulation.Stan.unordered.2$hparams$alpha0[,1],auto.layout=FALSE,
     main="Stan (unconstrained model)")
mtext(expression(alpha[0]),outer=TRUE,line=1)
mtext(paste("Rhat (JAGS) =",
         round(gelman.diag(K2.Simulation.JAGS.2$hparams$alpha0[,1])$psrf[1,1],2),
         "      effective sample size (JAGS) = ",
         round(effectiveSize(K2.Simulation.JAGS.2$hparams$alpha0[,1]))),
      outer=TRUE,line=0)
mtext(paste("Rhat (Stan) =",
            round(gelman.diag(K2.Simulation.Stan.unordered.2$hparams$alpha0[,1])$psrf[1,1],2),
            "     effective sample size (Stan) = ",
            round(effectiveSize(K2.Simulation.Stan.unordered.2$hparams$alpha0[,1]))),
      outer=TRUE,line=-1)
dev.off()


png("mu0Plots.png")
par(mfrow=c(2,2),oma=c(0,0,3,0))
plot(K2.Simulation.JAGS.2$hparams$mu0[,1],auto.layout=FALSE,
     main="JAGS")
plot(K2.Simulation.Stan.unordered.2$hparams$mu0[,1],auto.layout=FALSE,
     main="Stan (unconstrained model)")
mtext(expression(mu[0][1]),outer=TRUE,line=1)
mtext(paste("Rhat (JAGS) =",
         round(gelman.diag(K2.Simulation.JAGS.2$hparams$mu0[,1])$psrf[1,1],2),
         "      effective sample size (JAGS) = ",
         round(effectiveSize(K2.Simulation.JAGS.2$hparams$mu0[,1]))),
      outer=TRUE,line=0)
mtext(paste("Rhat (Stan) =",
            round(gelman.diag(K2.Simulation.Stan.unordered.2$hparams$mu0[,1])$psrf[1,1],2),
            "     effective sample size (Stan) = ",
            round(effectiveSize(K2.Simulation.Stan.unordered.2$hparams$mu0[,1]))),
      outer=TRUE,line=-1)
dev.off()

png("gamma0Plots.png")
par(mfrow=c(2,2),oma=c(0,0,3,0))
plot(K2.Simulation.JAGS.2$hparams$gamma0[,1],auto.layout=FALSE,
     main="JAGS")
plot(K2.Simulation.Stan.unordered.2$hparams$gamma0[,1],auto.layout=FALSE,
     main="Stan (unconstrained model)")
mtext(expression(gamma[0][1]),outer=TRUE,line=1)
mtext(paste("Rhat (JAGS) =",
         round(gelman.diag(K2.Simulation.JAGS.2$hparams$gamma0[,1])$psrf[1,1],2),
         "      effective sample size (JAGS) = ",
         round(effectiveSize(K2.Simulation.JAGS.2$hparams$gamma0[,1]))),
      outer=TRUE,line=0)
mtext(paste("Rhat (Stan) =",
            round(gelman.diag(K2.Simulation.Stan.unordered.2$hparams$gamma0[,1])$psrf[1,1],2),
            "     effective sample size (Stan) = ",
            round(effectiveSize(K2.Simulation.Stan.unordered.2$hparams$gamma0[,1]))),
      outer=TRUE,line=-1)
dev.off()