Creating and running a Markov chain Monte Carlo (MCMC) algorithm in NIMBLE

This example shows how to quickly build and run an MCMC algorithm in NIMBLE with default sampler choices.

Pump example

Let’s use another of the classic WinBUGS examples: the pump model.

A description can be found here. The original example can be found in our GitHub repository here.

Create the model

We could load the model using readBUGSmodel(), but instead we’ll show it fully here:

library(nimble, warn.conflicts = FALSE)
## nimble version 0.9.0 is loaded.
## For more information on NIMBLE and a User Manual,
## please visit http://R-nimble.org.
pumpCode <- nimbleCode({ 
# Define relationships between nodes
  for (i in 1:N){
      theta[i] ~ dgamma(alpha,beta)
      lambda[i] <- theta[i]*t[i]
      x[i] ~ dpois(lambda[i])
  }
# Set priors
  alpha ~ dexp(1.0)
  beta ~ dgamma(0.1,1.0)
})

# Create some contrants, data, and initial values to pass to the model builder
pumpConsts <- list(N = 10,
                   t = c(94.3, 15.7, 62.9, 126, 5.24,
                       31.4, 1.05, 1.05, 2.1, 10.5))

pumpData <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22))

pumpInits <- list(alpha = 1, beta = 1,
                  theta = rep(0.1, pumpConsts$N))

Three ways to run an MCMC

NIMBLE provides three different workflows for running an MCMC. These give you different levels of control over how the algorithm runs. From least to most control, these are: the nimbleMCMC() function, the runMCMC() function, and the mcmc$run() MCMC object function.

In general, we recommend the runMCMC() path, although nimbleMCMC() is useful for very straightforward cases.

nimbleMCMC()

The nimbleMCMC() function handles everything from model building to compilation to running the MCMC. It takes model code and model inputs (constants, data, and initial values) and returns MCMC samples. The nimbleMCMC() function is easy and quick to use but provides little control over the MCMC algorithm. Also, since the objects (the uncompiled and compiled model and MCMC) are used only in the function, they aren’t saved and can’t be manipulated for debugging.

nimbleMCMC_samples <- nimbleMCMC(code = pumpCode, 
                                 constants = pumpConsts, 
                                 data = pumpData, 
                                 inits = pumpInits,
                                 nburnin = 1000, niter = 10000)
## defining model...
## building model...
## setting data and initial values...
## running calculate on model (any error reports that follow may simply reflect missing values in model variables) ... 
## checking model sizes and dimensions...
## checking model calculations...
## model building finished.
## compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.
## compilation finished.
## running chain 1...
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
head(nimbleMCMC_samples)
##          alpha      beta
## [1,] 0.5309487 0.8042517
## [2,] 0.5309487 0.3869900
## [3,] 0.5309487 1.1717346
## [4,] 0.5309487 0.7636882
## [5,] 0.6799529 0.5995017
## [6,] 0.6799529 1.0561407

We can look at the results:

plot(nimbleMCMC_samples[ , 'alpha'], type = 'l', xlab = 'iteration',  
     ylab = expression(alpha))

plot(nimbleMCMC_samples[ , 'beta'], type = 'l', xlab = 'iteration',
     ylab = expression(beta))

runMCMC()

The runMCMC() function provides more control and a closer look at what NIMBLE is doing. runMCMC() takes a NIMBLE MCMC object, so we need to back up and build these things.

First, build the NIMBLE model.

pumpModel <- nimbleModel(code = pumpCode, name = 'pump', constants = pumpConsts,
                    data = pumpData, inits = pumpInits)
## defining model...
## building model...
## setting data and initial values...
## running calculate on model (any error reports that follow may simply reflect missing values in model variables) ... 
## checking model sizes and dimensions...
## model building finished.

Next, build an MCMC object for this model using buildMCMC().

pumpMCMC <- buildMCMC(pumpModel)

We can use runMCMC() on this uncompiled model if we want to, but we’ll get a warning reminding us that compiled MCMC will be much faster. Uncompiled MCMC will be really slow, but running in R allows easy testing and debugging. Here is how you would run it for 5 iterations:

runMCMC(pumpMCMC, niter = 5)
## Warning: running an uncompiled MCMC algorithm, use compileNimble() for faster execution.
## running chain 1...
##         alpha     beta
## [1,] 1.000000 2.155030
## [2,] 1.000000 1.428661
## [3,] 1.000000 3.548065
## [4,] 1.000000 2.545029
## [5,] 1.232827 1.778516

To compile the model and MCMC object, use compileNimble(). Note that the model object must be compiled first for the MCMC to compile.

CpumpModel <- compileNimble(pumpModel)
## compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.
## compilation finished.
CpumpMCMC <- compileNimble(pumpMCMC, project = pumpModel)
## compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.
## compilation finished.

Now we can use runMCMC() on the compiled model object.

runMCMC_samples <- runMCMC(CpumpMCMC, nburnin = 1000, niter = 10000)
## running chain 1...
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|

We can look at the results:

plot(runMCMC_samples[ , 'alpha'], type = 'l', xlab = 'iteration',  ylab = expression(alpha))