jags_sample: Run JAGS models with each MCMC chain computed in serial or...
In michaelmalick/r-jagstools: jagstools

Description Usage Arguments Details Value Author(s) See Also Examples

View source: R/jags_sample.R

jags_sample is a wrapper for functions in the rjags package that allows for running JAGS from R with each MCMC chain computed serially or in parallel using the parallel package.

jags_sample(data, inits, file, variable.names, n.chains = 1, n.adapt = 100,
  burnin = 100, n.iter = 1000, thin = 1, load.modules = NULL,
  method = "serial", progress.bar = "text", parallel = list(n.clusters =
  n.chains, RNGseed = 123, type = "PSOCK", verbose = FALSE))

jags_parallel(data, inits, file, variable.names, n.chains = 1,
  n.adapt = 100, burnin = 100, n.iter = 1000, thin = 1,
  load.modules = NULL, progress.bar = "text", parallel = list(n.clusters =
  n.chains, RNGseed = 123, type = "PSOCK", verbose = FALSE))

jags_serial(data, inits, file, variable.names, n.chains = 1, n.adapt = 100,
  burnin = 100, n.iter = 1000, thin = 1, load.modules = NULL,
  progress.bar = "text")

`data`	a list containing the data for the model; see `?rjags::jags.model` for details.
`inits`	a list or function giving the initial values for each variable; see `?rjags::jags.model` for details.
`file`	path to the file containing the JAGS model; see `?rjags::jags.model` for details.
`variable.names`	character vector of names of variables to monitor in JAGS; see `?rjags::coda.samples` for details.
`n.chains`	number of MCMC chains to run; see `?rjags::jags.model` for details.
`n.adapt`	number of iterations for adaption; see `?rjags::jags.model` for details.
`burnin`	number of iterations to use for burnin of the MCMC chains. The burnin iterations are not saved.
`n.iter`	number of samples for each Markov chain; see `?rjags::coda.samples` for details.
`thin`	thinning interval; see `?coda.samples` for details.
`load.modules`	JAGS modules to load before running model; see `?rjags::load.modules` for details.
`method`	string identifying method for computation; either 'serial' or 'parallel'; defaults to 'serial'.
`progress.bar`	type of progress bar to use; see `rjags::update.jags` for details
`parallel`	a list of arguments for parallel computing: `n.clusters` number of cluster to use for computation; currently needs to be set equal to `n.chains`. `RNGseed` seed for the RNG stream sent to the clusters; setting this to a positive integer allows the chains to be reproducible. `type` type of cluster to create; either "PSOCK" or "FORK"; see `?parallel::makeCluster`. `verbose` logical, should output from each cluster be sent to terminal. Can be useful for debugging.

If method = 'serial' the function runs JAGS using the rjags::coda.samples function. If method = 'parallel' the function uses the parallel::parLapply function to send each chain to different clusters. Unfortunately, when the parallel type = 'PSOCK' no progress bars are available. If parallel type = 'FORK' and verbose = TRUE, progress bars should be printed to the terminal among other potentially useful output.

Because the funtion returns an mcmc.list object, the model (as defined by rjags::jags.model) is not retained and therefore, updating or extending the MCMC runs is currently unavailable. Adding this functionality would require saving the model object, which would make the return not an mcmc.list object and would likely require the creation of a new class and additional functions to deal with this new class.

An mcmc.list

Michael Malick

jags.model coda.samples mcmc.list load.module unload.module

# ----------------------------
# Simulate data
# ----------------------------

library(rjags)

## True parameter values
alpha <- 3
beta  <- 0.5
sigma <- 5.0
N     <- 100

## Simulate data
set.seed(129)
x       <- runif(N, 5, 95)
epsilon <- rnorm(N, 0, sigma)
y       <- alpha + beta*x + epsilon

# Return data in a list
dat <- list(N = N, x = x, y = y)



# ----------------------------
# Setup JAGS model and write
# to file
# ----------------------------
model.string <-
    'model {
        for (i in 1:N){
            y[i] ~ dnorm(mu.y[i], tau)
            mu.y[i] <- alpha + beta * x[i]
        }
        # Priors
        alpha ~ dnorm(0, .0001)
        beta  ~ dnorm(0, .0001)
        tau   <- pow(sigma, -2)
        sigma ~ dunif(0, 100)
    }'
writeLines(model.string, con = "example_jags.bug")


# ----------------------------
# Setup monitors and inits
# ----------------------------
## Parameters to monitor
params <- c("alpha", "beta", "sigma", "deviance")

## Setup function for initial values
inits <- function() {
    list(.RNG.name  = "lecuyer::RngStream",
         .RNG.seed  = runif(1, 0, 2^31),
         alpha = rnorm(1), beta = rnorm(1), sigma = runif(1))
}


# ----------------------------
# Sample all chains in serial
# ----------------------------
fit.s <- jags_sample(data = dat,
                     inits = inits,
                     file = "example_jags.bug",
                     variable.names = params,
                     n.chains = 2,
                     n.adapt = 100,
                     burnin = 100,
                     n.iter = 1000,
                     thin = 1,
                     load.modules = c("dic", "lecuyer"),
                     method = "serial",
                     progress.bar = "text")

# ----------------------------
# Sample all chains in parallel
# ----------------------------

## type = "PSOCK"
fit.p <- jags_sample(data = dat,
                     inits = inits,
                     file = "example_jags.bug",
                     variable.names = params,
                     n.chains = 2,
                     n.adapt = 100,
                     burnin = 100,
                     n.iter = 1000,
                     thin = 1,
                     load.modules = c("dic", "lecuyer"),
                     method = "parallel",
                     progress.bar = "text",
                     parallel = list(n.clusters = 2,
                                     RNGseed = 123,
                                     type = "PSOCK",
                                     verbose = TRUE))

## type = "FORK"
fit.p <- jags_sample(data = dat,
                     inits = inits,
                     file = "example_jags.bug",
                     variable.names = params,
                     n.chains = 2,
                     n.adapt = 100,
                     burnin = 100,
                     n.iter = 1000,
                     thin = 1,
                     load.modules = c("dic", "lecuyer"),
                     method = "parallel",
                     progress.bar = "text",
                     parallel = list(n.clusters = 2,
                                     RNGseed = 123,
                                     type = "FORK",
                                     verbose = TRUE))