Back to NIMBLE Vignettes

 

NIMBLE’s RW and RW_block samplers contain methods for tracking and accessing the complete history during the course of an MCMC chain of the proposal scale, acceptance rate, and proposal covariance (RW_block sampler only).

This vignette demonstrates how to access these these histories.

 

Specify the two NIMBLE system options:

library(nimble)

nimbleOptions(MCMCsaveHistory = TRUE)
nimbleOptions(buildInterfacesForCompiledNestedNimbleFunctions = TRUE)

 

Build a NIMBLE model and MCMC algorithm.

We’ll used the venerated pump model, and manually configure the MCMC to also include one RW_block sampler.

code <- nimbleCode({ 
    alpha ~ dexp(1)
    beta ~ dgamma(0.1, 1)
    for(i in 1:N) {
        theta[i] ~ dgamma(alpha, beta)
        lambda[i] <- theta[i] * t[i]
        x[i] ~ dpois(lambda[i])
    }
})

N <- 10
constants <- list(N = N, t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5))
data <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22))
inits <- list(alpha = 1, beta = 1, theta = rep(0.1, N))

Rmodel <- nimbleModel(code, constants, data, inits)

conf <- configureMCMC(Rmodel)
## ===== Monitors =====
## thin = 1: alpha, beta
## ===== Samplers =====
## RW sampler (1)
##   - alpha
## conjugate sampler (11)
##   - beta
##   - theta[]  (10 elements)
conf$addSampler(target = c('alpha', 'beta'), type = 'RW_block', silent = TRUE)

Rmcmc <- buildMCMC(conf)

 

Now compile the model and the MCMC, and run the MCMC algorithm for a modest 1000 iterations.

Cmodel <- compileNimble(Rmodel)
Cmcmc <- compileNimble(Rmcmc, project = Rmodel)

set.seed(0)
samples <- runMCMC(Cmcmc, 1000)
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|

 

Now, we’ll access the proposal and acceptance rate histories of the RW and RW_block samplers.

Use conf$printSamplers() to see a printout of all the MCMC samplers.

conf$printSamplers()
## [1]  RW sampler: alpha
## [2]  conjugate_dgamma_dgamma_identity sampler: beta
## [3]  conjugate_dgamma_dpois_multiplicative sampler: theta[1]
## [4]  conjugate_dgamma_dpois_multiplicative sampler: theta[2]
## [5]  conjugate_dgamma_dpois_multiplicative sampler: theta[3]
## [6]  conjugate_dgamma_dpois_multiplicative sampler: theta[4]
## [7]  conjugate_dgamma_dpois_multiplicative sampler: theta[5]
## [8]  conjugate_dgamma_dpois_multiplicative sampler: theta[6]
## [9]  conjugate_dgamma_dpois_multiplicative sampler: theta[7]
## [10] conjugate_dgamma_dpois_multiplicative sampler: theta[8]
## [11] conjugate_dgamma_dpois_multiplicative sampler: theta[9]
## [12] conjugate_dgamma_dpois_multiplicative sampler: theta[10]
## [13] RW_block sampler: alpha, beta

More directly, use the type = 'RW' argument to print only RW and RW_block samplers.

conf$printSamplers(type = 'RW')
## [1]  RW sampler: alpha
## [13] RW_block sampler: alpha, beta

We see the RW sampler on alpha is sampler #1, and the bivariate RW_block operating jointly on alpha and beta is sampler #13.

The MCMC was run for 1000 iterations, and the default value of the adapatation interval for NIMBLE’s MCMC samplers is every 200 iterations. Thus, the samplers underwent five adaptations, and the empirical acceptance rate was recorded five times.

To access the scale and acceptance rate of the RW sampler for each of these adaptive phases:

Cmcmc$samplerFunctions$contentsList[[1]]$getScaleHistory()
## [1] 1.0000000 0.4606066 0.4012447 0.4158587 0.3986886
Cmcmc$samplerFunctions$contentsList[[1]]$getAcceptanceHistory()
## [1] 0.205 0.390 0.455 0.420 0.495

Similarly, the RW_block sampler underwent five adaptations. In addition to a scalar scale and acceptance rate, the RW_block sampler also has (in this case a 2x2) proposal covariance matrix. To access the history of these variables:

