diff --git a/nimbleHMC/DESCRIPTION b/nimbleHMC/DESCRIPTION index 1896948..d9e2663 100644 --- a/nimbleHMC/DESCRIPTION +++ b/nimbleHMC/DESCRIPTION @@ -1,13 +1,13 @@ Package: nimbleHMC Title: Hamiltonian Monte Carlo and Other Gradient-Based MCMC Sampling Algorithms for 'nimble' -Version: 0.2.3 -Date: 2024-12-18 +Version: 0.2.4 +Date: 2025-12-15 Authors@R: c(person("Daniel", "Turek", role = c("aut", "cre"), email = "danielturek@gmail.com"), person("Perry", "de Valpine", role = "aut"), person("Christopher", "Paciorek", role = "aut")) Maintainer: Daniel Turek Description: Provides gradient-based MCMC sampling algorithms for use with the MCMC engine provided by the 'nimble' package. This includes two versions of Hamiltonian Monte Carlo (HMC) No-U-Turn (NUTS) sampling, and (under development) Langevin samplers. The `NUTS_classic` sampler implements the original HMC-NUTS algorithm as described in Hoffman and Gelman (2014) . The `NUTS` sampler is a modern version of HMC-NUTS sampling matching the HMC sampler available in version 2.32.2 of Stan (Stan Development Team, 2023). In addition, convenience functions are provided for generating and modifying MCMC configuration objects which employ HMC sampling. Functionality of the 'nimbleHMC' package is described further in Turek, et al (2024) . -Depends: R (>= 3.5.0), nimble (>= 1.0.0) +Depends: R (>= 3.5.0), nimble (>= 1.4.0) Imports: methods Suggests: testthat License: BSD_3_clause + file LICENSE | GPL (>= 2) diff --git a/nimbleHMC/R/HMC_configuration.R b/nimbleHMC/R/HMC_configuration.R index 29836f3..52e26d1 100644 --- a/nimbleHMC/R/HMC_configuration.R +++ b/nimbleHMC/R/HMC_configuration.R @@ -256,9 +256,11 @@ buildHMC <- function(model, nodes = character(), type = 'NUTS', control = list() #' @param samplesAsCodaMCMC Logical argument. If \code{TRUE}, then a \code{coda} \code{mcmc} object is returned instead of an R matrix of samples, or when \code{nchains > 1} a \code{coda} \code{mcmc.list} object is returned containing \code{nchains} \code{mcmc} objects. This argument is only used when \code{samples} is \code{TRUE}. Default value is \code{FALSE}. See details. #' #' @param summary Logical argument. When \code{TRUE}, summary statistics for the posterior samples of each parameter are also returned, for each MCMC chain. This may be returned in addition to the posterior samples themselves. Default value is \code{FALSE}. See details. -#'z +#' #' @param WAIC Logical argument. When \code{TRUE}, the WAIC (Watanabe, 2010) of the model is calculated and returned. If multiple chains are run, then a single WAIC value is calculated using the posterior samples from all chains. Default value is \code{FALSE}. Note that the version of WAIC used is the default WAIC conditional on random effects/latent states and without any grouping of data nodes. See \code{help(waic)} for more details. #' +#' @param userEnv Environment in which if-then-else statements in model code will be evaluated if needed information not found in \code{constants}; intended primarily for internal use only. +#' #' @return A list is returned with named elements depending on the arguments, unless only one among samples, summary, and WAIC are requested, in which case only that element is returned. These elements may include \code{samples}, \code{summary}, and \code{WAIC}. When \code{nchains = 1}, posterior samples are returned as a single matrix, and summary statistics as a single matrix. When \code{nchains > 1}, posterior samples are returned as a list of matrices, one matrix for each chain, and summary statistics are returned as a list containing \code{nchains+1} matrices: one matrix corresponding to each chain, and the final element providing a summary of all chains, combined. If \code{samplesAsCodaMCMC} is \code{TRUE}, then posterior samples are provided as \code{coda} \code{mcmc} and \code{mcmc.list} objects. When \code{WAIC} is \code{TRUE}, a WAIC summary object is returned. #' #' @details @@ -320,11 +322,12 @@ nimbleHMC <- function(code, samples = TRUE, samplesAsCodaMCMC = FALSE, summary = FALSE, - WAIC = FALSE) { + WAIC = FALSE, + userEnv = parent.frame()) { if(missing(code) && missing(model)) stop('must provide either code or model argument') if(!samples && !summary && !WAIC) stop('no output specified, use samples = TRUE, summary = TRUE, or WAIC = TRUE') if(!missing(code) && inherits(code, 'modelBaseClass')) model <- code ## let's handle it, if model object is provided as un-named first argument - Rmodel <- mcmc_createModelObject(model, inits, nchains, setSeed, code, constants, data, dimensions, check, buildDerivs = TRUE) + Rmodel <- mcmc_createModelObject(model, inits, nchains, setSeed, code, constants, data, dimensions, check, buildDerivs = TRUE, userEnv = userEnv) Rmcmc <- buildHMC(Rmodel, type = type, monitors = monitors, thin = thin, enableWAIC = WAIC, print = FALSE) compiledList <- compileNimble(Rmodel, Rmcmc) ## only one compileNimble() call Cmcmc <- compiledList$Rmcmc diff --git a/nimbleHMC/R/HMC_samplers.R b/nimbleHMC/R/HMC_samplers.R index e720564..9da05c5 100644 --- a/nimbleHMC/R/HMC_samplers.R +++ b/nimbleHMC/R/HMC_samplers.R @@ -21,7 +21,6 @@ ## #' } ## #' ## #' @import nimble -## #' @import methods ## #' ## #' @aliases Langevin langevin ## #' diff --git a/nimbleHMC/man/buildHMC.Rd b/nimbleHMC/man/buildHMC.Rd index cac1f83..cc37c60 100644 --- a/nimbleHMC/man/buildHMC.Rd +++ b/nimbleHMC/man/buildHMC.Rd @@ -14,7 +14,7 @@ buildHMC( ) } \arguments{ -\item{model}{A nimble model, as returned by `nimbleModel`} +\item{model}{A nimble model, as returned by `nimbleModel`. Alternatively, an MCMC configuration object, as returned by either `configureHMC` or `configureHMC`. See details.} \item{nodes}{A character vector of stochastic node names to be sampled. If an empty character vector is provided (the default), then all stochastic non-data nodes will be sampled. An HMC sampler will be applied to all continuous-valued non-data nodes, and nimble's default sampler will be assigned for all discrete-valued nodes.} @@ -35,7 +35,9 @@ Build an MCMC algorithm which applies HMC sampling to continuous-valued dimensio \details{ This is the most direct way to create an MCMC algorithm using HMC sampling in nimble. This will create a compilable, executable MCMC algorithm, with HMC sampling assigned to all continuous-valued model dimensions, and nimble's default sampler assigned to all discrete-valued dimensions. The `nodes` argument can be used to control which model nodes are assigned samplers. Use this if you don't otherwise need to modify the MCMC configuration. -Either the `NUTS_classic` or the `NUTS` samplin can be applied. Both implement variants of No-U-Turn HMC sampling, however the `NUTS` sampler uses more modern adapatation techniques. See `help(NUTS)` or `help(NUTS_classic)` for details. +Either the `NUTS_classic` or the `NUTS` sampling can be applied, which is controled by the `type` argument. Both implement variants of No-U-Turn HMC sampling, however the `NUTS` sampler uses more modern adapatation techniques. See `help(NUTS)` or `help(NUTS_classic)` for details. + +Note that when an MCMC configuration object is provided as the `model` argument, then an executable MCMC algorithm is generated using the MCMC configuration that was provided - regardless of whether or not it specifies any HMC samplers. } \examples{ code <- nimbleCode({ diff --git a/nimbleHMC/man/nimbleHMC.Rd b/nimbleHMC/man/nimbleHMC.Rd index 0224d9f..e01511c 100644 --- a/nimbleHMC/man/nimbleHMC.Rd +++ b/nimbleHMC/man/nimbleHMC.Rd @@ -23,7 +23,8 @@ nimbleHMC( samples = TRUE, samplesAsCodaMCMC = FALSE, summary = FALSE, - WAIC = FALSE + WAIC = FALSE, + userEnv = parent.frame() ) } \arguments{ @@ -61,10 +62,11 @@ nimbleHMC( \item{samplesAsCodaMCMC}{Logical argument. If \code{TRUE}, then a \code{coda} \code{mcmc} object is returned instead of an R matrix of samples, or when \code{nchains > 1} a \code{coda} \code{mcmc.list} object is returned containing \code{nchains} \code{mcmc} objects. This argument is only used when \code{samples} is \code{TRUE}. Default value is \code{FALSE}. See details.} -\item{summary}{Logical argument. When \code{TRUE}, summary statistics for the posterior samples of each parameter are also returned, for each MCMC chain. This may be returned in addition to the posterior samples themselves. Default value is \code{FALSE}. See details. -z} +\item{summary}{Logical argument. When \code{TRUE}, summary statistics for the posterior samples of each parameter are also returned, for each MCMC chain. This may be returned in addition to the posterior samples themselves. Default value is \code{FALSE}. See details.} \item{WAIC}{Logical argument. When \code{TRUE}, the WAIC (Watanabe, 2010) of the model is calculated and returned. If multiple chains are run, then a single WAIC value is calculated using the posterior samples from all chains. Default value is \code{FALSE}. Note that the version of WAIC used is the default WAIC conditional on random effects/latent states and without any grouping of data nodes. See \code{help(waic)} for more details.} + +\item{userEnv}{Environment in which if-then-else statements in model code will be evaluated if needed information not found in \code{constants}; intended primarily for internal use only.} } \value{ A list is returned with named elements depending on the arguments, unless only one among samples, summary, and WAIC are requested, in which case only that element is returned. These elements may include \code{samples}, \code{summary}, and \code{WAIC}. When \code{nchains = 1}, posterior samples are returned as a single matrix, and summary statistics as a single matrix. When \code{nchains > 1}, posterior samples are returned as a list of matrices, one matrix for each chain, and summary statistics are returned as a list containing \code{nchains+1} matrices: one matrix corresponding to each chain, and the final element providing a summary of all chains, combined. If \code{samplesAsCodaMCMC} is \code{TRUE}, then posterior samples are provided as \code{coda} \code{mcmc} and \code{mcmc.list} objects. When \code{WAIC} is \code{TRUE}, a WAIC summary object is returned. diff --git a/nimbleHMC/tests/testthat/test-HMC.R b/nimbleHMC/tests/testthat/test-HMC.R index 1f57843..687dc8f 100644 --- a/nimbleHMC/tests/testthat/test-HMC.R +++ b/nimbleHMC/tests/testthat/test-HMC.R @@ -316,7 +316,7 @@ test_that('HMC on conjugate Wishart', { OmegaSimTrueSDs <- apply(wishRV, c(1,2), sd) ## expect_equal(as.numeric(apply(samples, 2, mean)), as.numeric(OmegaTrueMean), tol = 0.1) - expect_equal(as.numeric(apply(samples, 2, sd)), as.numeric(OmegaSimTrueSDs), tol = 0.02) + expect_equal(as.numeric(apply(samples, 2, sd)), as.numeric(OmegaSimTrueSDs), tol = 0.1) }) test_that('HMC on LKJ', { @@ -564,8 +564,8 @@ test_that('HMC runs with various non-differentiable constructs', { }) test_that('HMC results for CAR match non-HMC', { - set.seed(1) - code <- nimbleCode({ + set.seed(1) + code <- nimbleCode({ S[1:N] ~ dcar_normal(adj[1:L], weights[1:L], num[1:N], tau) tau ~ dunif(0, 5) for(i in 1:N) @@ -575,7 +575,7 @@ test_that('HMC results for CAR match non-HMC', { Y[i] ~ dpois(lambda[i]) } }) - + ## constants <- list(N = 6, num = c(1,2,2,2,2,1), adj = c(2, 1,3, 2,4, 3,5, 4,6, 5), @@ -583,24 +583,23 @@ test_that('HMC results for CAR match non-HMC', { L = 10) data <- list(Y = c(1,0,2,1,4,3)) inits <- list(tau = 1, S = c(0,0,0,0,0,0)) - + ## Rmodel <- nimbleModel(code, constants, data, inits, buildDerivs = TRUE) conf <- configureMCMC(Rmodel, monitors = c('tau','S')) Rmcmc <- buildMCMC(conf) Cmodel <- compileNimble(Rmodel) Cmcmc <- compileNimble(Rmcmc, project = Rmodel) out <- runMCMC(Cmcmc, niter = 505000, nburnin = 5000, thin = 500) - + ## Rmodel <- nimbleModel(code, constants, data, inits, buildDerivs = TRUE) conf <- configureHMC(Rmodel, monitors = c('tau','S')) Rmcmc <- buildMCMC(conf) Cmodel <- compileNimble(Rmodel) Cmcmc <- compileNimble(Rmcmc, project = Rmodel) outHMC <- runMCMC(Cmcmc, niter = 22000, nburnin=2000, thin=20) - + ## expect_equal(apply(out[,1:6],2,mean), apply(outHMC[,1:6],2,mean), tolerance = .06) expect_equal(mean(out[,7]),mean(outHMC[,7]), tolerance = .15) - expect_equal(apply(out,2,quantile,c(.1,.9)), apply(outHMC,2,quantile,c(.1,.9)), tolerance = 0.15) }) @@ -624,7 +623,7 @@ test_that('HMC results for mixture model match non-HMC', { mu <- c(0,2,4) data <- list(y=sample(c(rnorm(50,mu[1],.35), rnorm(250,mu[2],.35), rnorm(200,mu[3],.35)), n, replace=FALSE)) inits <- list(k=sample(1:K,n,replace=T),mu=c(-1,2,6),mu0=1) - + ## m <- nimbleModel(code, constants = constants, data = data, inits = inits, buildDerivs = TRUE) conf <- configureMCMC(m, monitors=c('mu','mu0','p')) @@ -632,27 +631,25 @@ test_that('HMC results for mixture model match non-HMC', { cm <- compileNimble(m) cmcmc <- compileNimble(mcmc) out <- runMCMC(cmcmc, niter=50000, nburnin=10000, thin=40) - + ## m <- nimbleModel(code, constants = constants, data = data, inits = inits, buildDerivs = TRUE) conf <- configureMCMC(m, nodes=c('k'), monitors=c('mu','mu0','p')) - conf$addSampler(c('mu0','mu','p'),'NUTS') + conf$addSampler(c('mu0','mu','p'), 'NUTS') mcmc <- buildMCMC(conf) cm <- compileNimble(m) cmcmc <- compileNimble(mcmc) outHMC <- runMCMC(cmcmc, niter=22000, nburnin=2000, thin=20) - - ## Deal with label switching. + ## deal with label switching sorter <- function(row) { ord <- order(row[1:3]) return(c(row[1:3][ord], row[4], row[5:7][ord])) } out <- t(apply(out, 1, sorter)) outHMC <- t(apply(outHMC, 1, sorter)) - - expect_equal(apply(out,2,mean), apply(outHMC,2,mean), tolerance = .1) + ## + expect_equal(apply(out,2,mean), apply(outHMC,2,mean), tolerance = 0.1) expect_equal(apply(out,2,quantile,c(.1,.9)), apply(outHMC,2,quantile,c(.1,.9)), tolerance = 0.1) - })