Back to NIMBLE Vignettes

 

 

 

This document defines a sampler RW_record that will record the scaleHistory, and acceptanceRateHistory. It then shows how to access those using a simple model.

 

Note that you need to specify the NIMBLE system option:

nimbleOptions(buildInterfacesForCompiledNestedNimbleFunctions = TRUE)

 

library(nimble)

## need to specify this option:
nimbleOptions(buildInterfacesForCompiledNestedNimbleFunctions = TRUE)

sampler_RW_record <- nimbleFunction(
    contains = sampler_BASE,
    setup = function(model, mvSaved, target, control) {
        ## control list extraction
        logScale      <- control$log
        reflective    <- control$reflective
        adaptive      <- control$adaptive
        adaptInterval <- control$adaptInterval
        scale         <- control$scale
        ## node list generation
        targetAsScalar <- model$expandNodeNames(target, returnScalarComponents = TRUE)
        calcNodes  <- model$getDependencies(target)
        ## numeric value generation
        scaleOriginal <- scale
        timesRan      <- 0
        timesAccepted <- 0
        timesAdapted  <- 0
        scaleHistory  <- c(0, 0)   ## scaleHistory
        acceptanceRateHistory <- c(0, 0)   ## scaleHistory
        optimalAR     <- 0.44
        gamma1        <- 0
        ## checks
        if(length(targetAsScalar) > 1)   stop('cannot use RW sampler on more than one target; try RW_block sampler')
        if(model$isDiscrete(target))     stop('cannot use RW sampler on discrete-valued target; try slice sampler')
        if(logScale & reflective)        stop('cannot use reflective RW sampler on a log scale (i.e. with options log=TRUE and reflective=TRUE')
    },
    run = function() {
        currentValue <- model[[target]]
        propLogScale <- 0
        if(logScale) { propLogScale <- rnorm(1, mean = 0, sd = scale)
                       propValue <- currentValue * exp(propLogScale)
                   } else         propValue <- rnorm(1, mean = currentValue,  sd = scale)
        if(reflective) {
            lower <- model$getBound(target, 'lower')
            upper <- model$getBound(target, 'upper')
            while(propValue < lower | propValue > upper) {
                if(propValue < lower) propValue <- 2*lower - propValue
                if(propValue > upper) propValue <- 2*upper - propValue
            }
        }
        model[[target]] <<- propValue
        logMHR <- calculateDiff(model, calcNodes) + propLogScale
        jump <- decide(logMHR)
        if(jump) nimCopy(from = model, to = mvSaved, row = 1, nodes = calcNodes, logProb = TRUE)
        else     nimCopy(from = mvSaved, to = model, row = 1, nodes = calcNodes, logProb = TRUE)
        if(adaptive)     adaptiveProcedure(jump)
    },
    methods = list(
        adaptiveProcedure = function(jump = logical()) {
            timesRan <<- timesRan + 1
            if(jump)     timesAccepted <<- timesAccepted + 1
            if(timesRan %% adaptInterval == 0) {
                acceptanceRate <- timesAccepted / timesRan
                timesAdapted <<- timesAdapted + 1
                setSize(scaleHistory,          timesAdapted)         ## scaleHistory
                setSize(acceptanceRateHistory, timesAdapted)         ## scaleHistory
                scaleHistory[timesAdapted]          <<- scale        ## scaleHistory
                acceptanceRateHistory[timesAdapted] <<- acceptanceRate        ## scaleHistory
                gamma1 <<- 1/((timesAdapted + 3)^0.8)
                gamma2 <- 10 * gamma1
                adaptFactor <- exp(gamma2 * (acceptanceRate - optimalAR))
                scale <<- scale * adaptFactor
                timesRan <<- 0
                timesAccepted <<- 0
            }
        },
        getScaleHistory = function()          { returnType(double(1)); return(scaleHistory) },          ## scaleHistory
        getAcceptanceRateHistory = function() { returnType(double(1)); return(acceptanceRateHistory) },          ## scaleHistory
        reset = function() {
            scale <<- scaleOriginal
            timesRan      <<- 0
            timesAccepted <<- 0
            timesAdapted  <<- 0
            scaleHistory           <<- scaleHistory * 0    ## scaleHistory
            acceptanceRateHistory  <<- acceptanceRateHistory * 0    ## scaleHistory
            gamma1 <<- 0
        }
    )
)
## Warning in nf_checkDSLcode(code): For this nimbleFunction to compile, these
## functions must be defined as nimbleFunctions or nimbleFunction methods:
## adaptiveProcedure.
code <- nimbleCode({
    a ~ dnorm(0, sd=1)
    b ~ dnorm(0, sd=10)
})
constants <- list()
data <- list()
inits <- list(a=0, b=0)

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

