diff --git a/packages/nimble/R/MCMC_samplers.R b/packages/nimble/R/MCMC_samplers.R index 92b8e9d31..0e5628b74 100644 --- a/packages/nimble/R/MCMC_samplers.R +++ b/packages/nimble/R/MCMC_samplers.R @@ -32,45 +32,39 @@ sampler_prior_samples <- nimbleFunction( contains = sampler_BASE, setup = function(model, mvSaved, target, control) { ## control list extraction - samples <- extractControlElement(control, 'samples', error = ' [Error] prior_samples sampler missing required control argument: samples') - randomDraws <- extractControlElement(control, 'randomDraws', FALSE) ## default is taking sequential draws + samples <- extractControlElement(control, 'samples', error = ' [Error] prior_samples sampler missing required control argument: samples') ## node list generation targetExpanded <- model$expandNodeNames(target) targetAsScalar <- model$expandNodeNames(targetExpanded, returnScalarComponents = TRUE) ccList <- mcmc_determineCalcAndCopyNodes(model, targetExpanded) - calcNodes <- ccList$calcNodes; calcNodesNoSelf <- ccList$calcNodesNoSelf; copyNodesDeterm <- ccList$copyNodesDeterm; copyNodesStoch <- ccList$copyNodesStoch + calcNodesNoSelf <- ccList$calcNodesNoSelf; copyNodesDeterm <- ccList$copyNodesDeterm; copyNodesStoch <- ccList$copyNodesStoch ## numeric value generation k <- length(targetAsScalar) if(is.null(dim(samples))) samples <- matrix(samples, ncol = 1) ## make vectors into 1-column array nSamples <- dim(samples)[1] - ind <- 0 ## checks - if(length(dim(samples)) != 2) stop('prior_samples sampler \'samples\' control argument must be a 2-dimensional array, but value provided was a ', length(dim(samples)), '-dimensional array') + if(length(dim(samples)) != 2) stop('prior_samples sampler \'samples\' control argument must be a vector or matrix, but value provided was a ', length(dim(samples)), '-dimensional array') if(!(storage.mode(samples) %in% c('integer', 'double'))) stop('prior_samples sampler \'samples\' control argument must be numeric or integer type') if(dim(samples)[2] != k) stop('prior_samples sampler \'samples\' control argument had ', dim(samples)[2], ' columns, but target nodes have ', k, ' scalar elements. These numbers must be equal.') - if(any(model$getNodeType(target) == 'stoch')) messageIfVerbose(' [Note] \'prior_samples\' sampler has been assigned to one or more stochastic nodes. The prior distribution for these nodes will be overridden by the prior samples.') + if(any(model$getNodeType(targetExpanded) == 'stoch')) messageIfVerbose(' [Note] \'prior_samples\' sampler has been assigned to one or more stochastic nodes. The prior distribution for these nodes will be overridden by the prior samples.') }, run = function() { - if(randomDraws) { - ind <<- ceiling(runif(1, 0, nSamples)) ## random draws + ind <- ceiling(runif(1, 0, nSamples)) ## random draw for proposal + values(model, targetExpanded) <<- samples[ind, 1:k] + logMHR <- model$calculateDiff(calcNodesNoSelf) + jump <- decide(logMHR) + if(jump) { + nimCopy(from = model, to = mvSaved, row = 1, nodes = targetExpanded, logProb = TRUE) + nimCopy(from = model, to = mvSaved, row = 1, nodes = copyNodesDeterm, logProb = FALSE) + nimCopy(from = model, to = mvSaved, row = 1, nodes = copyNodesStoch, logProbOnly = TRUE) } else { - ind <<- ind + 1 ## sequential draws (the default) - if(ind > nSamples) ind <<- 1 ## recycle sequential draws, if necessary + nimCopy(from = mvSaved, to = model, row = 1, nodes = targetExpanded, logProb = TRUE) + nimCopy(from = mvSaved, to = model, row = 1, nodes = copyNodesDeterm, logProb = FALSE) + nimCopy(from = mvSaved, to = model, row = 1, nodes = copyNodesStoch, logProbOnly = TRUE) } - values(model, targetExpanded) <<- samples[ind, 1:k] - model$calculate(calcNodes) - nimCopy(from = model, to = mvSaved, row = 1, nodes = target, logProb = FALSE) - nimCopy(from = model, to = mvSaved, row = 1, nodes = copyNodesDeterm, logProb = FALSE) - nimCopy(from = model, to = mvSaved, row = 1, nodes = copyNodesStoch, logProbOnly = TRUE) }, methods = list( - before_chain = function(MCMCniter = double(), MCMCnburnin = double(), MCMCchain = double()) { - ## issue a note if sequential draws from prior samples will have to be recycled - if((!randomDraws) & (MCMCniter > nSamples)) cat(' [Note] prior_samples sampler will recycle sequential draws from prior samples, since ', nSamples, ' samples were provided, and ', MCMCniter, ' MCMC iterations will be run.\n') - }, - reset = function() { - ind <<- 0 - } + reset = function() { } ) ) @@ -3990,20 +3984,17 @@ sampler_partial_mvn_pp <- nimbleFunction( #' #' @section prior_samples sampler: #' -#' The prior_samples sampler uses a provided set of numeric values (\code{samples}) to define the prior distribution of one or more model nodes. One every MCMC iteration, the prior_samples sampler takes value(s) from the numeric values provided, and stores these value(s) into the target model node(s). This allows one to define the prior distribution of model parameters empirically, using a set of numeric \code{samples}, presumably obtained previously using MCMC. The \code{target} node may be either a single scalar node (scalar case), or a collection of model nodes. -#' -#' The prior_samples sampler provides two options for selection of the value to use on each MCMC iteration. The default behaviour is to take sequential values from the \code{samples} vector (scalar case), or in the case of multiple dimensions, sequential rows of the \code{samples} matrix are used. The alternative behaviour, by setting the control argument \code{randomDraws = TRUE}, will instead use random draws from the \code{samples} vector (scalar case), or randomly selected rows of the \code{samples} matrix in the multidimensional case. +#' The prior_samples sampler uses a provided set of numeric values (\code{samples}) to define the prior distribution of one or more model nodes. Once every MCMC iteration, the prior_samples sampler randomly selects new value(s) from the numeric values provided, and probabilistically accepts these as the new values for the specific target node(s) using a standard Metropolis-Hastings accept/reject step. This allows one to define the prior distribution of model parameters empirically, using a set of numeric \code{samples}, presumably obtained previously using MCMC. #' -#' If the default of sequential selection of values is used, and the number of MCMC iterations exceeds the length of the \code{samples} vector (scalar case) or the number of rows of the \code{samples} matrix, then \code{samples} will be recycled as necessary for the number of MCMC iterations. A message to this effect is also printed at the beginning of the MCMC chain. +#' The \code{target} node may be either a single scalar node (scalar case), or a collection of model nodes (multidimensional case). In the scalar case, a propsed value is randomely selected from the \code{samples} vector. In the multidimensional case, a randomly selected rows of the \code{samples} matrix is used as the proposed values. #' #' Logically, prior_samples samplers might want to operate first, in advance of other samplers, on every MCMC iteration. By default, at the time of MCMC building, all prior_samples samplers are re-ordered to appear first in the list of MCMC samplers. This behavior can be subverted, however, by setting \code{nimbleOptions(MCMCorderPriorSamplesSamplersFirst = FALSE)}. #' -#' The prior_samples sampler can be assigned to non-stochastic model nodes (nodes which are not assigned a prior distribution in the model). In fact, it is recommended that nodes being assigned a prior_samples are not provided with a prior distribution in the model, and rather, that these nodes only appear on the right-hand-side of model declaration lines. In such case that a prior_samples sampler is assigned to a nodes with a prior distribution, the prior distribution will be overridden by the sample values provided to the sampler; however, the node will still be a stochastic node for other purposes, and will contribute to the model joint-density (using the sample values provided relative to the prior distribution), will have an MCMC sampler assigned to it by default, and also may introduce potential for confusion. In this case, a message is issued at the time of MCMC building. +#' The prior_samples sampler can be assigned to non-stochastic model nodes (nodes which are not assigned a prior distribution in the model). In fact, it is recommended that nodes being assigned a prior_samples are not provided with a prior distribution in the model, and rather, that these nodes only appear on the right-hand-side of model declaration lines. In such case that a prior_samples sampler is assigned to a nodes with a prior distribution, the prior distribution will be overridden by the \code{samples} provided to the sampler; however, the node will still be a stochastic node for other purposes, and will contribute to the model joint-density (using the current values relative to the prior distribution), will have an MCMC sampler assigned to it by default, and also may introduce potential for confusion. In this case, a message is issued at the time of MCMC building. #' #' The prior_samples sampler accepts the following control list elements: #' \itemize{ -#' \item \code{samples}. A numeric vector or matrix. When the \code{target} node is a single scalar-valued node, \code{samples} should be a numeric vector. When the \code{target} node specifies d > 2 model dimensions, \code{samples} should be a matrix containing d columns. The \code{samples} control argument is required. -#' \item \code{randomDraws}. A logical argument, specifying whether to use a random draw from \code{samples} on each iteration. If \code{samples} is a matrix, then a randomly-selected row of the \code{samples} matrix is used. When \code{FALSE}, sequential values (or sequential matrix rows) are used (default = \code{FALSE}). +#' \item \code{samples}. A numeric vector or matrix. When the \code{target} node is a single scalar-valued node, \code{samples} should be a numeric vector. When the \code{target} node specifies more \code{d} > 2 dimensions, \code{samples} should be a matrix containing \code{d} columns. The \code{samples} control argument is required. #' } #' #' @section posterior_predictive sampler: diff --git a/packages/nimble/tests/testthat/test-mcmc.R b/packages/nimble/tests/testthat/test-mcmc.R index ed6de077c..d6bca2e6c 100644 --- a/packages/nimble/tests/testthat/test-mcmc.R +++ b/packages/nimble/tests/testthat/test-mcmc.R @@ -2890,10 +2890,10 @@ test_that('prior_samples sampler operates correctly', { a ~ dnorm(0, 1) b ~ dnorm(0, 1) for(i in 1:4) { - y[i] <- c[i]^2 + y[i] ~ dnorm(c[i], sd = 1) } }) - Rmodel <- nimbleModel(code, inits = list(a=0, b=0, c=rep(0,4))) + Rmodel <- nimbleModel(code, inits = list(a=0, b=0, c=rep(0,4)), data = list(y = c(1,4,7,10))) ## conf <- configureMCMC(Rmodel) conf$addSampler('a', 'prior_samples') @@ -2923,25 +2923,11 @@ test_that('prior_samples sampler operates correctly', { nimbleOptions(MCMCorderPriorSamplesSamplersFirst = TRUE) nimbleOptions(MCMCorderPosteriorPredictiveSamplersLast = TRUE) ## - conf <- configureMCMC(Rmodel, nodes = NULL, monitors = c('a', 'c', 'y')) + conf <- configureMCMC(Rmodel, nodes = NULL, monitors = c('a', 'c')) conf$addSampler('a', 'prior_samples', samples = 1:10) conf$addSampler('c', 'prior_samples', samples = array(1:12, c(3,4))) expect_no_error(Rmcmc <- buildMCMC(conf)) - samples <- runMCMC(Rmcmc, 11) - expect_true(all(samples[,'a'] == c(1:10, 1))) - expect_true(all(samples[,'c[2]'] == c(rep(4:6, 3), 4, 5))) - expect_true(all(samples[,'y[2]'] == c(rep(4:6, 3), 4, 5)^2)) - ## - conf <- configureMCMC(Rmodel, nodes = NULL, monitors = c('a', 'c', 'y')) - conf$addSampler('a', 'prior_samples', samples = 1:10, randomDraws = TRUE) - conf$addSampler('c', 'prior_samples', samples = array(1:12, c(3,4)), randomDraws = TRUE) - expect_no_error(Rmcmc <- buildMCMC(conf)) - set.seed(0) - samples <- runMCMC(Rmcmc, 100) - expect_true(all(samples[,'a'] %in% 1:10)) - expect_true(all(range(samples[,'a']) == c(1,10))) - expect_true(all(samples[,'c[3]'] %in% 7:9)) - expect_true(all(samples[,'y[3]'] %in% c(49,64,81))) + expect_no_error(samples <- runMCMC(Rmcmc, 11)) }) test_that('assigning samplers to data and allowData argument', {