Skip to content

safe changing of MCMC thinning intervals when reset=FALSE #1609

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

Merged
danielturek merged 6 commits into from
Jan 23, 2026
Merged

safe changing of MCMC thinning intervals when reset=FALSE #1609

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
7cf4264
continuing to change book-keeping of MCMC iterations
danielturek Jan 14, 2026
0b250cc
wrong indexing for thin1
danielturek Jan 15, 2026
f8369fd
changed hard stop to a warning
danielturek Jan 15, 2026
48bd8df
corrected calculation for mvSamples resize when reset = FALSE
danielturek Jan 15, 2026
ffd58f6
updates to logic for changing thin when reset=FALSE
danielturek Jan 16, 2026
654bb88
Add tests for changing thin when reset is FALSE.
paciorek Jan 22, 2026
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
52 changes: 35 additions & 17 deletions packages/nimble/R/MCMC_build.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -245,6 +245,8 @@ buildMCMC <- nimbleFunction(
nburnin_extraWAIC <- 0
}
iterationTotal <- 0
iterUntilSaveVec <- c(0, 0)
warnUnfinishedThinningInterval <- getNimbleOption('MCMCwarnUnfinishedThinningInterval')
firstRun <- TRUE
setupOutputs(derivedTypes)
},
Expand All @@ -268,15 +270,13 @@ buildMCMC <- nimbleFunction(
if(nburnin > niter) stop('cannot specify nburnin > niter')
if(firstRun) reset <- TRUE ## compulsory reset on first run of MCMC
firstRun <<- FALSE
mvSamples_copyRow <- 0
mvSamples2_copyRow <- 0
if(reset) {
if(initializeModel) my_initializeModel$run()
thinToUseVec <<- thinFromConfVec
if(thin != -1) thinToUseVec[1] <<- thin
if(thin2 != -1) thinToUseVec[2] <<- thin2
for(iThin in 1:2) {
if(thinToUseVec[iThin] < 1) stop('cannot use thin < 1')
if(thinToUseVec[iThin] != floor(thinToUseVec[iThin])) stop('cannot use non-integer thin')
}
for(i in seq_along(derivedFunctions)) {
if(derivedIntervals[i] == 0) {
derivedIntervals[i] <<- thinToUseVec[1]
Expand All @@ -289,25 +289,40 @@ buildMCMC <- nimbleFunction(
for(i in seq_along(derivedFunctions)) derivedFunctions[[i]]$before_chain(niter-nburnin, nburnin, thinToUseVec, chain)
samplerTimes <<- numeric(length(samplerFunctions) + 1) ## default inititialization to zero
iterationTotal <<- 0
mvSamples_copyRow <- 0
mvSamples2_copyRow <- 0
iterUntilSaveVec <<- thinToUseVec
} else {
if(nburnin != 0) stop('cannot specify nburnin when using reset = FALSE.')
if((thin != -1) & (thin != thinToUseVec[1])) stop('cannot alter the value of thin, when using reset = FALSE.')
if((thin2 != -1) & (thin2 != thinToUseVec[2])) stop('cannot alter the value of thin2, when using reset = FALSE.')
if(nburnin != 0) stop('cannot specify nburnin when using reset = FALSE.')
## now allowing changing of thin and thin2, when reset = FALSE; Jan 2026.
## when reset = FALSE and a new value of thin (or thin2) is provided,
## issue a warning if the previous MCMC run didn't exactly finish a thinning/save cycle.
## this is checked using iterUntilSaveVec[1] == thinToUseVec[1] or iterUntilSaveVec[2] == thinToUseVec[2],
## which means saving on the previous thinning interval just took place.
if((thin != -1) & (thin != thinToUseVec[1])) {
if((iterUntilSaveVec[1] != thinToUseVec[1]) & warnUnfinishedThinningInterval)
cat('Warning: value of thin is being changed, in the midst of a thinning interval.\n')
thinToUseVec[1] <<- thin
iterUntilSaveVec[1] <<- thin
}
if((thin2 != -1) & (thin2 != thinToUseVec[2])) {
if((iterUntilSaveVec[2] != thinToUseVec[2]) & warnUnfinishedThinningInterval)
cat('Warning: value of thin2 is being changed, in the midst of a thinning interval.\n')
thinToUseVec[2] <<- thin2
iterUntilSaveVec[2] <<- thin2
}
if(dim(samplerTimes)[1] != length(samplerFunctions) + 1) samplerTimes <<- numeric(length(samplerFunctions) + 1) ## first run: default inititialization to zero
if (resetMV) {
mvSamples_copyRow <- 0
mvSamples2_copyRow <- 0
} else {
if(!resetMV) {
mvSamples_copyRow <- getsize(mvSamples)
mvSamples2_copyRow <- getsize(mvSamples2)
}
}
for(iThin in 1:2) {
if(thinToUseVec[iThin] < 1) stop('cannot use thin < 1')
if(thinToUseVec[iThin] != floor(thinToUseVec[iThin])) stop('cannot use non-integer thin')
}
nimCopy(from = model, to = mvSaved, row = 1, logProb = TRUE)
if(onlineWAIC & resetWAIC) waicFun[[1]]$reset()
resize(mvSamples, floor((iterationTotal+niter-nburnin) / thinToUseVec[1]))
resize(mvSamples2, floor((iterationTotal+niter-nburnin) / thinToUseVec[2]))
resize(mvSamples, mvSamples_copyRow + 1 + floor((niter-nburnin-iterUntilSaveVec[1]) / thinToUseVec[1]))
resize(mvSamples2, mvSamples2_copyRow + 1 + floor((niter-nburnin-iterUntilSaveVec[2]) / thinToUseVec[2]))
## reinstate samplerExecutionOrder as a runtime argument, once we support non-scalar default values for runtime arguments:
##if(dim(samplerExecutionOrder)[1] > 0 & samplerExecutionOrder[1] == -1) { ## runtime argument samplerExecutionOrder was not provided
## lengthSamplerExecutionOrderFromConf <- dim(samplerExecutionOrderFromConfPlusTwoZeros)[1] - 2
Expand Down Expand Up @@ -339,15 +354,18 @@ buildMCMC <- nimbleFunction(
}
}
if(iter > nburnin) {
iterUntilSaveVec <<- iterUntilSaveVec - 1
## save samples
iterPostBurnin <- iterationTotal - nburnin
if(iterPostBurnin %% thinToUseVec[1] == 0) {
if(iterUntilSaveVec[1] == 0) {
mvSamples_copyRow <- mvSamples_copyRow + 1
nimCopy(from = model, to = mvSamples, row = mvSamples_copyRow, nodes = monitors)
iterUntilSaveVec[1] <<- thinToUseVec[1]
}
if(iterPostBurnin %% thinToUseVec[2] == 0) {
if(iterUntilSaveVec[2] == 0) {
mvSamples2_copyRow <- mvSamples2_copyRow + 1
nimCopy(from = model, to = mvSamples2, row = mvSamples2_copyRow, nodes = monitors2)
iterUntilSaveVec[2] <<- thinToUseVec[2]
}
## save WAIC
if(enableWAIC & onlineWAIC & iter > nburnin + nburnin_extraWAIC) {
Expand Down
1 change: 1 addition & 0 deletions packages/nimble/R/options.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -213,6 +213,7 @@ nimOptimMethod("bobyqa",
MCMCassignSamplersToPosteriorPredictiveNodes = TRUE, ## whether any samplers are assigned (by default) to PP nodes
MCMCusePosteriorPredictiveSampler = TRUE, ## for PP nodes being sampled, use post_pred (or otherwise RW, etc)
MCMCwarnUnsampledStochasticNodes = TRUE,
MCMCwarnUnfinishedThinningInterval = TRUE,
MCMCRJcheckHyperparam = TRUE,
MCMCenableWAIC = FALSE,
MCMCuseBarkerAsDefaultMV = FALSE,
Expand Down
100 changes: 100 additions & 0 deletions packages/nimble/tests/testthat/test-mcmc.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3058,6 +3058,106 @@ test_that('asymptotic correct results from conjugate gamma - CAR_normal sampler'
expect_true(all(abs(summary_RW - summary_conj) < 0.04))
})

test_that("correct handling of changing `thin` value with `reset=FALSE`", {
code <- nimbleCode({
y ~ dnorm(mu, sd = sigma)
mu ~ dnorm(0,1)
sigma ~ dunif(0,5)
})

m <- nimbleModel(code, data = list(y=0), inits = list(mu = 0, sigma = 1))
cm <- compileNimble(m)

mcmc <- buildMCMC(m, thin = 1)
cmcmc <- compileNimble(mcmc, project = m)

set.seed(1)
cmcmc$run(niter = 20)
gold <- as.matrix(cmcmc$mvSamples)

## thin=1 -> thin=2
set.seed(1)
cm$mu <- 0
cm$sigma <- 1
cmcmc$run(niter = 8)
expect_silent(cmcmc$run(niter = 12, thin = 2, reset = FALSE))
expect_identical(as.matrix(cmcmc$mvSamples),
gold[c(1:8,seq(10,20,by=2)),])

## thin=2 -> thin=1
set.seed(1)
cm$mu <- 0
cm$sigma <- 1
cmcmc$run(niter = 8, thin = 2)
expect_silent(cmcmc$run(niter = 12, thin = 1, reset = FALSE))
expect_identical(as.matrix(cmcmc$mvSamples),
gold[c(seq(2,8,by=2),9:20),])

set.seed(1)
cm$mu <- 0
cm$sigma <- 1
cmcmc$run(niter = 9, thin = 2)
expect_output(cmcmc$run(niter = 11, thin = 1, reset = FALSE),
"thin is being changed")
expect_identical(as.matrix(cmcmc$mvSamples),
gold[c(seq(2,8,by=2),10:20),])

## thin=5 -> thin=1
set.seed(1)
cm$mu <- 0
cm$sigma <- 1
cmcmc$run(niter = 10, thin = 5)
expect_silent(cmcmc$run(niter = 10, thin = 1, reset = FALSE))
expect_identical(as.matrix(cmcmc$mvSamples),
gold[c(5,10,11:20),])

set.seed(1)
cm$mu <- 0
cm$sigma <- 1
cmcmc$run(niter = 8, thin = 5)
expect_output(cmcmc$run(niter = 12, thin = 1, reset = FALSE),
"thin is being changed")
expect_identical(as.matrix(cmcmc$mvSamples),
gold[c(5,9:20),])

## thin=4 -> thin=3
set.seed(1)
cm$mu <- 0
cm$sigma <- 1
cmcmc$run(niter = 8, thin = 4)
expect_silent(cmcmc$run(niter = 12, thin = 3, reset = FALSE))
expect_identical(as.matrix(cmcmc$mvSamples),
gold[c(4,8,seq(11,20,by=3)),])

set.seed(1)
cm$mu <- 0
cm$sigma <- 1
cmcmc$run(niter = 9, thin = 4)
expect_output(cmcmc$run(niter = 11, thin = 3, reset = FALSE),
"thin is being changed")
expect_identical(as.matrix(cmcmc$mvSamples),
gold[c(4,8,seq(12,18,by=3)),])

## thin=4 -> thin=5
set.seed(1)
cm$mu <- 0
cm$sigma <- 1
cmcmc$run(niter = 8, thin = 4)
expect_silent(cmcmc$run(niter = 12, thin = 5, reset = FALSE))
expect_identical(as.matrix(cmcmc$mvSamples),
gold[c(4,8,13,18),])

set.seed(1)
cm$mu <- 0
cm$sigma <- 1
cmcmc$run(niter = 9, thin = 4)
expect_output(cmcmc$run(niter = 11, thin = 5, reset = FALSE),
"thin is being changed")
expect_identical(as.matrix(cmcmc$mvSamples),
gold[c(4,8,14,19),])

})

sink(NULL)

if(!generatingGoldFile) {
Expand Down
Loading