conf <- configureMCMC(Rmodel, nodes=NULL)
conf$addSampler('a', 'RW_record')
conf$addSampler('b', 'RW_record')
conf$printSamplers()
## [1] RW_record sampler: a
## [2] RW_record sampler: b
Rmcmc <- buildMCMC(conf)

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

Cmcmc$run(10000)
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## NULL
## scale history and acceptance rate history for sampler #1
Cmcmc$samplerFunctions$contentsList[[1]]$getScaleHistory()
##  [1] 1.000000 2.477281 2.344266 2.488291 2.385554 2.340781 2.242022
##  [8] 2.565347 2.436823 2.591756 2.558670 2.253188 2.115570 2.047621
## [15] 2.179028 2.370397 2.260615 2.169885 2.188964 2.235612 2.395768
## [22] 2.376993 2.469238 2.460144 2.348273 2.228928 2.305577 2.406344
## [29] 2.460947 2.554985 2.531721 2.539269 2.629470 2.614556 2.651183
## [36] 2.572932 2.579805 2.600116 2.606789 2.574226 2.567883 2.488284
## [43] 2.512076 2.454035 2.414881 2.398568 2.335448 2.361123 2.422901
## [50] 2.453905
Cmcmc$samplerFunctions$contentsList[[1]]$getAcceptanceRateHistory()
##  [1] 0.715 0.420 0.465 0.420 0.430 0.415 0.525 0.405 0.485 0.430 0.335
## [12] 0.385 0.410 0.500 0.525 0.390 0.395 0.450 0.465 0.525 0.430 0.490
## [23] 0.435 0.375 0.365 0.490 0.505 0.475 0.500 0.425 0.445 0.500 0.430
## [34] 0.465 0.385 0.445 0.455 0.445 0.415 0.435 0.375 0.460 0.390 0.405
## [45] 0.425 0.380 0.465 0.500 0.470 0.355
## scale history and acceptance rate history for sampler #2
Cmcmc$samplerFunctions$contentsList[[2]]$getScaleHistory()
##  [1]  1.000000  6.238983 16.616984 22.123029 24.069585 23.842645 23.034436
##  [8] 22.852622 23.533776 24.521014 24.838099 23.952019 24.789696 24.924943
## [15] 24.796082 25.167178 24.929618 24.258040 24.578686 24.578686 24.880627
## [22] 24.204863 22.515742 22.766356 24.022318 23.363263 23.521783 23.444510
## [29] 23.294693 23.294693 23.652543 23.442238 22.969742 22.073150 22.382369
## [36] 23.000353 23.433876 24.307381 24.557859 24.681648 25.173633 25.480363
## [43] 25.000001 24.422387 24.032723 24.087082 24.140671 25.054973 25.380635
## [50] 24.743339
Cmcmc$samplerFunctions$contentsList[[2]]$getAcceptanceRateHistory()
##  [1] 0.995 0.795 0.560 0.480 0.435 0.420 0.435 0.460 0.470 0.450 0.410
## [12] 0.470 0.445 0.435 0.455 0.430 0.410 0.455 0.440 0.455 0.405 0.345
## [23] 0.455 0.515 0.400 0.450 0.435 0.430 0.440 0.465 0.425 0.405 0.370
## [34] 0.465 0.490 0.475 0.510 0.460 0.450 0.480 0.465 0.400 0.390 0.405
## [45] 0.445 0.445 0.525 0.470 0.380 0.385