Cmcmc$samplerFunctions$contentsList[[13]]$getScaleHistory()
## [1] 1.0000000 0.4530718 0.4657482 0.5636542 0.5072625
Cmcmc$samplerFunctions$contentsList[[13]]$getPropCovHistory()
## , , 1
## 
##           [,1]       [,2]
## [1,] 1.0000000 0.00000000
## [2,] 0.6869277 0.02271129
## [3,] 0.5203157 0.04614878
## [4,] 0.4093119 0.05372164
## [5,] 0.3321212 0.05496567
## 
## , , 2
## 
##            [,1]      [,2]
## [1,] 0.00000000 1.0000000
## [2,] 0.02271129 0.7767706
## [3,] 0.04614878 0.6472710
## [4,] 0.05372164 0.5543444
## [5,] 0.05496567 0.4807200
Cmcmc$samplerFunctions$contentsList[[13]]$getAcceptanceHistory()
## [1] 0.11 0.36 0.43 0.30 0.46

The proposal covariance history is a 5x2x2 array. That is, the initial default proposal covariance matrix was the 2x2 identity matrix:

Cmcmc$samplerFunctions$contentsList[[13]]$getPropCovHistory()[1,,]
##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1

and the proposal covariance after the final adaptation phase was:

Cmcmc$samplerFunctions$contentsList[[13]]$getPropCovHistory()[5,,]
##            [,1]       [,2]
## [1,] 0.33212119 0.05496567
## [2,] 0.05496567 0.48072001

 

We may wish to use these adapated parameters as the starting values of the sampler tuning parameters in a subsequent (perhaps longer) MCMC run.

Let’s extract the final scale value of the RW sampler, and the final proposal covariance of the RW_block sampler:

(RW_scale <- Cmcmc$samplerFunctions$contentsList[[1]]$getScaleHistory()[5])
## [1] 0.3986886
(RW_block_cov <- Cmcmc$samplerFunctions$contentsList[[13]]$getPropCovHistory()[5,,])
##            [,1]       [,2]
## [1,] 0.33212119 0.05496567
## [2,] 0.05496567 0.48072001

 

Finally, we’ll configure and build a new MCMC algorithm, using this scale and proposal covariance as the starting values for the RW and RW_block samplers.

conf <- configureMCMC(Rmodel)
## ===== Monitors =====
## thin = 1: alpha, beta
## ===== Samplers =====
## RW sampler (1)
##   - alpha
## conjugate sampler (11)
##   - beta
##   - theta[]  (10 elements)
conf$removeSampler('alpha', print = TRUE)
## [1]  conjugate_dgamma_dgamma_identity sampler: beta
## [2]  conjugate_dgamma_dpois_multiplicative sampler: theta[1]
## [3]  conjugate_dgamma_dpois_multiplicative sampler: theta[2]
## [4]  conjugate_dgamma_dpois_multiplicative sampler: theta[3]
## [5]  conjugate_dgamma_dpois_multiplicative sampler: theta[4]
## [6]  conjugate_dgamma_dpois_multiplicative sampler: theta[5]
## [7]  conjugate_dgamma_dpois_multiplicative sampler: theta[6]
## [8]  conjugate_dgamma_dpois_multiplicative sampler: theta[7]
## [9]  conjugate_dgamma_dpois_multiplicative sampler: theta[8]
## [10] conjugate_dgamma_dpois_multiplicative sampler: theta[9]
## [11] conjugate_dgamma_dpois_multiplicative sampler: theta[10]
conf$addSampler(target = 'alpha', type = 'RW', scale = RW_scale, print = TRUE)
## [12] RW sampler: alpha,  scale: 0.39868861345930101
conf$addSampler(target = c('alpha', 'beta'), type = 'RW_block',
                propCov = RW_block_cov, silent = TRUE, print = TRUE)
## [13] RW_block sampler: alpha, beta,  propCov: custom array
conf$printSamplers()
## [1]  conjugate_dgamma_dgamma_identity sampler: beta
## [2]  conjugate_dgamma_dpois_multiplicative sampler: theta[1]
## [3]  conjugate_dgamma_dpois_multiplicative sampler: theta[2]
## [4]  conjugate_dgamma_dpois_multiplicative sampler: theta[3]
## [5]  conjugate_dgamma_dpois_multiplicative sampler: theta[4]
## [6]  conjugate_dgamma_dpois_multiplicative sampler: theta[5]
## [7]  conjugate_dgamma_dpois_multiplicative sampler: theta[6]
## [8]  conjugate_dgamma_dpois_multiplicative sampler: theta[7]
## [9]  conjugate_dgamma_dpois_multiplicative sampler: theta[8]
## [10] conjugate_dgamma_dpois_multiplicative sampler: theta[9]
## [11] conjugate_dgamma_dpois_multiplicative sampler: theta[10]
## [12] RW sampler: alpha,  scale: 0.39868861345930101
## [13] RW_block sampler: alpha, beta,  propCov: custom array
Rmcmc <- buildMCMC(conf)

 

Voila!