Back to NIMBLE Vignettes
This will demonstrate how to toggle MCMC samplers on/off. It’s a fairly custom approach, not really intended for mainstream use, but the techniques could be refined. In any case, it should get you going.
I would be happy to be kept in the loop regarding this research, and possibly continue to support the NIMBLE/MCMC side as an active contributor.
Load NIMBLE library.
library(nimble)You’ll need this sampler definition.
It will generate a warning when you define it. You can thank Chris for that ;)
sampler_toggle <- nimbleFunction(
contains = sampler_BASE,
setup = function(model, mvSaved, target, control) {
type <- control$type
nested_sampler_name <- paste0('sampler_', type)
control_new <- nimbleOptions('MCMCcontrolDefaultList')
control_new[[names(control)]] <- control
nested_sampler_list <- nimbleFunctionList(sampler_BASE)
nested_sampler_list[[1]] <- do.call(nested_sampler_name, list(model, mvSaved, target, control_new))
toggle <- 1
},
run = function() {
if(toggle == 1)
nested_sampler_list[[1]]$run()
},
methods = list(
reset = function()
nested_sampler_list[[1]]$reset()
)
)## Warning in nf_checkDSLcode(code): For this nimbleFunction to compile, these
## functions must be defined as nimbleFunctions or nimbleFunction methods:
## reset.This next section only needs to be done once. There are several important specifics here, and generally all steps must be done, in order.
The model is a simplification of yours, but it should work fine adding additional layers above y.censored and X.
code <- nimbleCode({
##q[1:N,1:N] ~ dwish(R = aq[1:N,1:N], df = bq)
##Q[1:N,1:N] <- inverse(q[1:N,1:N])
##X.mod[1:N] ~ dmnorm(muf[1:N],prec = pf[1:N,1:N])
## add process error
X[1:N] ~ dmnorm(X.mod[1:N],prec = q[1:N,1:N])
## Analysis
y.censored[1:N] ~ dmnorm(X[1:N], prec = r[1:N,1:N])
for(i in 1:N){
y.ind[i] ~ dconstraint(y.censored[i] > 0)
}
})
N <- 5 ## works with other values of N
constants <- list(N=N, X.mod=rep(10,N), q=diag(N), r=diag(N))
## the values of y.ind and y.censored are just placeholders...
## BUT, important, y.censored *must* be specified as data at this point.
## it's a long story why, has to do with the initializeModel routine
## at the beginning of MCMC execution
data <- list(y.ind=rep(1,N), y.censored=rep(10,N))
inits <- list(X=rep(10,N))
Rmodel <- nimbleModel(code, constants, data, inits)## Adding X.mod,q,r as data for building model.conf <- configureMCMC(Rmodel, print=TRUE)## [1] conjugate_dmnorm_dmnorm sampler: X[1:5]## important!
## this is needed for correct indexing later
samplerNumberOffset <- length(conf$getSamplers())
for(i in 1:N) {
node <- paste0('y.censored[',i,']')
conf$addSampler(node, 'toggle', control=list(type='RW'))
## could instead use slice samplers, or any combination thereof, e.g.:
##conf$addSampler(node, 'toggle', control=list(type='slice'))
}
conf$printSamplers()## [1] conjugate_dmnorm_dmnorm sampler: X[1:5]
## [2] toggle sampler: y.censored[1], type: RW
## [3] toggle sampler: y.censored[2], type: RW
## [4] toggle sampler: y.censored[3], type: RW
## [5] toggle sampler: y.censored[4], type: RW
## [6] toggle sampler: y.censored[5], type: RWconf$printMonitors()## thin = 1: X## can monitor y.censored, if you wish, to verify correct behaviour
conf$addMonitors('y.censored')## thin = 1: X, y.censoredRmcmc <- buildMCMC(conf)
Cmodel <- compileNimble(Rmodel)
Cmcmc <- compileNimble(Rmcmc, project = Rmodel)The steps hereafter can be repeated with new datasets (as demonstrated)
Mix of data and non-data.
y.ind <- c(1, 1, 1, 0, 0)
## important!
## change: rather than providing NA for the non-data values (those to be sampled),
## you'll have to provide some values here.
## that's because we effectively disabled the model initialization routine earlier
y.censored <- c(9, 9, 11, 20, 20)
Cmodel$y.ind <- y.ind
Cmodel$y.censored <- y.censored
for(i in 1:N) {
## ironically, here we have to "toggle" the value of y.ind[i]
## this specifies that when y.ind[i] = 1,
## indicator variable is set to 0, which specifies *not* to sample
valueInCompiledNimbleFunction(Cmcmc$samplerFunctions[[samplerNumberOffset+i]], 'toggle', 1-y.ind[i])
}
niter <- 1000
set.seed(0)
samples <- runMCMC(Cmcmc, niter, progressBar=FALSE)## running chain 1...head(samples)## X[1] X[2] X[3] X[4] X[5] y.censored[1]
## [1,] 10.393044 9.269318 11.44031 15.89974 15.29320 9
## [2,] 9.291601 9.495922 12.20035 14.76997 13.97073 9
## [3,] 9.288423 9.209018 10.67835 13.02551 14.69906 9
## [4,] 9.766859 9.594283 11.06865 12.99705 14.63496 9
## [5,] 9.161897 9.795529 11.46073 12.96475 13.81760 9
## [6,] 9.206555 9.221195 10.45806 13.78836 13.39459 9
## y.censored[2] y.censored[3] y.censored[4] y.censored[5]
## [1,] 9 11 18.46005 19.07143
## [2,] 9 11 17.31239 18.78197
## [3,] 9 11 16.07485 18.55770
## [4,] 9 11 16.07485 17.08695
## [5,] 9 11 16.02105 15.70989
## [6,] 9 11 15.85653 15.45653tail(samples)## X[1] X[2] X[3] X[4] X[5] y.censored[1]
## [995,] 9.140770 10.012700 10.816557 8.565837 9.749762 9
## [996,] 10.680715 8.858806 9.891825 9.818678 10.221955 9
## [997,] 9.342837 10.067078 10.584596 8.954470 11.456985 9
## [998,] 10.109610 9.518324 9.783697 8.952032 9.452207 9
## [999,] 10.294098 11.135964 9.750794 10.530146 10.559196 9
## [1000,] 10.150436 10.261331 10.312753 9.534064 10.413439 9
## y.censored[2] y.censored[3] y.censored[4] y.censored[5]
## [995,] 9 11 7.474586 10.79963
## [996,] 9 11 7.474586 10.79963
## [997,] 9 11 7.474586 10.80651
## [998,] 9 11 7.474586 10.80651
## [999,] 9 11 7.474586 10.80651
## [1000,] 9 11 7.658921 10.80651
Everything is data.
y.ind <- c(1, 1, 1, 1, 1) ## everything is data
y.censored <- c(9, 9, 11, 11, 12)
Cmodel$y.ind <- y.ind
Cmodel$y.censored <- y.censored
for(i in 1:N) {
valueInCompiledNimbleFunction(Cmcmc$samplerFunctions[[samplerNumberOffset+i]], 'toggle', 1-y.ind[i])
}
niter <- 1000
set.seed(0)
samples <- runMCMC(Cmcmc, niter, progressBar=FALSE)## running chain 1...head(samples)## X[1] X[2] X[3] X[4] X[5] y.censored[1]
## [1,] 10.393044 9.269318 11.440310 11.399743 11.29320 9
## [2,] 8.411091 8.843404 10.291601 10.495922 12.70035 9
## [3,] 10.039942 8.935015 9.688484 10.295320 10.78842 9
## [4,] 9.209018 9.678349 9.869317 10.808075 10.12493 9
## [5,] 9.341419 9.766859 10.594283 11.068648 10.95962 9
## [6,] 9.856105 10.267755 10.011422 9.591651 11.03304 9
## y.censored[2] y.censored[3] y.censored[4] y.censored[5]
## [1,] 9 11 11 12
## [2,] 9 11 11 12
## [3,] 9 11 11 12
## [4,] 9 11 11 12
## [5,] 9 11 11 12
## [6,] 9 11 11 12tail(samples)## X[1] X[2] X[3] X[4] X[5] y.censored[1]
## [995,] 8.999090 9.868799 9.927915 9.973051 10.10205 9
## [996,] 8.657415 9.527744 11.268144 10.624727 10.57735 9
## [997,] 10.641860 8.836730 10.513638 10.610404 11.63634 9
## [998,] 9.608817 8.616277 9.820056 9.819228 10.14049 9
## [999,] 9.858187 10.723524 10.635293 8.844648 10.57828 9
## [1000,] 10.817599 11.029022 11.331527 10.937754 11.00598 9
## y.censored[2] y.censored[3] y.censored[4] y.censored[5]
## [995,] 9 11 11 12
## [996,] 9 11 11 12
## [997,] 9 11 11 12
## [998,] 9 11 11 12
## [999,] 9 11 11 12
## [1000,] 9 11 11 12
Nothing is data.
y.ind <- c(0, 0, 0, 0, 0) ## nothing is data
## again, you'll have to provide actual numeric (initial) values
## for the nodes that will be sampled:
y.censored <- c(20, 20, 20, 20, 20)
Cmodel$y.ind <- y.ind
Cmodel$y.censored <- y.censored
for(i in 1:N) {
valueInCompiledNimbleFunction(Cmcmc$samplerFunctions[[samplerNumberOffset+i]], 'toggle', 1-y.ind[i])
}
niter <- 1000
set.seed(0)
samples <- runMCMC(Cmcmc, niter, progressBar=FALSE)## running chain 1...head(samples)## X[1] X[2] X[3] X[4] X[5] y.censored[1]
## [1,] 15.89304 14.76932 15.94031 15.89974 15.29320 18.46005
## [2,] 15.29902 14.81138 14.41336 13.43109 15.79545 18.41512
## [3,] 14.47442 14.62190 15.42129 14.33797 15.24397 18.41512
## [4,] 13.90116 13.33300 15.12746 15.19287 15.31799 17.98560
## [5,] 13.16791 13.03876 13.50780 15.19614 15.20477 17.75827
## [6,] 14.05601 14.22941 14.18549 12.97637 13.69547 17.75827
## y.censored[2] y.censored[3] y.censored[4] y.censored[5]
## [1,] 19.07143 19.70528 19.99423 20.00000
## [2,] 19.05524 19.70528 18.75669 19.77573
## [3,] 17.58449 19.22713 18.75669 19.23284
## [4,] 17.58449 19.22713 18.75669 19.23284
## [5,] 17.58449 18.61510 19.09781 19.17797
## [6,] 17.54525 18.61510 18.28285 19.42023tail(samples)## X[1] X[2] X[3] X[4] X[5] y.censored[1]
## [995,] 8.522043 10.684843 12.17742 9.298193 10.349728 8.733287
## [996,] 9.100242 11.307015 10.51391 11.102941 10.690697 6.358070
## [997,] 7.150281 9.572054 10.95131 11.090109 9.906660 6.358070
## [998,] 8.959874 11.368493 12.95964 11.229054 11.374903 7.040653
## [999,] 8.923126 10.969560 11.09850 9.453420 11.220409 8.252536
## [1000,] 8.701234 12.654686 11.55302 8.821897 9.974707 8.252536
## y.censored[2] y.censored[3] y.censored[4] y.censored[5]
## [995,] 11.32094 12.44488 9.590339 10.17209
## [996,] 11.11355 12.44488 10.691195 10.17209
## [997,] 11.11355 14.02627 13.026282 11.82663
## [998,] 12.04485 14.02627 11.401560 13.67533
## [999,] 12.04485 12.52238 7.680039 13.05367
## [1000,] 12.04485 12.52238 7.680039 12.97148Send along any questions or comments.
-Daniel