Creating and running a Markov chain Monte Carlo algorithm in NIMBLE

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

pump example

Let’s use another of the classic WinBUGS examples: the pump model (because it is really easy to show the model and the data).

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)
pumpCode <- nimbleCode({ 
  for (i in 1:N){
      theta[i] ~ dgamma(alpha,beta)
      lambda[i] <- theta[i]*t[i]
      x[i] ~ dpois(lambda[i])
  }
  alpha ~ dexp(1.0)
  beta ~ dgamma(0.1,1.0)
})

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))

pump <- nimbleModel(code = pumpCode, name = 'pump', constants = pumpConsts,
                    data = pumpData, inits = pumpInits)
## defining model...
## building model...
## setting data and initial values...
## checking model...   (use nimbleModel(..., check = FALSE) to skip model check)
## model building finished

Create the MCMC algorithm

pumpMCMC <- buildMCMC(pump)

Run the uncompiled MCMC if you want to

The model and algorithm can be used completely in R. They will be really slow, but running in R allows easy testing and debugging. I’m not going to actually run it for this document because it generates a bunch of warnings that are harmless but annoying. Here is how you would run it for 10 iterations:

pumpMCMC$run(10)

Compile the model and MCMC

It’s more exciting to compile both the model and algorithm and run them.

Cpump <- compileNimble(pump)
CpumpMCMC <- compileNimble(pumpMCMC, project = pump)

Run it and look at the results

CpumpMCMC$run(10000)
## NULL
MCMCsamples <- as.matrix(CpumpMCMC$mvSamples)
plot(MCMCsamples[ , 'alpha'], type = 'l', xlab = 'iteration',  ylab = expression(alpha))

plot of chunk unnamed-chunk-4

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

plot of chunk unnamed-chunk-4