R/Juxtapose.R

In benmarwick/LaplacesDemon: Complete Environment for Bayesian Inference

###########################################################################
# Juxtapose                                                               #
#                                                                         #
# The purpose of the Juxtapose function is to compare the inefficiency of #
# multiple updates of LaplacesDemon, each with a different algorithm, but #
# the same model, initial values, and data. The internal ar.act and       #
# ar.act1 functions are slightly modified versions of those in the        #
# SamplerCompare package (modified to use the rmvn function), and the     #
# method of assessing inefficiency differs.                               #
###########################################################################

Juxtapose <- function(x)
     {
     ### Initial Checks
     if(missing(x)) stop("The x argument is required.")
     if(!is.list(x)) stop("The x argument must be a list.")
     ### IAT with probability intervals
     ar.act1  <- function(y)
          {
          stopifnot(NCOL(y) == 1)
          if(length(unique(y)) < 5)
               return(list(act=NA, act.025=NA, act.975=NA, se=NA,
                    order=NA))
          order.max <- NULL
          repeat {
               A <- ar.yw(y, demean=FALSE, order.max=order.max)
               if(A$order == 0)
                    A <- ar.yw(y, demean=FALSE, order.max=1, aic=FALSE)
               pi <- A$ar
               pi.var <- A$asy.var.coef
               if(kappa(pi.var) < 1 / sqrt(.Machine$double.eps) ||
                    isTRUE(order.max == 1))
                    break
               order.max <- floor(A$order * 0.7)}
          acf <- matrix(ARMAacf(ar=pi)[2:(A$order+1)])
          act <- (1-sum(pi*acf))/(1-sum(pi))^2
          simulation.length <- min(max(40, length(y)), 5000)
          pi.var2 <- as.symmetric.matrix(pi.var)
          AX <- rmvn(simulation.length, pi, pi.var2)
          act.sim <- numeric(simulation.length)
          for (i in 1:simulation.length) {
               pi.sim <- AX[i,]
               acf.sim <- ARMAacf(ar=pi.sim)[2:(A$order+1)]
               if(any(abs(polyroot(c(-1,pi.sim)))<1)) act.sim[i] <- Inf
               else act.sim[i] <- (1-sum(pi.sim*acf.sim)) /
                    (1-sum(pi.sim))^2}
          act.sim[is.na(act.sim)] <- Inf
          act.025 <- as.numeric(quantile(act.sim, 0.025))
          act.975 <- as.numeric(quantile(act.sim, 0.975))
          se <- (act.975-act.025)/(2*1.96)
          return(list(act=act, se=se, act.025=act.025, act.975=act.975,
               order=A$order))
          }
     ar.act <- function(Y, true.mean=NULL)
          {
          Y <- as.matrix(Y)
          stopifnot(is.null(true.mean) || ncol(Y)==length(true.mean))
          if(is.null(true.mean)) mu <- colMeans(Y)
          else mu <- true.mean
          acts <- sapply(1:ncol(Y), function(i) ar.act1(Y[,i]-mu[i]))
          max.i <- which.max(unlist(acts['act',]))
          if(length(max.i)!=1) {
               return(list(act=NA, act.025=NA, act.975=NA, se=NA,
                    order=NA))}
          else return(acts[,max.i])
          }
     ### Back to Juxtapose
     lenx <- length(x)
     ### Set up Output
     out <- matrix(NA, 9, lenx)
     algs <- rep(NA, lenx)
     for (i in 1:lenx) {
          if(!identical(class(x[[i]]), "demonoid"))
               stop("A component of x was found not be of class demonoid.")
          ### Use the abbreviated name of the algorithm
          alg <- x[[i]][["Algorithm"]]
          if(alg == "Adaptive Directional Metropolis-within-Gibbs")
               alg <- "ADMG"
          else if(alg == "Adaptive Griddy-Gibbs") alg <- "AGG"
          else if(alg == "Adaptive Hamiltonian Monte Carlo") alg <- "AHMC"
          else if(alg == "Adaptive Metropolis") alg <- "AM"
          else if(alg == "Adaptive-Mixture Metropolis") alg <- "AMM"
          else if(alg == "Adaptive Metropolis-within-Gibbs") alg <- "AMWG"
          else if(alg == "Affine-Invariant Ensemble Sampler") alg <- "AIES"
          else if(alg == "Componentwise Hit-And-Run Metropolis")
               alg <- "CHARM"
          else if(alg == "Delayed Rejection Adaptive Metropolis")
               alg <- "DRAM"
          else if(alg == "Delayed Rejection Metropolis") alg <- "DRM"
          else if(alg == "Differential Evolution Markov Chain")
               alg <- "DEMC"
          else if(alg == "Elliptical Slice Sampling") alg <- "ESS"
          else if(alg == "Experimental") alg <- "Exper"
          else if(alg == "Griddy-Gibbs") alg <- "GG"
          else if(alg == "Hamiltonian Monte Carlo") alg <- "HMC"
          else if(alg == "Hamiltonian Monte Carlo with Dual-Averaging")
               alg <- "HMCDA"
          else if(alg == "Hit-And-Run Metropolis") alg <- "HARM"
          else if(alg == "Independence Metropolis") alg <- "IM"
          else if(alg == "Interchain Adaptation") alg <- "INCA"
          else if(alg == "Metropolis-Adjusted Langevin Algorithm")
               alg <- "MALA"
          else if(alg == "Metropolis-within-Gibbs") alg <- "MWG"
          else if(alg == "No-U-Turn Sampler") alg <- "NUTS"
          else if(alg == "Reversible-Jump") alg <- "RJ"
          else if(alg == "Robust Adaptive Metropolis") alg <- "RAM"
          else if(alg == "Random-Walk Metropolis") alg <- "RWM"
          else if(alg == "Sequential Adaptive Metropolis-within-Gibbs")
               alg <- "SAMWG"
          else if(alg == "Sequential Metropolis-within-Gibbs")
               alg <- "SMWG"
          else if(alg == "Slice Sampler") alg <- "Slice"
          else if(alg == "Stochastic Gradient Langevin Dynamics")
               alg <- "SGLD"
          else if(alg == "Tempered Hamiltonian Monte Carlo") alg <- "THMC"
          else if(alg == "t-walk") alg <- "twalk"
          else if(alg == "Updating Sequential Adaptive Metropolis-within-Gibbs")
               alg <- "USAMWG"
          else if(alg == "Updating Sequential Metropolis-within-Gibbs")
               alg <- "USMWG"
          algs[i] <- alg
          }
     colnames(out) <- algs
     rownames(out) <- c("iter.min","t.iter.min","prop.stat","IAT.025",
          "IAT.500","IAT.975","ISM.025","ISM.500","ISM.975")
     class(out) <- "juxtapose"
     ### Juxtapose Algorithms
     for (i in 1:lenx) {
          iter.min <- x[[i]][["Iterations"]] / x[[i]][["Minutes"]]
          t.iter.min <- x[[i]][["Iterations"]] / x[[i]][["Thinning"]] /
               x[[i]][["Minutes"]]
          if(x[[i]][["Rec.BurnIn.Thinned"]] >= x[[i]][["Thinned.Samples"]])
               prop.stat <- 0
          else prop.stat <- 1 - (x[[i]][["Rec.BurnIn.Thinned"]] /
               x[[i]][["Thinned.Samples"]])
          if(all(is.na(x[[i]][["Summary2"]]))) {
               iat.500 <- iat.025 <- iat.975 <- Inf}
          else {
               iat.temp <- ar.act(x[[i]][["Posterior2"]])
               iat.500 <- iat.temp$act
               iat.025 <- iat.temp$act.025
               iat.975 <- iat.temp$act.975}
          ism.025 <- prop.stat * t.iter.min / iat.975
          ism.500 <- prop.stat * t.iter.min / iat.500
          ism.975 <- prop.stat * t.iter.min / iat.025
          out[1,i] <- round(iter.min,2)
          out[2,i] <- round(t.iter.min,2)
          out[3,i] <- round(prop.stat,2)
          out[4,i] <- round(iat.025,2)
          out[5,i] <- round(iat.500,2)
          out[6,i] <- round(iat.975,2)
          out[7,i] <- round(ism.025,2)
          out[8,i] <- round(ism.500,2)
          out[9,i] <- round(ism.975,2)
          }
     return(out)
     }

#End