Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

added maxDimCovHistory control option for RW_block sampler #1329

Merged
merged 3 commits into from
Dec 15, 2023
Merged

added maxDimCovHistory control option for RW_block sampler #1329

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
8df4bd7
added maxDimCovHistory control option for RW_block sampler
danielturek Jun 25, 2023
26b0e03
updated NEWS
danielturek Jun 25, 2023
a36a886
Merge branch 'devel' into RWblock_max_dim
paciorek Dec 15, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
57 changes: 29 additions & 28 deletions packages/nimble/R/MCMC_samplers.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -388,6 +388,7 @@ sampler_RW_block <- nimbleFunction(
scale <- extractControlElement(control, 'scale', 1)
propCov <- extractControlElement(control, 'propCov', 'identity')
tries <- extractControlElement(control, 'tries', 1)
maxDimCovHistory <- extractControlElement(control, 'maxDimCovHistory', 10)
## node list generation
targetAsScalar <- model$expandNodeNames(target, returnScalarComponents = TRUE)
ccList <- mcmc_determineCalcAndCopyNodes(model, target)
Expand All @@ -402,9 +403,9 @@ sampler_RW_block <- nimbleFunction(
timesAccepted <- 0
timesAdapted <- 0
d <- length(targetAsScalar)
scaleHistory <- c(0, 0) ## scaleHistory
acceptanceHistory <- c(0, 0) ## scaleHistory
propCovHistory <- if(d<=10) array(0, c(2,d,d)) else array(0, c(2,2,2)) ## scaleHistory
scaleHistory <- c(0, 0) ## scaleHistory
acceptanceHistory <- c(0, 0) ## scaleHistory
propCovHistory <- if(d <= maxDimCovHistory) array(0, c(2,d,d)) else array(0, c(2,2,2)) ## scaleHistory
saveMCMChistory <- if(nimbleOptions('MCMCsaveHistory')) TRUE else FALSE
if(is.character(propCov) && propCov == 'identity') propCov <- diag(d)
propCovOriginal <- propCov
Expand Down Expand Up @@ -463,7 +464,7 @@ sampler_RW_block <- nimbleFunction(
scaleHistory[timesAdapted] <<- scale ## scaleHistory
setSize(acceptanceHistory, timesAdapted) ## scaleHistory
acceptanceHistory[timesAdapted] <<- acceptanceRate ## scaleHistory
if(d <= 10) {
if(d <= maxDimCovHistory) {
propCovTemp <- propCovHistory ## scaleHistory
setSize(propCovHistory, timesAdapted, d, d) ## scaleHistory
if(timesAdapted > 1) ## scaleHistory
Expand Down Expand Up @@ -509,7 +510,7 @@ sampler_RW_block <- nimbleFunction(
return(acceptanceHistory)
},
getPropCovHistory = function() { ## scaleHistory
if(!saveMCMChistory | d > 10) print("Please set 'nimbleOptions(MCMCsaveHistory = TRUE)' before building the MCMC and note that to reduce memory use we only save the proposal covariance history for parameter vectors of length 10 or less.")
if(!saveMCMChistory | d > maxDimCovHistory) print("Please set 'nimbleOptions(MCMCsaveHistory = TRUE)' before building the MCMC. Note that to reduce memory use, proposal covariance histories are only saved for parameter vectors of length <= 10; this value can be modified using the 'maxDimCovHistory' control list element.")
returnType(double(3))
return(propCovHistory)
},
Expand All @@ -536,7 +537,7 @@ sampler_RW_block <- nimbleFunction(
if(saveMCMChistory) {
scaleHistory <<- c(0, 0) ## scaleHistory
acceptanceHistory <<- c(0, 0)
if(d <= 10)
if(d <= maxDimCovHistory)
propCovHistory <<- nimArray(0, dim = c(2,d,d))
}
my_calcAdaptationFactor$reset()
Expand Down Expand Up @@ -1548,18 +1549,18 @@ sampler_RW_lkj_corr_cholesky <- nimbleFunction(
name = 'sampler_RW_lkj_corr_cholesky',
contains = sampler_BASE,
setup = function(model, mvSaved, target, control) {
## control list extraction
adaptive <- extractControlElement(control, 'adaptive', TRUE)
adaptInterval <- extractControlElement(control, 'adaptInterval', 200)
adaptFactorExponent <- extractControlElement(control, 'adaptFactorExponent', 0.8)
scale <- extractControlElement(control, 'scale', 1)
## node list generation
target <- model$expandNodeNames(target)
targetAsScalar <- model$expandNodeNames(target, returnScalarComponents = TRUE)
calcNodesNoSelf <- model$getDependencies(target, self = FALSE)
##
## numeric value generation
d <- sqrt(length(targetAsScalar))
nTheta <- d*(d-1)/2 # this many unconstrained elements to be sampled
## control list extraction
adaptive <- if(!is.null(control$adaptive)) control$adaptive else TRUE
adaptInterval <- if(!is.null(control$adaptInterval)) control$adaptInterval else 200
adaptFactorExponent <- if(!is.null(control$adaptFactorExponent)) control$adaptFactorExponent else 0.8
scaleOriginal <- if(!is.null(control$scale)) control$scale else 1
if(length(adaptInterval) > 1 || length(adaptFactorExponent) > 1 || length(scaleOriginal) > 1)
stop("RW_lkj_corr_cholesky: 'adaptInterval', 'adaptFactorExponent', and 'scaleOriginal' should be single values.")
if(scaleOriginal < 0)
Expand All @@ -1575,8 +1576,8 @@ sampler_RW_lkj_corr_cholesky <- nimbleFunction(
z <- array(0, c(d, d)) # canonical partial correlations
diag(z) <- 1
partialSums <- array(0, c(d, d)) # 1-x_{13}^2, 1-x_{14}^2, 1-x_{14}^2-x_{24}^2, ...
partialSums[1, ] <- 1
## Temporary vectors for current colum calculations
partialSums[1, ] <- 1
## temporary vectors for current column calculations
partialSumsProp <- numeric(d)
currentValue <- numeric(d)
propValue <- numeric(d)
Expand All @@ -1589,7 +1590,6 @@ sampler_RW_lkj_corr_cholesky <- nimbleFunction(
run = function() {
## calculate transformed values (in unconstrained space) and partial sums in each column
transform(model[[target]]) # compute z and partialSums
##
## Individual univariate RW on the nTheta elements:
## advantage: probably better movement through space than with block update
## (plus note only column values of target matrix are recalculated for a given scalar update in a given column)
Expand Down Expand Up @@ -1688,13 +1688,14 @@ sampler_RW_block_lkj_corr_cholesky <- nimbleFunction(
contains = sampler_BASE,
setup = function(model, mvSaved, target, control) {
## control list extraction
adaptive <- if(!is.null(control$adaptive)) control$adaptive else TRUE
adaptScaleOnly <- if(!is.null(control$adaptScaleOnly)) control$adaptScaleOnly else FALSE
adaptInterval <- if(!is.null(control$adaptInterval)) control$adaptInterval else 200
adaptFactorExponent <- if(!is.null(control$adaptFactorExponent)) control$adaptFactorExponent else 0.8
scale <- if(!is.null(control$scale)) control$scale else 1
propCov <- if(!is.null(control$propCov)) control$propCov else 'identity'
tries <- if(!is.null(control$tries)) control$tries else 1
adaptive <- extractControlElement(control, 'adaptive', TRUE)
adaptScaleOnly <- extractControlElement(control, 'adaptScaleOnly', FALSE)
adaptInterval <- extractControlElement(control, 'adaptInterval', 200)
adaptFactorExponent <- extractControlElement(control, 'adaptFactorExponent', 0.8)
scale <- extractControlElement(control, 'scale', 1)
propCov <- extractControlElement(control, 'propCov', 'identity')
tries <- extractControlElement(control, 'tries', 1)
maxDimCovHistory <- extractControlElement(control, 'maxDimCovHistory', 10)
## node list generation
target <- model$expandNodeNames(target)
targetAsScalar <- model$expandNodeNames(target, returnScalarComponents = TRUE)
Expand All @@ -1714,9 +1715,9 @@ sampler_RW_block_lkj_corr_cholesky <- nimbleFunction(
timesAdapted <- 0
p <- sqrt(length(targetAsScalar))
d <- p*(p-1)/2 # this many unconstrained elements to be sampled
scaleHistory <- c(0, 0) ## scaleHistory
acceptanceHistory <- c(0, 0) ## scaleHistory
propCovHistory <- if(d<=10) array(0, c(2,d,d)) else array(0, c(2,2,2)) ## scaleHistory
scaleHistory <- c(0, 0) ## scaleHistory
acceptanceHistory <- c(0, 0) ## scaleHistory
propCovHistory <- if(d <= maxDimCovHistory) array(0, c(2,d,d)) else array(0, c(2,2,2)) ## scaleHistory
saveMCMChistory <- if(nimbleOptions('MCMCsaveHistory')) TRUE else FALSE
if(is.character(propCov) && propCov == 'identity') propCov <- diag(d)
propCovOriginal <- propCov
Expand Down Expand Up @@ -1839,7 +1840,7 @@ sampler_RW_block_lkj_corr_cholesky <- nimbleFunction(
scaleHistory[timesAdapted] <<- scale ## scaleHistory
setSize(acceptanceHistory, timesAdapted) ## scaleHistory
acceptanceHistory[timesAdapted] <<- acceptanceRate ## scaleHistory
if(d <= 10) {
if(d <= maxDimCovHistory) {
propCovTemp <- propCovHistory ## scaleHistory
setSize(propCovHistory, timesAdapted, d, d) ## scaleHistory
if(timesAdapted > 1) ## scaleHistory
Expand Down Expand Up @@ -1874,7 +1875,7 @@ sampler_RW_block_lkj_corr_cholesky <- nimbleFunction(
return(acceptanceHistory)
},
getPropCovHistory = function() { ## scaleHistory
if(!saveMCMChistory | d > 10) print("Please set 'nimbleOptions(MCMCsaveHistory = TRUE)' before building the MCMC and note that to reduce memory use we only save the proposal covariance history for parameter vectors of length 10 or less.")
if(!saveMCMChistory | d > maxDimCovHistory) print("Please set 'nimbleOptions(MCMCsaveHistory = TRUE)' before building the MCMC. Note that to reduce memory use, proposal covariance histories are only saved for parameter vectors of length <= 10; this value can be modified using the 'maxDimCovHistory' control list element.")
returnType(double(3))
return(propCovHistory)
},
Expand All @@ -1901,7 +1902,7 @@ sampler_RW_block_lkj_corr_cholesky <- nimbleFunction(
if(saveMCMChistory) {
scaleHistory <<- c(0, 0) ## scaleHistory
acceptanceHistory <<- c(0, 0)
if(d <= 10)
if(d <= maxDimCovHistory)
propCovHistory <<- nimArray(0, dim = c(2,d,d))
}
my_calcAdaptationFactor$reset()
Expand Down
3 changes: 3 additions & 0 deletions packages/nimble/inst/NEWS.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,9 @@

- Change argument of `besselK` in manual table to be `x` not `k`.

- Add new control list option, `maxDimCovHistory` to `RW_block` sampler
specify maximum dimension for saving proposal covariance history.

## DEVELOPER LEVEL CHANGES

- Make change to `nimble-package` documentation to use `"_PACKAGE"`
Expand Down
Loading