R/jagsposteriors.R

In RafaelSdeSouza/nuclear: Functions to estimate S-factor for nuclear reactions

#' Print results for specified rjags or mcmc.list parameters
#'
#' Prints results for \code{rjags} or \code{mcmc.list} parameters, which are
#' specified with a regular expression, or by exact name.
#'
#' @param x The \code{rjags}, \code{rjags.parallel}, or \code{mcmc.list} object
#'   for which results will be printed.
#' @param params Character vector or a regular expression pattern. The
#'   parameters for which results will be printed (unless \code{invert} is
#'   \code{FALSE}, in which case results for all parameters other than those
#'   given in \code{params} will be returned). If \code{regex} is \code{FALSE},
#'   only those parameters that match \code{params} exactly will be returned. If
#'   \code{regex} is \code{TRUE}, \code{param} should be a character string
#'   giving the regular expression pattern to be matched.
#' @param regex Logical. If \code{regex} is \code{TRUE}, then \code{param} is
#'   expected to be a single string giving a text pattern to be matched.
#'   Parameters with names matching the pattern will be returned (unless
#'   \code{invert} is \code{TRUE}, which results in all parameters that do not
#'   match the pattern being returned). Text pattern matching uses regular
#'   expressions (\code{\link{regex}}).
#' @param invert Logical. If \code{invert} is \code{TRUE}, only those parameters
#'   that do not match elements of \code{params} will be returned.
#' @param probs A numeric vector of probabilities within range [0, 1],
#'   representing the sample quantiles to be calculated and returned.
#' @param signif If supplied, all columns other than \code{n.eff} will have
#'   their values rounded such that the most extreme value has the specified
#'   number of significant digits.
#' @param ... Additional arguments accepted by \code{\link{grep}}, e.g.
#'   \code{perl=TRUE}, to allow look-around pattern matching.
#' @return A matrix with one row for each parameter that matches \code{param},
#'   and one column for each of \code{mean}, \code{sd}, percentiles \code{2.5},
#'   \code{25}, \code{50}, \code{75}, and \code{97.5}. In addition, if \code{x}
#'   is an \code{rjags} object, columns for \code{Rhat} and \code{neff} are
#'   returned.
#' @seealso \code{\link{rhats}} for simplified rhat output,
#'   \code{\link{rearray}} for recovering array structure of vector-valued
#'   parameters.
#' @export
#' @examples
#' ## Data
#' N <- 100
#' temp <- runif(N)
#' rain <- runif(N)
#' wind <- runif(N)
#' a <- 0.13
#' beta.temp <- 1.3
#' beta.rain <- 0.86
#' beta.wind <- -0.44
#' sd <- 0.16
#' y <- rnorm(N, a + beta.temp*temp + beta.rain*rain + beta.wind*wind, sd)
#' dat <- list(N=N, temp=temp, rain=rain, wind=wind, y=y)
#'
#' ### jags example
#' library(R2jags)
#'
#' ## Model
#' M <- function() {
#'   for (i in 1:N) {
#'     y[i] ~ dnorm(y.hat[i], sd^-2)
#'     y.hat[i] <- a + beta.temp*temp[i] + beta.rain*rain[i] + beta.wind*wind[i]
#'     resid[i] <- y[i] - y.hat[i]
#'   }
#'   sd ~ dunif(0, 100)
#'   a ~ dnorm(0, 0.0001)
#'   beta.temp ~ dnorm(0, 0.0001)
#'   beta.rain ~ dnorm(0, 0.0001)
#'   beta.wind ~ dnorm(0, 0.0001)
#' }
#'
#' ## Fit model
#' jagsfit <- jags(dat, inits=NULL,
#'                 parameters.to.save=c('a', 'beta.temp', 'beta.rain',
#'                                      'beta.wind', 'sd', 'resid'),
#'                 model.file=M, n.iter=10000)
#'
#' ## Output
#' # model summary
#' jagsfit
#'
#' # Results for beta.rain only
#' jagsposteriors(x=jagsfit, param='beta.rain')
#'
#' # Results for 'a' and 'sd' only
#' jagsposteriors(x=jagsfit, param=c('a', 'sd'))
#' jagsposteriors(x=jagsfit, param=c('a', 'sd'),
#'             probs=c(0.01, 0.025, 0.1, 0.25, 0.5, 0.75, 0.9, 0.975))
#'
#' # Results for all parameters including the string 'beta'
#' jagsposteriors(x=jagsfit, param='beta', regex=TRUE)
#'
jagsposteriors <- function(x, params, regex=FALSE, invert=FALSE,
                        probs=c(0.025, 0.25, 0.5, 0.75, 0.975), signif, ...) {
  if(!regex) {
    params <- paste0(gsub('(?=\\.|\\[|\\])', '\\1\\\\', params, perl=TRUE),
                     '(\\[.*\\])?', collapse='|')
    params <- paste("^", gsub("\\|", "$|^", params), '$', sep = "")
  } else if(length(params) > 1) {
    stop("If 'regex' is TRUE, 'params' must be a single regex string.",
         call.=FALSE)
  }
  if(any(is(x) %in% c('rjags.parallel', 'rjags'))) {
    nm <- dimnames(x$BUGSoutput$sims.array)[[3]]
    i <- grep(params, nm, invert=invert, ...)
    if(length(i) == 0) stop("No parameters match 'params'", call.=FALSE)
    samp <- x$BUGSoutput$sims.array[, , i, drop=FALSE]
    rhat_neff <- x$BUGSoutput$summary[i, c('Rhat', 'n.eff'), drop=FALSE]
    out <- cbind(t(apply(
      samp, 3, function(x)
        c(mean=mean(x), sd=sd(x), quantile(x, probs=probs)))), rhat_neff)
  } else if(any(is(x)=='mcmc.list')) {
    nm <- colnames(x[[1]])
    i <- grep(params, nm, invert=invert, ...)
    if(length(i) == 0) stop("No parameters match 'params'", call.=FALSE)
    out <- t(apply(do.call(rbind, x), 2, function(z) {
      c(mean=mean(z), sd=sd(z), quantile(z, probs))
    }))[i, , drop=FALSE]
  } else {
    stop("x must be an 'mcmc.list' or 'rjags'  object.")
  }
  if(!missing(signif)) {
    out[, colnames(out) != 'n.eff'] <-
      signif(out[, colnames(out) != 'n.eff'], signif)
    out
  } else out
}