R/plot.demonoid.R

In benmarwick/LaplacesDemon: Complete Environment for Bayesian Inference

###########################################################################
# plot.demonoid                                                           #
#                                                                         #
# The purpose of the plot.demonoid function is to plot an object of class #
# demonoid.                                                               #
###########################################################################

plot.demonoid <- function(x, BurnIn=0, Data=NULL, PDF=FALSE,
     Parms=NULL, ...)
     {
     ### Initial Checks
     if(missing(x)) stop("The x argument is required.")
     if(is.null(Data)) stop("The Data argument is NULL.")
     if(BurnIn >= nrow(x$Posterior1)) BurnIn <- 0
     Stat.at <- BurnIn + 1
     ### Selecting Parms
     if(is.null(Parms)) {Posterior <- x$Posterior1}
     else {
          Parms <- sub("\\[","\\\\[",Parms)
          Parms <- sub("\\]","\\\\]",Parms)
          Parms <- sub("\\.","\\\\.",Parms)
          if(length(grep(Parms[1], colnames(x$Posterior1))) == 0)
               stop("Parameter in Parms does not exist.")
          keepcols <- grep(Parms[1], colnames(x$Posterior1))
          if(length(Parms) > 1) {
               for (i in 2:length(Parms)) {
                    if(length(grep(Parms[i],
                         colnames(x$Posterior1))) == 0)
                         stop("Parameter in Parms does not exist.")
                    keepcols <- c(keepcols,
                         grep(Parms[i], colnames(x$Posterior1)))}}
          Posterior <- as.matrix(x$Posterior1[,keepcols])
          colnames(Posterior) <- colnames(x$Posterior1)[keepcols]}
     if(PDF == TRUE) {
          pdf("LaplacesDemon.Plots.pdf")
          par(mfrow=c(3,3))
          }
     else {par(mfrow=c(3,3), ask=TRUE)}
     ### Plot Parameters
     for (j in 1:ncol(Posterior))
          {
          plot(Stat.at:x$Thinned.Samples,
               Posterior[Stat.at:x$Thinned.Samples,j],
               type="l", xlab="Iterations", ylab="Value",
               main=colnames(Posterior)[j])
          panel.smooth(Stat.at:x$Thinned.Samples,
               Posterior[Stat.at:x$Thinned.Samples,j], pch="")
          plot(density(Posterior[Stat.at:x$Thinned.Samples,j]),
               xlab="Value", main=colnames(Posterior)[j])
          polygon(density(Posterior[Stat.at:x$Thinned.Samples,j]),
               col="black", border="black")
          abline(v=0, col="red", lty=2)
          ### Only plot an ACF if there's > 1 unique values
          if(!is.constant(Posterior[Stat.at:x$Thinned.Samples,j])) {
               z <- acf(Posterior[Stat.at:x$Thinned.Samples,j], plot=FALSE)
               se <- 1/sqrt(length(Posterior[Stat.at:x$Thinned.Samples,j]))
               plot(z$lag, z$acf, ylim=c(min(z$acf,-2*se),1), type="h",
                    main=colnames(Posterior)[j], xlab="Lag",
                    ylab="Correlation")
               abline(h=(2*se), col="red", lty=2)
               abline(h=(-2*se), col="red", lty=2)
               }
          else {plot(0,0, main=paste(colnames(Posterior)[j],
               "is a constant."))}
          }
     ### Plot Deviance
     plot(Stat.at:length(x$Deviance),
          x$Deviance[Stat.at:length(x$Deviance)],
          type="l", xlab="Iterations", ylab="Value", main="Deviance")
     panel.smooth(Stat.at:length(x$Deviance),
          x$Deviance[Stat.at:length(x$Deviance)], pch="")
     plot(density(x$Deviance[Stat.at:length(x$Deviance)]),
          xlab="Value", main="Deviance")
     polygon(density(x$Deviance[Stat.at:length(x$Deviance)]), col="black",
               border="black")
     abline(v=0, col="red", lty=2)
     ### Only plot an ACF if there's > 1 unique values
     if(!is.constant(x$Deviance[Stat.at:length(x$Deviance)])) {
          z <- acf(x$Deviance[Stat.at:length(x$Deviance)], plot=FALSE)
          se <- 1/sqrt(length(x$Deviance[Stat.at:length(x$Deviance)]))
          plot(z$lag, z$acf, ylim=c(min(z$acf,-2*se),1), type="h",
               main="Deviance", xlab="Lag", ylab="Correlation")
          abline(h=(2*se), col="red", lty=2)
          abline(h=(-2*se), col="red", lty=2)
          }
     else {plot(0,0, main="Deviance is a constant.")}
     ### Plot Monitored Variables
     if(is.vector(x$Monitor)) {J <- 1; nn <- length(x$Monitor)}
     else if(is.matrix(x$Monitor)) {
          J <- ncol(x$Monitor); nn <- nrow(x$Monitor)}
     for (j in 1:J)
          {
          plot(Stat.at:nn, x$Monitor[Stat.at:nn,j],
               type="l", xlab="Iterations", ylab="Value",
               main=Data$mon.names[j])
          panel.smooth(Stat.at:nn, x$Monitor[Stat.at:nn,j], pch="")
          plot(density(x$Monitor[Stat.at:nn,j]),
               xlab="Value", main=Data$mon.names[j])
          polygon(density(x$Monitor[Stat.at:nn,j]), col="black",
               border="black")
          abline(v=0, col="red", lty=2)
          ### Only plot an ACF if there's > 1 unique values
          if(!is.constant(x$Monitor[Stat.at:nn,j])) {
               z <- acf(x$Monitor[Stat.at:nn,j], plot=FALSE)
               se <- 1/sqrt(length(x$Monitor[Stat.at:nn,j]))
               plot(z$lag, z$acf, ylim=c(min(z$acf,-2*se),1), type="h",
                    main=Data$mon.names[j], xlab="Lag", ylab="Correlation")
               abline(h=(2*se), col="red", lty=2)
               abline(h=(-2*se), col="red", lty=2)
               }
          else {plot(0,0, main=paste(Data$mon.names[j], "is a constant."))}
          }
     ### Diminishing Adaptation
     if(nrow(x$CovarDHis) > 1) {
          if(x$Algorithm %in% c("Adaptive Hamiltonian Monte Carlo",
               "Hamiltonian Monte Carlo with Dual-Averaging",
               "No-U-Turn Sampler")) {
               plot(x$CovarDHis[,1], type="l", xlab="Adaptations",
                    ylab=expression(epsilon))}
          else {
               Diff <- abs(diff(x$CovarDHis))
               adaptchange <- matrix(NA, nrow(Diff), 3)
               for (i in 1:nrow(Diff)) {
                    adaptchange[i,1:3] <- as.vector(quantile(Diff[i,],
                         probs=c(0.025, 0.500, 0.975)))}
               plot(adaptchange[,2], ylim=c(min(adaptchange), max(adaptchange)),
                    type="l", col="red", xlab="Adaptations",
                    ylab="Absolute Difference", main="Proposal Variance",
                    sub="Median=Red, 95% Bounds=Gray")
               lines(adaptchange[,1], col="gray")
               lines(adaptchange[,3], col="gray")}
          }
     if(PDF == TRUE) dev.off()
     }

#End