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)
bonesModel <- readBUGSmodel('bones', dir = getBUGSexampleDir('bones'))
## defining model...
## Detected grade as data within 'constants'.
## Adding grade as data for building model.
## building model...
## setting data and initial values...
## checking model...   (use nimbleModel(..., check = FALSE) to skip model check)
## model building finished

Make an MCMC configuration object

An MCMC configuration holds the information on which samplers are included in the MCMC, which nodes they operate on, and any parameters needed for them. It also has a set of nodes to include (monitor) in the MCMC output. Actually it allows two different sets of nodes to be monitored, each with its own thinning interval. 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)
bonesMCMCconfiguration$getSamplers()
## [1]  RW sampler: theta[1],  adaptive: TRUE,  adaptInterval: 200,  scale: 1
## [2]  RW sampler: theta[2],  adaptive: TRUE,  adaptInterval: 200,  scale: 1
## [3]  RW sampler: theta[3],  adaptive: TRUE,  adaptInterval: 200,  scale: 1
## [4]  RW sampler: theta[4],  adaptive: TRUE,  adaptInterval: 200,  scale: 1
## [5]  RW sampler: theta[5],  adaptive: TRUE,  adaptInterval: 200,  scale: 1
## [6]  RW sampler: theta[6],  adaptive: TRUE,  adaptInterval: 200,  scale: 1
## [7]  RW sampler: theta[7],  adaptive: TRUE,  adaptInterval: 200,  scale: 1
## [8]  RW sampler: theta[8],  adaptive: TRUE,  adaptInterval: 200,  scale: 1
## [9]  RW sampler: theta[9],  adaptive: TRUE,  adaptInterval: 200,  scale: 1
## [10] RW sampler: theta[10],  adaptive: TRUE,  adaptInterval: 200,  scale: 1
## [11] RW sampler: theta[11],  adaptive: TRUE,  adaptInterval: 200,  scale: 1
## [12] RW sampler: theta[12],  adaptive: TRUE,  adaptInterval: 200,  scale: 1
## [13] RW sampler: theta[13],  adaptive: TRUE,  adaptInterval: 200,  scale: 1
## [14] end sampler: grade[4, 13]
## [15] end sampler: grade[6, 18]
## [16] end sampler: grade[7, 10]
## [17] end sampler: grade[7, 11]
## [18] end sampler: grade[7, 18]
## [19] end sampler: grade[8, 10]
## [20] end sampler: grade[8, 11]
## [21] end sampler: grade[10, 18]
## [22] end sampler: grade[11, 3]
## [23] end sampler: grade[11, 7]
## [24] end sampler: grade[11, 11]
## [25] end sampler: grade[11, 12]
## [26] end sampler: grade[11, 15]
## [27] end sampler: grade[11, 16]
## [28] end sampler: grade[11, 22]
## [29] end sampler: grade[11, 24]
## [30] end sampler: grade[12, 24]
## [31] end sampler: grade[13, 11]
## [32] end sampler: grade[13, 23]
## [33] end sampler: grade[13, 25]

Now we can see that theta[1]-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 essentially posterior predictive nodes – so they have been assigned end samplers.

Note that if we had called buildMCMC(bonesModel), it would have made the default MCMC configuration and then built the MCMC algorithm in one step.

Customize by replacing the univariate samplers for theta dimensions with a block sampler.

Let’s say we want to compare the efficiency of the univariate samplers to 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')
## [21] RW_block sampler: theta[1:13],  adaptive: TRUE,  adaptScaleOnly: FALSE,  adaptInterval: 200,  scale: 1,  propCov: identity

Build the customized MCMC

bonesMCMC <- buildMCMC(bonesMCMCconfiguration)

Compile the model and MCMC

Cbones <- compileNimble(bonesModel, bonesMCMC)

Run the MCMC

Cbones$bonesMCMC$run(10000)
## NULL
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 = '.')

plot of chunk unnamed-chunk-7

Of course we haven’t actually compared the efficiencies of the default to the customized MCMC. We’ll do that kind of thing in other examples.

Writing your own samplers and running different MCMCs side by side

You can learn to write your own samplers and include them in an MCMC in Chapter 7 of the User Manual. You can also learn to run multiple MCMCs side by side, including from WinBUGS, JAGS, and Stan (if you provide the necessary files) using NIMBLE’s MCMCsuite function.