Changing the samplers used in an MCMC

This example shows how you can control which samplers are included in an MCMC.

Bones example

Let’s use another of the classic WinBUGS examples: the bones example.

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

The BUGS code looks like this:

{
   for (i in 1:nChild) {
      theta[i] ~ dnorm(0.0, 0.001);

      for (j in 1:nInd) { 
         # Cumulative probability of > grade k given theta
         for (k in 1:(ncat[j]-1)) {
            logit(Q[i,j,k]) <- delta[j]*(theta[i] - gamma[j,k]);
         }
         Q[i,j,ncat[j]] <- 0;
      }

      for (j in 1:nInd) {
         # Probability of observing grade k given theta
         p[i,j,1] <- 1 - Q[i,j,1];
         for (k in 2:ncat[j]) {
            p[i,j,k] <- Q[i,j,(k-1)] - Q[i,j,k];
         }
         grade[i,j] ~ dcat(p[i,j,1:ncat[j]]);
      }
   }
}   

We will load it this way to avoid showing a bunch of data here:

library(nimble, warn.conflicts = FALSE)

## nimble version 1.3.0 is loaded.
## For more information on NIMBLE and a User Manual,
## please visit https://R-nimble.org.
## 
## Note for advanced users who have written their own MCMC samplers:
##   As of version 0.13.0, NIMBLE's protocol for handling posterior
##   predictive nodes has changed in a way that could affect user-defined
##   samplers in some situations. Please see Section 15.5.1 of the User Manual.

bonesModel <- readBUGSmodel('bones', dir = getBUGSexampleDir('bones'))

## Defining model

##   [Note] 'nGrade' is provided in 'constants' but not used in the model code and is being ignored.

##   [Note] Using 'grade' (given within 'constants') as data.

## Building model

## Setting data and initial values

## Running calculate on model
##   [Note] Any error reports that follow may simply reflect missing values in model variables.

## Checking model sizes and dimensions

##   [Note] This model is not fully initialized. This is not an error.
##          To see which variables are not initialized, use model$initializeInfo().
##          For more information on model initialization, see help(modelInitialization).

Make an MCMC configuration object

An MCMC configuration holds the information on which samplers are included in the MCMC, which nodes the samplers operate on, and any parameters they need. We can modify the MCMC configuration before we build the MCMC algorithm from it.

Here is how to make the configuration and look at the default samplers:

bonesMCMCconfiguration <- configureMCMC(bonesModel)

## ===== Monitors =====
## thin = 1: theta
## ===== Samplers =====
## RW sampler (13)
##   - theta[]  (13 elements)
## posterior_predictive sampler (20)
##   - grade[]  (20 elements)

bonesMCMCconfiguration$printSamplers()

## [1]  RW sampler: theta[1]
## [2]  RW sampler: theta[2]
## [3]  RW sampler: theta[3]
## [4]  RW sampler: theta[4]
## [5]  RW sampler: theta[5]
## [6]  RW sampler: theta[6]
## [7]  RW sampler: theta[7]
## [8]  RW sampler: theta[8]
## [9]  RW sampler: theta[9]
## [10] RW sampler: theta[10]
## [11] RW sampler: theta[11]
## [12] RW sampler: theta[12]
## [13] RW sampler: theta[13]
## [14] posterior_predictive sampler: grade[4, 13]
## [15] posterior_predictive sampler: grade[6, 18]
## [16] posterior_predictive sampler: grade[7, 10]
## [17] posterior_predictive sampler: grade[7, 11]
## [18] posterior_predictive sampler: grade[7, 18]
## [19] posterior_predictive sampler: grade[8, 10]
## [20] posterior_predictive sampler: grade[8, 11]
## [21] posterior_predictive sampler: grade[10, 18]
## [22] posterior_predictive sampler: grade[11, 3]
## [23] posterior_predictive sampler: grade[11, 7]
## [24] posterior_predictive sampler: grade[11, 11]
## [25] posterior_predictive sampler: grade[11, 12]
## [26] posterior_predictive sampler: grade[11, 15]
## [27] posterior_predictive sampler: grade[11, 16]
## [28] posterior_predictive sampler: grade[11, 22]
## [29] posterior_predictive sampler: grade[11, 24]
## [30] posterior_predictive sampler: grade[12, 24]
## [31] posterior_predictive sampler: grade[13, 11]
## [32] posterior_predictive sampler: grade[13, 23]
## [33] posterior_predictive sampler: grade[13, 25]

Now we can see that theta[1] through theta[13] have each been assigned adaptive random walk Metropolis-Hastings samplers. A smattering of entries in the grade matrix are missing. Those have no dependencies – they are posterior predictive nodes – so they have been assigned samplers that simply sample from the predictive distribution. As of nimble version 0.13.0, this sampling does not affect the sampling of the model parameters, because the samplers for the parameters do not condition on the value(s) of posterior predictive nodes. This new behavior generally improves MCMC mixing.

The MCMC configuration also has a set of nodes to include (or monitor) in the MCMC output, that can be accessed as in the following.

bonesMCMCconfiguration$monitors

## [1] "theta"

In addition it allows a second sets of nodes to be monitored with its own thinning interval, which is empty by default.

bonesMCMCconfiguration$monitors2

## character(0)

Note that if we call buildMCMC(bonesModel) at this stage, it would make the default MCMC configuration and then build the MCMC algorithm in one step.

Customize the MCMC

Let’s say we want to replace the univariate samplers with a block sampler. We can remove the univariate samplers and insert a block sampler like this:

bonesMCMCconfiguration$removeSamplers('theta', print = FALSE)
bonesMCMCconfiguration$addSampler(target = 'theta[1:13]', type = 'RW_block')

##   [Note] Assigning an RW_block sampler to nodes with very different scales can result in low MCMC efficiency.  If all nodes assigned to RW_block are not on a similar scale, we recommend providing an informed value for the "propCov" control list argument, or using the "barker" sampler instead.

Build the customized MCMC

bonesMCMC <- buildMCMC(bonesMCMCconfiguration)

##   [Note] Reordering posterior_predictive samplers to execute last

Compile the model and MCMC

Cbones <- compileNimble(bonesModel, bonesMCMC)

## Compiling
##   [Note] This may take a minute.
##   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.

Run the MCMC

Cbones$bonesMCMC$run(10000)

## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|

MCMCsamples <- as.matrix(Cbones$bonesMCMC$mvSamples)

Look at samples from theta[1:4] because that fits conveniently on one plot:

pairs(MCMCsamples[,1:4], pch = '.')

Writing your own samplers

You can find more on the MCMC configuration and its customization in Chapter 7.2 of the User Manual, while you can learn how to write your own samplers and include them in an MCMC in Chapter 15 of the User Manual.