From c331871dd8a6834e98e39150c4c63f334ad10e6f Mon Sep 17 00:00:00 2001 From: Daniel Turek Date: Tue, 24 Sep 2024 15:34:18 -0400 Subject: [PATCH 01/16] added getMixedDataNodeNames method for models --- packages/nimble/R/BUGS_model.R | 34 ++++++++++++++++++++++++++-------- packages/nimble/inst/NEWS.md | 6 ++++++ 2 files changed, 32 insertions(+), 8 deletions(-) diff --git a/packages/nimble/R/BUGS_model.R b/packages/nimble/R/BUGS_model.R index 1fade697a..0ebdf1afd 100644 --- a/packages/nimble/R/BUGS_model.R +++ b/packages/nimble/R/BUGS_model.R @@ -499,6 +499,16 @@ Details: Multiple logical input arguments may be used simultaneously. For examp } return(ans) }, + getMixedDataNodeNames = function(returnType = 'names') { + multivariateStochBool <- sapply(modelDef$declInfo, + function(di) di$type == 'stoch' && grepl(':', deparse(di$targetExpr))) + multivariateStochIDs <- sapply(modelDef$declInfo[multivariateStochBool], `[[`, 'graphIDs') + isDataResult <- isDataFromGraphID(multivariateStochIDs, includeMixed = TRUE) ## values are in {0, 1, 2} + mixedDataIDs <- multivariateStochIDs[isDataResult == 2] + if(returnType == 'ids') return(mixedDataIDs) + if(returnType == 'names') return(modelDef$maps$graphID_2_nodeName[mixedDataIDs]) + stop('returnType argument to getMixedDataNodeNames was invalid') + }, safeUpdateValidValues = function(validValues, idsVec_only, idsVec_exclude) { if(!missing(idsVec_only) && !missing(idsVec_exclude)) stop() if(!missing(idsVec_only)) { @@ -783,15 +793,23 @@ Details: The variable or node names specified is expanded into a vector of model return(isDataFromGraphID(g_id)) }, - isDataFromGraphID = function(g_id){ - ## returns TRUE if any elements are flagged as data + isDataFromGraphID = function(g_id, includeMixed = FALSE) { + ## default behaviour: returns TRUE if any elements are flagged as data + ## when includeMixed=TRUE: 0 for FALSE, 1 for TRUE, or 2 for MIXED nodeNames <- modelDef$maps$graphID_2_nodeName[g_id] - ret <- unlist(lapply(as.list(nodeNames), - function(nn) - return(any(eval(parse(text=nn, keep.source = FALSE)[[1]], - envir=isDataEnv))))) - if(is.null(ret)) ret <- logical(0) - return(ret) + if(includeMixed) { + f <- function(nn) { + vals <- eval(parse(text=nn, keep.source = FALSE)[[1]], envir=isDataEnv) + if(!any(vals)) return(0) + if(all(vals)) return(1) + return(2) + } + } else { + f <- function(nn) return(any(eval(parse(text=nn, keep.source = FALSE)[[1]], envir=isDataEnv))) + } + ret <- unlist(lapply(as.list(nodeNames), f)) + if(is.null(ret)) ret <- logical(0) + return(ret) }, getDependenciesList = function(returnNames = TRUE, sort = TRUE) { diff --git a/packages/nimble/inst/NEWS.md b/packages/nimble/inst/NEWS.md index b8a9ac339..609292cff 100644 --- a/packages/nimble/inst/NEWS.md +++ b/packages/nimble/inst/NEWS.md @@ -1,3 +1,9 @@ +## USER LEVEL CHANGES + +- Added a new method `getMixedDataNodeNames` for model objects. This method + returns the node names (or, optionally, the graph ids) of any stochastic + multivariate model nodes which are partially, but not entirely, observed data. + # CHANGES IN VERSION 1.2.1 (July 2024) ## USER LEVEL CHANGES From a56dfd48df394d41f7271057707f260ea55dd903 Mon Sep 17 00:00:00 2001 From: Daniel Turek Date: Wed, 7 May 2025 13:01:53 -0400 Subject: [PATCH 02/16] removed comments from NEWS --- packages/nimble/inst/NEWS.md | 6 ------ 1 file changed, 6 deletions(-) diff --git a/packages/nimble/inst/NEWS.md b/packages/nimble/inst/NEWS.md index 609292cff..b8a9ac339 100644 --- a/packages/nimble/inst/NEWS.md +++ b/packages/nimble/inst/NEWS.md @@ -1,9 +1,3 @@ -## USER LEVEL CHANGES - -- Added a new method `getMixedDataNodeNames` for model objects. This method - returns the node names (or, optionally, the graph ids) of any stochastic - multivariate model nodes which are partially, but not entirely, observed data. - # CHANGES IN VERSION 1.2.1 (July 2024) ## USER LEVEL CHANGES From 1bf77cb7c8b321e7d8b04c838e66f414e5776208 Mon Sep 17 00:00:00 2001 From: Daniel Turek Date: Wed, 7 May 2025 13:06:49 -0400 Subject: [PATCH 03/16] partial_mvn (#1543) `partial_mvn` sampler, to have the code here for posterity. --- packages/nimble/DESCRIPTION | 2 +- packages/nimble/R/BUGS_model.R | 5 + packages/nimble/R/MCMC_configuration.R | 17 +- packages/nimble/R/MCMC_samplers.R | 59 +++- packages/nimble/R/initializeModel.R | 2 +- packages/nimble/man/initializeModel.Rd | 2 +- packages/nimble/man/samplers.Rd | 23 +- packages/nimble/tests/testthat/test-mcmc.R | 337 +++++++++++++++++++++ 8 files changed, 438 insertions(+), 9 deletions(-) diff --git a/packages/nimble/DESCRIPTION b/packages/nimble/DESCRIPTION index 95af3f19f..21d82cc54 100644 --- a/packages/nimble/DESCRIPTION +++ b/packages/nimble/DESCRIPTION @@ -150,4 +150,4 @@ Collate: nimble-package.r QuadratureGrids.R zzz.R -RoxygenNote: 7.2.3 +RoxygenNote: 7.3.2 diff --git a/packages/nimble/R/BUGS_model.R b/packages/nimble/R/BUGS_model.R index 0ebdf1afd..322485f4a 100644 --- a/packages/nimble/R/BUGS_model.R +++ b/packages/nimble/R/BUGS_model.R @@ -503,12 +503,17 @@ Details: Multiple logical input arguments may be used simultaneously. For examp multivariateStochBool <- sapply(modelDef$declInfo, function(di) di$type == 'stoch' && grepl(':', deparse(di$targetExpr))) multivariateStochIDs <- sapply(modelDef$declInfo[multivariateStochBool], `[[`, 'graphIDs') + if(length(multivariateStochIDs) == 0) multivariateStochIDs <- numeric() ## length=0 case, make a numeric vector isDataResult <- isDataFromGraphID(multivariateStochIDs, includeMixed = TRUE) ## values are in {0, 1, 2} mixedDataIDs <- multivariateStochIDs[isDataResult == 2] if(returnType == 'ids') return(mixedDataIDs) if(returnType == 'names') return(modelDef$maps$graphID_2_nodeName[mixedDataIDs]) stop('returnType argument to getMixedDataNodeNames was invalid') }, + isMixedData = function(nodeNames) { + ids <- expandNodeNames(nodeNames, returnType = 'ids') + return(isDataFromGraphID(ids, includeMixed = TRUE) == 2) + }, safeUpdateValidValues = function(validValues, idsVec_only, idsVec_exclude) { if(!missing(idsVec_only) && !missing(idsVec_exclude)) stop() if(!missing(idsVec_only)) { diff --git a/packages/nimble/R/MCMC_configuration.R b/packages/nimble/R/MCMC_configuration.R index c2958fc20..4b91bb0d5 100644 --- a/packages/nimble/R/MCMC_configuration.R +++ b/packages/nimble/R/MCMC_configuration.R @@ -42,6 +42,7 @@ samplerConf <- setRefClass( ) + ## NOTE: methods are documented as a "docstring" with each method - see 'removeSamplers' below. roxygen will automatically grab info from these docstrings and inject into the Rd in the Methods Section ## NOTE: including the name of the class in @aliases is important because by default we only get help("MCMCconf-class") and not help(MCMCconf) ## NOTE: the empty lines are important in the final formatting, so please don't remove any of them in your own help info @@ -190,6 +191,11 @@ print: A logical argument specifying whether to print the montiors and samplers. if(missing(nodes)) { nodes <- model$getNodeNames(stochOnly = TRUE, includeData = FALSE) # Check of all(model$isStoch(nodes)) is not needed in this case + # Check and adds any partially observed nodes + mixedDataNodeNames <- model$getMixedDataNodeNames() + if(length(mixedDataNodeNames)) { + nodes <- model$topologicallySortNodes(c(nodes, mixedDataNodeNames)) + } } else if(is.null(nodes) || length(nodes)==0) { nodes <- character(0) } else nodes <- filterOutDataNodes(nodes) ## configureMCMC *never* assigns samplers to data nodes @@ -369,6 +375,13 @@ For internal use. Adds default MCMC samplers to the specified nodes. ## for multivariate nodes, either add a conjugate sampler, RW_multinomial, or RW_block sampler if(nodeLength > 1) { + if(model$isMixedData(node)) { + if(nodeDist == 'dmnorm') { + thisControlList <- c(controlDefaultsArg, multivariateNodesAsScalars = multivariateNodesAsScalars) + addSampler(target = node, type = 'partial_mvn', control = thisControlList, allowData = TRUE) ; next + } + stop(paste0('The node ', node, ' is partially observed. NIMBLE only handles this case for multivariate normal distibutions.')) + } if(useConjugacy) { conjugacyResult <- conjugacyResultsAll[[node]] if(!is.null(conjugacyResult)) { @@ -729,7 +742,9 @@ For internal use only filterOutDataNodes = function(nodes) { nodes <- model$expandNodeNames(nodes) - return(nodes[!model$isData(nodes)]) + # We don't filter out partially observed data nodes + gIDs <- model$modelDef$nodeName2GraphIDs(nodes) + return(nodes[model$isDataFromGraphID(gIDs, includeMixed=TRUE) != 1]) }, removeSamplers = function(..., ind, print = FALSE) { diff --git a/packages/nimble/R/MCMC_samplers.R b/packages/nimble/R/MCMC_samplers.R index 9e55f2403..a65d48be4 100644 --- a/packages/nimble/R/MCMC_samplers.R +++ b/packages/nimble/R/MCMC_samplers.R @@ -3229,8 +3229,56 @@ sampler_polyagamma <- nimbleFunction( ) ) +#################################################################### +### partially observed multivariate normal sampler ################# +#################################################################### +#' @rdname samplers +#' @export +sampler_partial_mvn <- nimbleFunction( + name = 'sampler_partial_mvn', + contains = sampler_BASE, + setup = function(model, mvSaved, target, control) { + ## control list extraction + multivariateNodesAsScalars <- extractControlElement(control, 'multivariateNodesAsScalars', error = 'The control list must include the argument multivariateNodesAsScalars') + ## node list generation + targetAsScalar <- model$expandNodeNames(target, returnScalarComponents = TRUE) + partObsTru <- sapply(targetAsScalar, function(targetAsScalar) !eval(parse(text = targetAsScalar)[[1]], envir = model$isDataEnv)) + targetUnobservedComponents <- targetAsScalar[partObsTru] + ## nested function and function list definitions + samplerList <- nimbleFunctionList(sampler_BASE) + + if(multivariateNodesAsScalars) { + for(i in seq_along(targetUnobservedComponents)) { + samplerList[[i]] <- sampler_RW(model, mvSaved, targetUnobservedComponents[i], control) + } + } else { + if(length(targetUnobservedComponents) == 1) { + samplerList[[1]] <- sampler_RW(model, mvSaved, targetUnobservedComponents, control) + } else { + samplerList[[1]] <- sampler_RW_block(model, mvSaved, targetUnobservedComponents, control) + } + } + ## checks + if (model$getDistribution(target) != "dmnorm") stop(paste0('The node ', target, ' is parially observed. NIMBLE only handles this case for multivariate normal distibutions.')) + if (!model$isMixedData(target)) stop(paste0('The target node ', target, ' is not partially observed.')) + }, + + run = function() { + for (i in seq_along(samplerList)) { + samplerList[[i]]$run() + } + }, + methods = list( + reset = function() { + for (i in seq_along(samplerList)) { + samplerList[[i]]$reset() + } + } + ) +) + #' MCMC Sampling Algorithms #' #' Details of the MCMC sampling algorithms provided with the NIMBLE MCMC engine; HMC samplers are in the \code{nimbleHMC} package and particle filter samplers are in the \code{nimbleSMC} package. @@ -3618,6 +3666,15 @@ sampler_polyagamma <- nimbleFunction( #' The posterior_predictive sampler functions by simulating new values for all downstream (dependent) nodes using their conditional distributions, as well as updating the associated model probabilities. A posterior_predictive sampler will automatically be assigned to all trailing non-data stochastic nodes in a model, or when possible, to any node at a point in the model after which all downstream (dependent) stochastic nodes are non-data. #' #' The posterior_predictive sampler accepts no control list arguments. +#' +#' @section partial_mvn sampler: +#' +#' The partial_mvn sampler is designed to sample multivariate normal distributions that are partially observed. That is, some dimensions of the target node are observed data values, some dimensions are not data. Sampling is accomplished using either univariate or multivariate random walk Metropolis Hastings of the unobserved dimensions, as determined by the \code{multivariateNodesAsScalars} argument. +#' +#' The \code{partial_mvn} sampler accepts the following control list elements: +#' \itemize{ +#' \item multivariateNodesAsScalars. A logical argument, specifying whether the sampler should sample the unobserved parts of a partially observed node jointly or independently (default = FALSE). +#' } #' #' @section RJ_fixed_prior sampler: #' @@ -3633,7 +3690,7 @@ sampler_polyagamma <- nimbleFunction( #' #' @name samplers #' -#' @aliases sampler binary categorical prior_samples posterior_predictive RW RW_block RW_multinomial RW_dirichlet RW_wishart RW_llFunction slice AF_slice crossLevel RW_llFunction_block sampler_prior_samples sampler_posterior_predictive sampler_binary sampler_categorical sampler_RW sampler_RW_block sampler_RW_multinomial sampler_RW_dirichlet sampler_RW_wishart sampler_RW_llFunction sampler_slice sampler_AF_slice sampler_crossLevel sampler_RW_llFunction_block CRP CRP_concentration DPmeasure RJ_fixed_prior RJ_indicator RJ_toggled RW_PF RW_PF_block RW_lkj_corr_cholesky sampler_RW_lkj_corr_cholesky RW_block_lkj_corr_cholesky sampler_RW_block_lkj_corr_cholesky +#' @aliases sampler binary categorical prior_samples posterior_predictive RW RW_block RW_multinomial RW_dirichlet RW_wishart RW_llFunction slice AF_slice crossLevel RW_llFunction_block sampler_prior_samples sampler_posterior_predictive sampler_binary sampler_categorical sampler_RW sampler_RW_block sampler_RW_multinomial sampler_RW_dirichlet sampler_RW_wishart sampler_RW_llFunction sampler_slice sampler_AF_slice sampler_crossLevel sampler_RW_llFunction_block CRP CRP_concentration DPmeasure RJ_fixed_prior RJ_indicator RJ_toggled RW_PF RW_PF_block RW_lkj_corr_cholesky sampler_RW_lkj_corr_cholesky RW_block_lkj_corr_cholesky sampler_RW_block_lkj_corr_cholesky partial_mvn sampler_partial_mvn #' #' @examples #' ## y[1] ~ dbern() or dbinom(): diff --git a/packages/nimble/R/initializeModel.R b/packages/nimble/R/initializeModel.R index 316b076b9..70f92ca91 100644 --- a/packages/nimble/R/initializeModel.R +++ b/packages/nimble/R/initializeModel.R @@ -7,7 +7,7 @@ #' @author Daniel Turek #' @details This nimbleFunction may be used at the beginning of nimble algorithms to perform model initialization. #' The intended usage is to specialize an instance of this nimbleFunction in the setup function of an algorithm, -#' then execute that specialied function at the beginning of the algorithm run function. +#' then execute that specialized function at the beginning of the algorithm run function. #' The specialized function takes no arguments. #' #' Executing this function ensures that all right-hand-side only nodes have been assigned real values, diff --git a/packages/nimble/man/initializeModel.Rd b/packages/nimble/man/initializeModel.Rd index 04eea8d22..18a088d80 100644 --- a/packages/nimble/man/initializeModel.Rd +++ b/packages/nimble/man/initializeModel.Rd @@ -17,7 +17,7 @@ Performs initialization of nimble model node values and log probabilities \details{ This nimbleFunction may be used at the beginning of nimble algorithms to perform model initialization. The intended usage is to specialize an instance of this nimbleFunction in the setup function of an algorithm, -then execute that specialied function at the beginning of the algorithm run function. +then execute that specialized function at the beginning of the algorithm run function. The specialized function takes no arguments. Executing this function ensures that all right-hand-side only nodes have been assigned real values, diff --git a/packages/nimble/man/samplers.Rd b/packages/nimble/man/samplers.Rd index e34347840..db5cfeccb 100644 --- a/packages/nimble/man/samplers.Rd +++ b/packages/nimble/man/samplers.Rd @@ -23,6 +23,7 @@ \alias{sampler_CAR_normal} \alias{sampler_CAR_proper} \alias{sampler_polyagamma} +\alias{sampler_partial_mvn} \alias{samplers} \alias{sampler} \alias{binary} @@ -49,6 +50,7 @@ \alias{RW_PF_block} \alias{RW_lkj_corr_cholesky} \alias{RW_block_lkj_corr_cholesky} +\alias{partial_mvn} \alias{sampler_RJ_fixed_prior} \alias{sampler_RJ_indicator} \alias{sampler_RJ_toggled} @@ -101,6 +103,8 @@ sampler_CAR_proper(model, mvSaved, target, control) sampler_polyagamma(model, mvSaved, target, control) +sampler_partial_mvn(model, mvSaved, target, control) + sampler_RJ_fixed_prior(model, mvSaved, target, control) sampler_RJ_indicator(model, mvSaved, target, control) @@ -377,7 +381,7 @@ Note that this sampler is likely run much more slowly than the blocked sampler f \section{polyagamma sampler}{ -The polyagamma sampler uses Pólya-gamma data augmentation to do conjugate sampling for the parameters in the linear predictor of a logistic regression model (Polson et al., 2013), analogous to the Albert-Chib data augmentation scheme for probit regression. This sampler is not assigned as a default sampler by \code{configureMCMC} and so can only be used if manually added to an MCMC configuration. +The polyagamma sampler uses \enc{Pólya}{Polya}-gamma data augmentation to do conjugate sampling for the parameters in the linear predictor of a logistic regression model (Polson et al., 2013), analogous to the Albert-Chib data augmentation scheme for probit regression. This sampler is not assigned as a default sampler by \code{configureMCMC} and so can only be used if manually added to an MCMC configuration. As an example, consider model code containing: \preformatted{for(i in 1:n) { @@ -398,9 +402,9 @@ MCMCconf$addSampler(target = logistic_nodes, type = "polyagamma", control = list(fixedDesignColumns=TRUE)) } -As shown here, the stochastic dependencies (\code{y[i]} here) of the target nodes must follow \code{dbin} or \code{dbern} distributions. The logit transformation of their probability parameter must be a linear function (technically an affine function) of the target nodes, which themselves must have \code{dnorm} or \code{dmnorm} priors. Zero inflation to account for structural zeroes is also supported, allowed as discussed below. The stochastic dependencies will often but not always be the observations in the logistic regression and will be referred to as 'responses' henceforth. Internally, the sampler draws latent values from the Pólya-gamma distribution, one per response. These latent values are then used to draw from the multivariate normal conditional distribution of the target nodes. +As shown here, the stochastic dependencies (\code{y[i]} here) of the target nodes must follow \code{dbin} or \code{dbern} distributions. The logit transformation of their probability parameter must be a linear function (technically an affine function) of the target nodes, which themselves must have \code{dnorm} or \code{dmnorm} priors. Zero inflation to account for structural zeroes is also supported, allowed as discussed below. The stochastic dependencies will often but not always be the observations in the logistic regression and will be referred to as 'responses' henceforth. Internally, the sampler draws latent values from the \enc{Pólya}{Polya}-gamma distribution, one per response. These latent values are then used to draw from the multivariate normal conditional distribution of the target nodes. -Importantly, note that because the Pólya-gamma draws are not retained when an iteration of the sampler finishes, one generally wants to apply the sampler to all parameter nodes involved in the linear predictor of the logistic regression, to avoid duplicative Pólya-gamma draws of the latent values. If there are stochastic indices (e.g., if \code{group[i]} above is stochastic), the Pólya-gamma sampler can still be used, but the stochastic nodes cannot be sampled by it and must have separate sampler(s). It is also possible in some models that regression parameters can be split into (conditionally independent) groups that can be sampled independently, e.g., if one has distinct logistic regression specifications for different sets of responses in the model. +Importantly, note that because the \enc{Pólya}{Polya}-gamma draws are not retained when an iteration of the sampler finishes, one generally wants to apply the sampler to all parameter nodes involved in the linear predictor of the logistic regression, to avoid duplicative \enc{Pólya}{Polya}-gamma draws of the latent values. If there are stochastic indices (e.g., if \code{group[i]} above is stochastic), the \enc{Pólya}{Polya}-gamma sampler can still be used, but the stochastic nodes cannot be sampled by it and must have separate sampler(s). It is also possible in some models that regression parameters can be split into (conditionally independent) groups that can be sampled independently, e.g., if one has distinct logistic regression specifications for different sets of responses in the model. Sampling involves use of the design matrix. The design matrix includes one column corresponding to each regression covariate as well as one columnar block corresponding to each random effect. In the example above, the columns would include a vector of ones (to multiply \code{beta0}), the vectors \code{x1} and \code{x2}, and a vector of indicators for each \code{u[j]} (with a 1 in row \code{i} if \code{group[i]} is \code{j}), resulting in \code{3+num_groups} columns. Note that the polyagamma sampler can determine the design matrix from the model, even when written as above such that the design matrix is not explicitly in the model code. It is also possible to write model code for the linear prediction using matrix multiplication and an explicit design matrix, but that is not necessary. @@ -552,6 +556,17 @@ The posterior_predictive sampler functions by simulating new values for all down The posterior_predictive sampler accepts no control list arguments. } +\section{partial_mvn sampler}{ + + +The partial_mvn sampler is designed to sample multivariate normal distributions that are partially observed. That is, some dimensions of the target node are observed data values, some dimensions are not data. Sampling is accomplished using either univariate or multivariate random walk Metropolis Hastings of the unobserved dimensions, as determined by the \code{multivariateNodesAsScalars} argument. + +The \code{partial_mvn} sampler accepts the following control list elements: +\itemize{ +\item multivariateNodesAsScalars. A logical argument, specifying whether the sampler should sample the unobserved parts of a partially observed node jointly or independently (default = FALSE). +} +} + \section{RJ_fixed_prior sampler}{ @@ -621,7 +636,7 @@ Neal, R. M. (2011). MCMC Using Hamiltonian Dynamics. \emph{Handbook of Markov Ch Pitt, M. K. and Shephard, N. (1999). Filtering via simulation: Auxiliary particle filters. \emph{Journal of the American Statistical Association} 94(446), 590-599. -Polson, N.G., Scott, J.G., and J. Windle. (2013). Bayesian inference for logistic models using Pólya-gamma latent variables. Journal of the American Statistical Association, 108(504), 1339–1349. https://doi.org/10.1080/01621459.2013.829001 +Polson, N.G., Scott, J.G., and J. Windle. (2013). Bayesian inference for logistic models using \enc{Pólya}{Polya}-gamma latent variables. Journal of the American Statistical Association, 108(504), 1339–1349. https://doi.org/10.1080/01621459.2013.829001 Roberts, G. O. and S. K. Sahu (1997). Updating Schemes, Correlation Structure, Blocking and Parameterization for the Gibbs Sampler. \emph{Journal of the Royal Statistical Society: Series B (Statistical Methodology)}, 59(2), 291-317. diff --git a/packages/nimble/tests/testthat/test-mcmc.R b/packages/nimble/tests/testthat/test-mcmc.R index 6827bf59d..13fd0cf48 100644 --- a/packages/nimble/tests/testthat/test-mcmc.R +++ b/packages/nimble/tests/testthat/test-mcmc.R @@ -3022,6 +3022,343 @@ test_that('assigning samplers to data and allowData argument', { expect_true(samps[[4]]$name == 'posterior_predictive') }) +test_that('partial_mvn sampler was given to dmnorm distribution when dependent on a node without being depended on by another node', { + code <- nimbleCode({ for (i in 1:N){ + theta[i] ~ dgamma(alpha,beta) + lambda[i] <- theta[i]*t[i] + x[i] ~ dpois(lambda[i]) + } + alpha ~ dexp(1.0) + beta ~ dgamma(0.1,1.0) + for(i in 1:5) { + mu[i] <- alpha+i + } + y[1:5]~dmnorm(mu[1:5], Sigma[1:5,1:5]) + }) + Consts <- list(N = 10, t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5)) + Data <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22), y=c(rep(NA,3),4,5)) + + model <- nimbleModel(code = code, name = "model", constants = Consts, data = Data) + + conf <- configureMCMC(model, nodes = 'y[1:5]') + + expect_true(any(sapply(conf$getSamplers(), function(sc) sc$name)=='partial_mvn', info= "partial_mvn sampler not assigned to dmnorm distribution")) +}) + + +test_that('partial_mvn sampler was given to dmnorm distribution when not dependent on nodes in the model', { + code <- nimbleCode({ for (i in 1:N){ + theta[i] ~ dnorm(alpha,beta) + lambda[i] <- theta[i]*t[i] + x[i] ~ dexp(lambda[i]) + } + alpha ~ dpois(1.0) + beta ~ dunif(0.1,1.0) + y[1:10]~dmnorm(mu[1:10], Sigma[1:10,1:10]) + }) + consts <- list(N = 10, t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5)) + data <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22), y=c(rep(NA,3),4:10)) + + model <- nimbleModel(code = code, name = "model", constants = consts, data = data) + + conf <- configureMCMC(model, nodes = 'y[1:10]') + + expect_true(any(sapply(conf$getSamplers(), function(sc) sc$name)=='partial_mvn', info= "partial_mvn sampler not assigned to dmnorm distribution")) +}) + +test_that('partial_mvn sampler was given to dmnorm distribution when intermediate in a model', { + code <- nimbleCode({ for (i in 1:N){ + theta[i] ~ dgamma(alpha,beta) + lambda[i] <- theta[i]*t[i] + x[i] ~ dpois(lambda[i]) + } + alpha ~ dexp(1.0) + beta ~ dgamma(0.1,1.0) + for(i in 1:5) { + mu[i] <- alpha+i + } + y[1:5]~dmnorm(mu[1:5], Sigma[1:5,1:5]) + h[1:5]~dmvt(y[1:5], Sigma[1:5,1:5], df = 10) + }) + Consts <- list(N = 10, t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5)) + Data <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22), y=c(rep(NA,3),4,5)) + + model <- nimbleModel(code = code, name = "model", constants = Consts, data = Data) + + conf <- configureMCMC(model, nodes = 'y[1:5]') + + expect_true(any(sapply(conf$getSamplers(), function(sc) sc$name)=='partial_mvn', info= "partial_mvn sampler not assigned to dmnorm distribution")) +}) + +test_that('partial_mvn sampler was given to each dmnorm distribution', { + code <- nimbleCode({ for (i in 1:N){ + theta[i] ~ dgamma(alpha,beta) + lambda[i] <- theta[i]*t[i] + x[i] ~ dpois(lambda[i]) + } + alpha ~ dexp(1.0) + beta ~ dgamma(0.1,1.0) + for(i in 1:5) { + mu[i] <- alpha+i + me[i] <- beta-i + } + y[1:5]~dmnorm(mu[1:5], Sigma[1:5,1:5]) + z[1:5]~dmnorm(me[1:5], Sigmo[1:5,1:5]) + }) + Consts <- list(N = 10, t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5)) + Data <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22), y=c(rep(NA,3),4,5), z=c(1:3, rep(NA,2))) + + model <- nimbleModel(code = code, name = "model", constants = Consts, data = Data) + + conf <- configureMCMC(model, nodes = model$getNodeNames()[26:27]) + + expect_true(sum(sapply(conf$getSamplers(), function(sc) sc$name)=='partial_mvn')==2, info= "partial_mvn sampler not assigned to each dmnorm distribution") +}) + +test_that('partial_mvn sampler was given to each partially obs dmnorm distribution when both are at intermediate level in a model', { + code <- nimbleCode({ for (i in 1:N){ + theta[i] ~ dgamma(alpha,beta) + lambda[i] <- theta[i]*t[i] + x[i] ~ dpois(lambda[i]) + } + alpha ~ dexp(1.0) + beta ~ dgamma(0.1,1.0) + for(i in 1:5) { + mu[i] <- alpha+i + me[i] <- beta-i + } + y[1:5]~dmnorm(mu[1:5], Sigma[1:5,1:5]) + z[1:5]~dmnorm(me[1:5], Sigmo[1:5,1:5]) + + a[1:5]~dmvt(y[1:5], Sigma[1:5,1:5], df = 10) + b[1:5]~dmvt(z[1:5], Sigmo[1:5,1:5], df = 10) + }) + Consts <- list(N = 10, t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5)) + Data <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22), y=c(rep(NA,3),4,5), z=c(1:3, rep(NA,2))) + + model <- nimbleModel(code = code, name = "model", constants = Consts, data = Data) + + conf <- configureMCMC(model, nodes = model$getNodeNames()[26:27]) + + expect_true(sum(sapply(conf$getSamplers(), function(sc) sc$name)=='partial_mvn')==2, info= "partial_mvn sampler not assigned to each dmnorm distribution") +}) + +test_that('an error is given when trying to assign sampler to partially observed dmvt dist', { + code <- nimbleCode({ for (i in 1:N){ + theta[i] ~ dgamma(alpha,beta) + lambda[i] <- theta[i]*t[i] + x[i] ~ dpois(lambda[i]) + } + alpha ~ dexp(1.0) + beta ~ dgamma(0.1,1.0) + y[1:5]~dmvt(mu[1:5], Sigma[1:5,1:5], df=10) + }) + Consts <- list(N = 10, t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5)) + Data <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22), y=c(rep(NA,3),4,5)) + + model <- nimbleModel(code = code, name = "model", constants = Consts, data = Data) + + expect_error(configureMCMC(model, nodes = 'y[1:5]', print = FALSE), info = "partially observed dmvt node given a sampler") + expect_error(configureMCMC(model, print = FALSE), info = "model with partially observed dmvt given a sampler for every node" ) +}) + +test_that('an error is given when trying to assign a sampler to partially observed dmvt dist when it is intermediate in a model', { + code <- nimbleCode({ for (i in 1:N){ + theta[i] ~ dgamma(alpha,beta) + lambda[i] <- theta[i]*t[i] + x[i] ~ dpois(lambda[i]) + } + alpha ~ dexp(1.0) + beta ~ dgamma(0.1,1.0) + for(i in 1:5) { + mu[i] <- alpha+i + } + y[1:5]~dmvt(mu[1:5], Sigma[1:5,1:5], df = 10) + h[1:5]~dmvt(y[1:5], Sigma[1:5,1:5], df = 10) + }) + Consts <- list(N = 10, t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5)) + Data <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22), y=c(rep(NA,3),4,5)) + + model <- nimbleModel(code = code, name = "model", constants = Consts, data = Data) + + expect_error(configureMCMC(model, print = FALSE), info = "model with dmvt given a sampler for every node" ) +}) + +test_that('trying to give a sampler to an observed node doesnt yield an error from configureMCMC', { + code <- nimbleCode({ for (i in 1:N){ + theta[i] ~ dgamma(alpha,beta) + lambda[i] <- theta[i]*t[i] + x[i] ~ dpois(lambda[i]) + } + alpha ~ dexp(1.0) + beta ~ dgamma(0.1,1.0) + y[1:5]~dmvt(mu[1:5], Sigma[1:5,1:5], df=10) + }) + Consts <- list(N = 10, t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5)) + Data <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22), y=c(rep(NA,3),4,5)) + + model <- nimbleModel(code = code, name = "model", constants = Consts, data = Data) + + expect_error(configureMCMC(model, nodes = 'alpha', print = FALSE), NA, info = "observed node in model given a sampler") +}) + +test_that('partial_mvn sampler was not given to non dmnorm dist', { + code <- nimbleCode({ for (i in 1:N){ + theta[i] ~ dgamma(alpha,beta) + lambda[i] <- theta[i]*t[i] + x[i] ~ dpois(lambda[i]) + } + alpha ~ dexp(1.0) + beta ~ dgamma(0.1,1.0) + }) + Consts <- list(N = 10, t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5)) + Data <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22)) + + model <- nimbleModel(code = code, name = "model", constants = Consts, data = Data) + + conf <- configureMCMC(model) + + expect_true(!any(sapply(conf$getSamplers(), function(sc) sc$name)=='partial_mvn'), info = "partial_mvn sampler assigned as sampler to non dmnorm distribution") +}) + +test_that('partial_mvn sampler was not given to non dmnorm distribution from model that has dmnorm distribution', { + code <- nimbleCode({ for (i in 1:N){ + theta[i] ~ dnorm(alpha,beta) + lambda[i] <- theta[i]*t[i] + x[i] ~ dexp(lambda[i]) + } + alpha ~ dpois(1.0) + beta ~ dunif(0.1,1.0) + y[1:10]~dmnorm(mu[1:10], Sigma[1:10,1:10]) + }) + consts <- list(N = 10, t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5)) + data <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22), y=c(rep(NA,3),4:10)) + + model <- nimbleModel(code = code, name = "model", constants = consts, data = data) + + conf <- configureMCMC(model) + + expect_true(!any(sapply(conf$getSamplers(), function(sc) sc$name)[head(seq_along(sapply(conf$getSamplers(), function(sc) sc$name)), -1)]=='partial_mvn'), info="partial_mvn sampler assigned as sampler to non dmnorm node in model with dmnorm") +}) + +test_that('Values change at different rows for unobserved parts of partially observed node when multivariateNodesAsScalars = TRUE', { + Code <- nimbleCode({ for (i in 1:N){ + theta[i] ~ dgamma(alpha,beta) + lambda[i] <- theta[i]*t[i] + x[i] ~ dpois(lambda[i]) + } + alpha ~ dexp(1.0) + beta ~ dgamma(0.1,1.0) + y[1:5]~dmnorm(mu[1:5], Sigma[1:5,1:5]) + }) + Consts <- list(N = 10, t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5)) + Data <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22), y=c(rep(NA,3),4,5)) + Inits <- list(alpha = 1, beta = 1, theta = rep(0.1, Consts$N), mu=1:5, Sigma=diag(5), y=c(rep(1,3),rep(NA,2))) + model <- nimbleModel(code = Code, name = "model", constants = Consts, data = Data, inits = Inits) + + Cmodel <- compileNimble(model) + + conf <- configureMCMC(model, nodes = 'y[1:5]', multivariateNodesAsScalars = TRUE) + conf$addMonitors("y") + + modelMCMC <- buildMCMC(conf) + CmodelMCMC <- compileNimble(modelMCMC, project = model) + modrun <- runMCMC(CmodelMCMC, niter=100, setSeed = 0) + + + onezro=apply(modrun[,3:5], 2, diff)!=0 + + expect_true(any(apply(onezro, 1, function(x) any(x) && !all(x))), info= "Values of MCMC sampled nodes change at same time as other nodes") +}) + +test_that('Values change at same rows for unobserved parts of partially observed node when multivariateNodesAsScalars = FALSE', { + Code <- nimbleCode({ for (i in 1:N){ + theta[i] ~ dgamma(alpha,beta) + lambda[i] <- theta[i]*t[i] + x[i] ~ dpois(lambda[i]) + } + alpha ~ dexp(1.0) + beta ~ dgamma(0.1,1.0) + y[1:5]~dmnorm(mu[1:5], Sigma[1:5,1:5]) + }) + Consts <- list(N = 10, t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5)) + Data <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22), y=c(rep(NA,3),4,5)) + Inits <- list(alpha = 1, beta = 1, theta = rep(0.1, Consts$N), mu=1:5, Sigma=diag(5), y=c(rep(1,3),rep(NA,2))) + model <- nimbleModel(code = Code, name = "model", constants = Consts, data = Data, inits = Inits) + + Cmodel <- compileNimble(model) + + conf <- configureMCMC(model, nodes = 'y[1:5]', multivariateNodesAsScalars = FALSE) + conf$addMonitors("y") + + modelMCMC <- buildMCMC(conf) + CmodelMCMC <- compileNimble(modelMCMC, project = model) + modrun <- runMCMC(CmodelMCMC, niter=100, setSeed = 0) + + onezro <- apply(modrun[,3:5], 2, diff)!=0 + + expect_true(all(apply(onezro, 1, function(x) all(x) || all(!x))), info= "Values of MCMC sampled nodes change at same time as other nodes") +}) + + +test_that('Values do not change for data when multivariateNodesAsScalars = TRUE', { + Code <- nimbleCode({ for (i in 1:N){ + theta[i] ~ dgamma(alpha,beta) + lambda[i] <- theta[i]*t[i] + x[i] ~ dpois(lambda[i]) + } + alpha ~ dexp(1.0) + beta ~ dgamma(0.1,1.0) + y[1:5]~dmnorm(mu[1:5], Sigma[1:5,1:5]) + }) + Consts <- list(N = 10, t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5)) + Data <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22), y=c(rep(NA,3),4,5)) + Inits <- list(alpha = 1, beta = 1, theta = rep(0.1, Consts$N), mu=1:5, Sigma=diag(5), y=c(rep(1,3),rep(NA,2))) + model <- nimbleModel(code = Code, name = "model", constants = Consts, data = Data, inits = Inits) + + Cmodel <- compileNimble(model) + + conf <- configureMCMC(model, nodes = 'y[1:5]', multivariateNodesAsScalars = TRUE) + conf$addMonitors("y") + + modelMCMC <- buildMCMC(conf) + CmodelMCMC <- compileNimble(modelMCMC, project = model) + modrun <- runMCMC(CmodelMCMC, niter=100, setSeed = 0) + + onezro <- apply(modrun[,c(1,2,6,7)], 2, diff)!=0 + + expect_true(all(apply(onezro, 1, function(x) !any(x))), info= "Values of data do not change when MCMC sampling") +}) + + +test_that('Values do not change for data when multivariateNodesAsScalars = FALSE', { + Code <- nimbleCode({ for (i in 1:N){ + theta[i] ~ dgamma(alpha,beta) + lambda[i] <- theta[i]*t[i] + x[i] ~ dpois(lambda[i]) + } + alpha ~ dexp(1.0) + beta ~ dgamma(0.1,1.0) + y[1:5]~dmnorm(mu[1:5], Sigma[1:5,1:5]) + }) + Consts <- list(N = 10, t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5)) + Data <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22), y=c(rep(NA,3),4,5)) + Inits <- list(alpha = 1, beta = 1, theta = rep(0.1, Consts$N), mu=1:5, Sigma=diag(5), y=c(rep(1,3),rep(NA,2))) + model <- nimbleModel(code = Code, name = "model", constants = Consts, data = Data, inits = Inits) + + Cmodel <- compileNimble(model) + + conf <- configureMCMC(model, nodes = 'y[1:5]', multivariateNodesAsScalars = FALSE) + conf$addMonitors("y") + + modelMCMC <- buildMCMC(conf) + CmodelMCMC <- compileNimble(modelMCMC, project = model) + modrun <- runMCMC(CmodelMCMC, niter=100, setSeed = 0) + + onezro <- apply(modrun[,c(1,2,6,7)], 2, diff)!=0 + + expect_true(all(apply(onezro, 1, function(x) !any(x))), info= "Values of data do not change when MCMC sampling") +}) sink(NULL) From 536c62f6158b7b0405d740bf0b24dfcfbac28ca6 Mon Sep 17 00:00:00 2001 From: Daniel Turek Date: Mon, 26 Jan 2026 13:15:08 -0500 Subject: [PATCH 04/16] fixed unlist for graph node ids --- packages/nimble/R/BUGS_model.R | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/packages/nimble/R/BUGS_model.R b/packages/nimble/R/BUGS_model.R index fc0cc0536..2cf6e609b 100644 --- a/packages/nimble/R/BUGS_model.R +++ b/packages/nimble/R/BUGS_model.R @@ -502,9 +502,10 @@ Details: Multiple logical input arguments may be used simultaneously. For examp getMixedDataNodeNames = function(returnType = 'names') { multivariateStochBool <- sapply(modelDef$declInfo, function(di) di$type == 'stoch' && grepl(':', deparse(di$targetExpr))) - multivariateStochIDs <- sapply(modelDef$declInfo[multivariateStochBool], `[[`, 'graphIDs') + multivariateStochIDsList <- lapply(modelDef$declInfo[multivariateStochBool], `[[`, 'graphIDs') + multivariateStochIDs <- unlist(multivariateStochIDsList) if(length(multivariateStochIDs) == 0) multivariateStochIDs <- numeric() ## length=0 case, make a numeric vector - isDataResult <- isDataFromGraphID(multivariateStochIDs, includeMixed = TRUE) ## values are in {0, 1, 2} + isDataResult <- isDataFromGraphID(multivariateStochIDs, includeMixed = TRUE) ## values are in {0, 1, 2}, where 0=FALSE, 1=TRUE, 2=MIXED mixedDataIDs <- multivariateStochIDs[isDataResult == 2] if(returnType == 'ids') return(mixedDataIDs) if(returnType == 'names') return(modelDef$maps$graphID_2_nodeName[mixedDataIDs]) From 9a40789d9a9cc63877a268c2a86a0c86a0638f2b Mon Sep 17 00:00:00 2001 From: Daniel Turek Date: Tue, 10 Feb 2026 15:08:05 -0500 Subject: [PATCH 05/16] updates --- packages/nimble/R/MCMC_samplers.R | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/packages/nimble/R/MCMC_samplers.R b/packages/nimble/R/MCMC_samplers.R index 8885efbe6..07290c9a5 100644 --- a/packages/nimble/R/MCMC_samplers.R +++ b/packages/nimble/R/MCMC_samplers.R @@ -3467,11 +3467,15 @@ sampler_partial_mvn <- nimbleFunction( contains = sampler_BASE, setup = function(model, mvSaved, target, control) { ## control list extraction - multivariateNodesAsScalars <- extractControlElement(control, 'multivariateNodesAsScalars', error = 'The control list must include the argument multivariateNodesAsScalars') + browser() + multivariateNodesAsScalars <- extractControlElement(control, 'multivariateNodesAsScalars', default = getNimbleOption('MCMCmultivariateNodesAsScalars')) ## node list generation targetAsScalar <- model$expandNodeNames(target, returnScalarComponents = TRUE) - partObsTru <- sapply(targetAsScalar, function(targetAsScalar) !eval(parse(text = targetAsScalar)[[1]], envir = model$isDataEnv)) - targetUnobservedComponents <- targetAsScalar[partObsTru] + isdataBool <- sapply(targetAsScalar, function(targetAsScalar) eval(parse(text = targetAsScalar)[[1]], envir = model$isDataEnv)) + targetNonDataComponents <- targetAsScalar[!isdataBool] + predictiveBool <- sapply(targetNonDataComponents, function(n) length(model$getDependencies(n, self=FALSE, downstream=TRUE, dataOnly=TRUE)) == 0) + targetNonDataPP <- targetNonDataComponents[ predictiveBool] + targetNonDataNP <- targetNonDataComponents[!predictiveBool] ## nested function and function list definitions samplerList <- nimbleFunctionList(sampler_BASE) if(multivariateNodesAsScalars) { @@ -3486,8 +3490,8 @@ sampler_partial_mvn <- nimbleFunction( } } ## checks - if (model$getDistribution(target) != "dmnorm") stop(paste0('The node ', target, ' is parially observed. NIMBLE only handles this case for multivariate normal distibutions.')) - if (!model$isMixedData(target)) stop(paste0('The target node ', target, ' is not partially observed.')) + if (model$getDistribution(target) != "dmnorm") stop(paste0('The node ', target, ' is parially observed. NIMBLE only handles this case for multivariate normal distibutions.')) + if (!model$isMixedData(target)) stop(paste0('The target node ', target, ' is not partially observed.')) }, run = function() { for (i in seq_along(samplerList)) { From 1b34ab7be78538ae0fa58cc0eb5726ea2f0d7829 Mon Sep 17 00:00:00 2001 From: Daniel Turek Date: Tue, 10 Feb 2026 22:02:25 -0500 Subject: [PATCH 06/16] updates --- packages/nimble/R/MCMC_samplers.R | 86 +++++++++++++++++++++++++------ 1 file changed, 71 insertions(+), 15 deletions(-) diff --git a/packages/nimble/R/MCMC_samplers.R b/packages/nimble/R/MCMC_samplers.R index 07290c9a5..707f65c67 100644 --- a/packages/nimble/R/MCMC_samplers.R +++ b/packages/nimble/R/MCMC_samplers.R @@ -3467,46 +3467,102 @@ sampler_partial_mvn <- nimbleFunction( contains = sampler_BASE, setup = function(model, mvSaved, target, control) { ## control list extraction - browser() multivariateNodesAsScalars <- extractControlElement(control, 'multivariateNodesAsScalars', default = getNimbleOption('MCMCmultivariateNodesAsScalars')) ## node list generation targetAsScalar <- model$expandNodeNames(target, returnScalarComponents = TRUE) isdataBool <- sapply(targetAsScalar, function(targetAsScalar) eval(parse(text = targetAsScalar)[[1]], envir = model$isDataEnv)) targetNonDataComponents <- targetAsScalar[!isdataBool] predictiveBool <- sapply(targetNonDataComponents, function(n) length(model$getDependencies(n, self=FALSE, downstream=TRUE, dataOnly=TRUE)) == 0) - targetNonDataPP <- targetNonDataComponents[ predictiveBool] - targetNonDataNP <- targetNonDataComponents[!predictiveBool] + targetNonDataPP <- targetNonDataComponents[ predictiveBool] ## predictive + targetNonDataNP <- targetNonDataComponents[!predictiveBool] ## not predictive ## nested function and function list definitions samplerList <- nimbleFunctionList(sampler_BASE) - if(multivariateNodesAsScalars) { - for(i in seq_along(targetUnobservedComponents)) { - samplerList[[i]] <- sampler_RW(model, mvSaved, targetUnobservedComponents[i], control) - } - } else { - if(length(targetUnobservedComponents) == 1) { - samplerList[[1]] <- sampler_RW(model, mvSaved, targetUnobservedComponents, control) + if(length(targetNonDataNP) > 0) { + if(multivariateNodesAsScalars) { + for(i in seq_along(targetNonDataNP)) { + samplerList[[i]] <- sampler_RW(model, mvSaved, targetNonDataNP[i], control) + } } else { - samplerList[[1]] <- sampler_RW_block(model, mvSaved, targetUnobservedComponents, control) + if(length(targetNonDataNP) == 1) { + samplerList[[1]] <- sampler_RW(model, mvSaved, targetNonDataNP, control) + } else { + samplerList[[1]] <- sampler_RW_block(model, mvSaved, targetNonDataNP, control) + } } } + if(length(targetNonDataPP) > 0) { + samplerList[[ length(samplerList)+1 ]] <- sampler_partial_mvn_pp(model, mvSaved, targetNonDataPP) + } ## checks - if (model$getDistribution(target) != "dmnorm") stop(paste0('The node ', target, ' is parially observed. NIMBLE only handles this case for multivariate normal distibutions.')) - if (!model$isMixedData(target)) stop(paste0('The target node ', target, ' is not partially observed.')) + if(model$getDistribution(target) != 'dmnorm') stop(paste0('The node ', target, ' is parially observed. NIMBLE only handles this case for multivariate normal distibutions.')) + if(!model$isMixedData(target)) stop(paste0('The target node ', target, ' is not partially observed.')) }, run = function() { - for (i in seq_along(samplerList)) { + for(i in seq_along(samplerList)) { samplerList[[i]]$run() } }, methods = list( reset = function() { - for (i in seq_along(samplerList)) { + for(i in seq_along(samplerList)) { samplerList[[i]]$reset() } } ) ) +sampler_partial_mvn_pp <- nimbleFunction( + name = 'sampler_partial_mvn_pp', + contains = sampler_BASE, + setup = function(model, mvSaved, target) { + browser() + ## node list generation + mvNode <- model$expandNodeNames(target) + mvNodeComponents <- model$expandNodeNames(mvNode, returnScalarComponents = TRUE) + given <- setdiff(mvNodeComponents, target) + calcNodes <- model$getDependencies(target, downstream = TRUE, includePredictive = TRUE) + ## numeric value generation + nTarget <- length(target) + nGiven <- length(given) + nTotal <- nTarget + nGiven + mu <- numeric(nTotal ) + muT <- numeric(nTarget) + muG <- numeric(nGiven ) + sigma <- array(0, c(nTotal, nTotal )) + sigmaTT <- array(0, c(nTarget, nTarget)) + sigmaTG <- array(0, c(nTarget, nGiven )) + sigmaGT <- array(0, c(nGiven, nTarget)) + sigmaGG <- array(0, c(nGiven, nGiven )) + sigmaPP <- array(0, c(nTarget, nTarget)) + indT <- match(target, mvNodeComponents) + indG <- match(given, mvNodeComponents) + ## checks + if(length(mvNode) != 1) stop('something went wrong in sampler_partial_mvn_pp setup') + if(nTarget*nGiven == 0) stop('something went wrong in sampler_partial_mvn_pp setup') + }, + run = function() { + mu <<- model$getParam(mvNode, 'mean') + muT <<- mu[indT] + muG <<- mu[indG] + sigma <<- model$getParam(mvNode, 'cov') + sigmaTT <<- sigma[indT, indT] + sigmaTG <<- sigma[indT, indG] + sigmaGT <<- sigma[indG, indT] + sigmaGG <<- sigma[indG, indG] + sigmaTG <<- sigmaTG %*% solve(sigmaGG) + muT <<- muT + sigmaTG %*% (values(model,given) - muG) + sigmaTT <<- sigmaTT - sigmaTG %*% sigmaGT + sigmaTT <<- chol(sigmaTT) + values(model, target) <<- rmnorm_chol(1, muT, sigmaTT, prec_param = 0) + model$calculate(calcNodes) + nimCopy(from = model, to = mvSaved, row = 1, nodes = calcNodes, logProb = TRUE) + }, + methods = list( + reset = function() { } + ) +) + + #' MCMC Sampling Algorithms #' #' Details of the MCMC sampling algorithms provided with the NIMBLE MCMC engine; HMC samplers are in the \code{nimbleHMC} package and particle filter samplers are in the \code{nimbleSMC} package. Additional details, including some recommendations for samplers that may perform better than the samplers that NIMBLE assigns by default are provided in Section 7.11 of the User Manual. From 73790990926bd94a7e6f6edd566fec2aef52e3bd Mon Sep 17 00:00:00 2001 From: Daniel Turek Date: Wed, 11 Feb 2026 05:58:00 -0500 Subject: [PATCH 07/16] updates --- packages/nimble/R/MCMC_samplers.R | 31 +++++++++++++++---------------- 1 file changed, 15 insertions(+), 16 deletions(-) diff --git a/packages/nimble/R/MCMC_samplers.R b/packages/nimble/R/MCMC_samplers.R index 707f65c67..f854fc7f3 100644 --- a/packages/nimble/R/MCMC_samplers.R +++ b/packages/nimble/R/MCMC_samplers.R @@ -3528,12 +3528,11 @@ sampler_partial_mvn_pp <- nimbleFunction( mu <- numeric(nTotal ) muT <- numeric(nTarget) muG <- numeric(nGiven ) - sigma <- array(0, c(nTotal, nTotal )) - sigmaTT <- array(0, c(nTarget, nTarget)) - sigmaTG <- array(0, c(nTarget, nGiven )) - sigmaGT <- array(0, c(nGiven, nTarget)) - sigmaGG <- array(0, c(nGiven, nGiven )) - sigmaPP <- array(0, c(nTarget, nTarget)) + Sigma <- array(0, c(nTotal, nTotal )) + SigmaTT <- array(0, c(nTarget, nTarget)) + SigmaTG <- array(0, c(nTarget, nGiven )) + SigmaGT <- array(0, c(nGiven, nTarget)) + SigmaGG <- array(0, c(nGiven, nGiven )) indT <- match(target, mvNodeComponents) indG <- match(given, mvNodeComponents) ## checks @@ -3544,16 +3543,16 @@ sampler_partial_mvn_pp <- nimbleFunction( mu <<- model$getParam(mvNode, 'mean') muT <<- mu[indT] muG <<- mu[indG] - sigma <<- model$getParam(mvNode, 'cov') - sigmaTT <<- sigma[indT, indT] - sigmaTG <<- sigma[indT, indG] - sigmaGT <<- sigma[indG, indT] - sigmaGG <<- sigma[indG, indG] - sigmaTG <<- sigmaTG %*% solve(sigmaGG) - muT <<- muT + sigmaTG %*% (values(model,given) - muG) - sigmaTT <<- sigmaTT - sigmaTG %*% sigmaGT - sigmaTT <<- chol(sigmaTT) - values(model, target) <<- rmnorm_chol(1, muT, sigmaTT, prec_param = 0) + Sigma <<- model$getParam(mvNode, 'cov') + SigmaTT <<- Sigma[indT, indT] + SigmaTG <<- Sigma[indT, indG] + SigmaGT <<- Sigma[indG, indT] + SigmaGG <<- Sigma[indG, indG] + SigmaTG <<- SigmaTG %*% inverse(SigmaGG) + muT <<- muT + (SigmaTG %*% (values(model,given) - muG))[,1] + SigmaTT <<- SigmaTT - SigmaTG %*% SigmaGT + SigmaTT <<- chol(SigmaTT) + values(model, target) <<- rmnorm_chol(1, muT, SigmaTT, prec_param = 0) model$calculate(calcNodes) nimCopy(from = model, to = mvSaved, row = 1, nodes = calcNodes, logProb = TRUE) }, From c53f181569d9a8fc3f425f99304a064a82cbb632 Mon Sep 17 00:00:00 2001 From: Daniel Turek Date: Wed, 11 Feb 2026 11:12:55 -0500 Subject: [PATCH 08/16] updates --- packages/nimble/R/MCMC_samplers.R | 59 ++++++++++++++++--------------- 1 file changed, 30 insertions(+), 29 deletions(-) diff --git a/packages/nimble/R/MCMC_samplers.R b/packages/nimble/R/MCMC_samplers.R index f854fc7f3..fef3c75ee 100644 --- a/packages/nimble/R/MCMC_samplers.R +++ b/packages/nimble/R/MCMC_samplers.R @@ -3515,44 +3515,45 @@ sampler_partial_mvn_pp <- nimbleFunction( name = 'sampler_partial_mvn_pp', contains = sampler_BASE, setup = function(model, mvSaved, target) { - browser() ## node list generation mvNode <- model$expandNodeNames(target) mvNodeComponents <- model$expandNodeNames(mvNode, returnScalarComponents = TRUE) given <- setdiff(mvNodeComponents, target) calcNodes <- model$getDependencies(target, downstream = TRUE, includePredictive = TRUE) ## numeric value generation - nTarget <- length(target) - nGiven <- length(given) - nTotal <- nTarget + nGiven - mu <- numeric(nTotal ) - muT <- numeric(nTarget) - muG <- numeric(nGiven ) - Sigma <- array(0, c(nTotal, nTotal )) - SigmaTT <- array(0, c(nTarget, nTarget)) - SigmaTG <- array(0, c(nTarget, nGiven )) - SigmaGT <- array(0, c(nGiven, nTarget)) - SigmaGG <- array(0, c(nGiven, nGiven )) - indT <- match(target, mvNodeComponents) - indG <- match(given, mvNodeComponents) + n1 <- length(target) + n2 <- length(given) + n <- n1 + n2 + mu <- array(0, c(n , 1)) + mu1 <- array(0, c(n1, 1)) + mu2 <- array(0, c(n2, 1)) + Sigma <- array(0, c(n, n )) + Sigma11 <- array(0, c(n1, n1)) + Sigma12 <- array(0, c(n1, n2)) + Sigma21 <- array(0, c(n2, n1)) + Sigma22 <- array(0, c(n2, n2)) + ind1 <- match(target, mvNodeComponents) + ind2 <- match(given, mvNodeComponents) ## checks - if(length(mvNode) != 1) stop('something went wrong in sampler_partial_mvn_pp setup') - if(nTarget*nGiven == 0) stop('something went wrong in sampler_partial_mvn_pp setup') + if(length(mvNode) != 1) stop('something went wrong in sampler_partial_mvn_pp setup') + if(model$getDistribution(mvNode) != 'dmnorm') stop('something went wrong in sampler_partial_mvn_pp setup') + if(n != length(mvNodeComponents)) stop('something went wrong in sampler_partial_mvn_pp setup') + if(n1*n2 == 0) stop('something went wrong in sampler_partial_mvn_pp setup') }, run = function() { - mu <<- model$getParam(mvNode, 'mean') - muT <<- mu[indT] - muG <<- mu[indG] - Sigma <<- model$getParam(mvNode, 'cov') - SigmaTT <<- Sigma[indT, indT] - SigmaTG <<- Sigma[indT, indG] - SigmaGT <<- Sigma[indG, indT] - SigmaGG <<- Sigma[indG, indG] - SigmaTG <<- SigmaTG %*% inverse(SigmaGG) - muT <<- muT + (SigmaTG %*% (values(model,given) - muG))[,1] - SigmaTT <<- SigmaTT - SigmaTG %*% SigmaGT - SigmaTT <<- chol(SigmaTT) - values(model, target) <<- rmnorm_chol(1, muT, SigmaTT, prec_param = 0) + mu[,1] <<- model$getParam(mvNode, 'mean') + Sigma <<- model$getParam(mvNode, 'cov' ) + mu1[,1] <<- mu[ind1,1] + mu2[,1] <<- mu[ind2,1] + Sigma11[1:n1,1:n1] <<- Sigma[ind1, ind1] + Sigma12[1:n1,1:n2] <<- Sigma[ind1, ind2] + Sigma21[1:n2,1:n1] <<- Sigma[ind2, ind1] + Sigma22[1:n2,1:n2] <<- Sigma[ind2, ind2] + Sigma12 <<- Sigma12 %*% inverse(Sigma22) + mu1[,1] <<- mu1[,1] + (Sigma12 %*% (values(model,given) - mu2[,1]))[,1] + Sigma11 <<- Sigma11 - Sigma12 %*% Sigma21 + Sigma11 <<- chol(Sigma11) + values(model, target) <<- rmnorm_chol(1, mu1[,1], Sigma11, prec_param = 0) model$calculate(calcNodes) nimCopy(from = model, to = mvSaved, row = 1, nodes = calcNodes, logProb = TRUE) }, From 84653645b0a872ba8eea8589fcff80125deac0bf Mon Sep 17 00:00:00 2001 From: Daniel Turek Date: Wed, 11 Feb 2026 15:41:37 -0500 Subject: [PATCH 09/16] removed extraneous tests of partial_mvn sampler --- packages/nimble/tests/testthat/test-mcmc.R | 89 ---------------------- 1 file changed, 89 deletions(-) diff --git a/packages/nimble/tests/testthat/test-mcmc.R b/packages/nimble/tests/testthat/test-mcmc.R index ec61c049a..bfee35ba0 100644 --- a/packages/nimble/tests/testthat/test-mcmc.R +++ b/packages/nimble/tests/testthat/test-mcmc.R @@ -3241,36 +3241,6 @@ test_that('partial_mvn sampler was not given to non dmnorm distribution from mod expect_true(!any(sapply(conf$getSamplers(), function(sc) sc$name)[head(seq_along(sapply(conf$getSamplers(), function(sc) sc$name)), -1)]=='partial_mvn'), info="partial_mvn sampler assigned as sampler to non dmnorm node in model with dmnorm") }) -test_that('Values change at different rows for unobserved parts of partially observed node when multivariateNodesAsScalars = TRUE', { - Code <- nimbleCode({ for (i in 1:N){ - theta[i] ~ dgamma(alpha,beta) - lambda[i] <- theta[i]*t[i] - x[i] ~ dpois(lambda[i]) - } - alpha ~ dexp(1.0) - beta ~ dgamma(0.1,1.0) - y[1:5]~dmnorm(mu[1:5], Sigma[1:5,1:5]) - }) - Consts <- list(N = 10, t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5)) - Data <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22), y=c(rep(NA,3),4,5)) - Inits <- list(alpha = 1, beta = 1, theta = rep(0.1, Consts$N), mu=1:5, Sigma=diag(5), y=c(rep(1,3),rep(NA,2))) - model <- nimbleModel(code = Code, name = "model", constants = Consts, data = Data, inits = Inits) - - Cmodel <- compileNimble(model) - - conf <- configureMCMC(model, nodes = 'y[1:5]', multivariateNodesAsScalars = TRUE) - conf$addMonitors("y") - - modelMCMC <- buildMCMC(conf) - CmodelMCMC <- compileNimble(modelMCMC, project = model) - modrun <- runMCMC(CmodelMCMC, niter=100, setSeed = 0) - - - onezro=apply(modrun[,3:5], 2, diff)!=0 - - expect_true(any(apply(onezro, 1, function(x) any(x) && !all(x))), info= "Values of MCMC sampled nodes change at same time as other nodes") -}) - test_that('Values change at same rows for unobserved parts of partially observed node when multivariateNodesAsScalars = FALSE', { Code <- nimbleCode({ for (i in 1:N){ theta[i] ~ dgamma(alpha,beta) @@ -3301,65 +3271,6 @@ test_that('Values change at same rows for unobserved parts of partially observed }) -test_that('Values do not change for data when multivariateNodesAsScalars = TRUE', { - Code <- nimbleCode({ for (i in 1:N){ - theta[i] ~ dgamma(alpha,beta) - lambda[i] <- theta[i]*t[i] - x[i] ~ dpois(lambda[i]) - } - alpha ~ dexp(1.0) - beta ~ dgamma(0.1,1.0) - y[1:5]~dmnorm(mu[1:5], Sigma[1:5,1:5]) - }) - Consts <- list(N = 10, t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5)) - Data <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22), y=c(rep(NA,3),4,5)) - Inits <- list(alpha = 1, beta = 1, theta = rep(0.1, Consts$N), mu=1:5, Sigma=diag(5), y=c(rep(1,3),rep(NA,2))) - model <- nimbleModel(code = Code, name = "model", constants = Consts, data = Data, inits = Inits) - - Cmodel <- compileNimble(model) - - conf <- configureMCMC(model, nodes = 'y[1:5]', multivariateNodesAsScalars = TRUE) - conf$addMonitors("y") - - modelMCMC <- buildMCMC(conf) - CmodelMCMC <- compileNimble(modelMCMC, project = model) - modrun <- runMCMC(CmodelMCMC, niter=100, setSeed = 0) - - onezro <- apply(modrun[,c(1,2,6,7)], 2, diff)!=0 - - expect_true(all(apply(onezro, 1, function(x) !any(x))), info= "Values of data do not change when MCMC sampling") -}) - - -test_that('Values do not change for data when multivariateNodesAsScalars = FALSE', { - Code <- nimbleCode({ for (i in 1:N){ - theta[i] ~ dgamma(alpha,beta) - lambda[i] <- theta[i]*t[i] - x[i] ~ dpois(lambda[i]) - } - alpha ~ dexp(1.0) - beta ~ dgamma(0.1,1.0) - y[1:5]~dmnorm(mu[1:5], Sigma[1:5,1:5]) - }) - Consts <- list(N = 10, t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5)) - Data <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22), y=c(rep(NA,3),4,5)) - Inits <- list(alpha = 1, beta = 1, theta = rep(0.1, Consts$N), mu=1:5, Sigma=diag(5), y=c(rep(1,3),rep(NA,2))) - model <- nimbleModel(code = Code, name = "model", constants = Consts, data = Data, inits = Inits) - - Cmodel <- compileNimble(model) - - conf <- configureMCMC(model, nodes = 'y[1:5]', multivariateNodesAsScalars = FALSE) - conf$addMonitors("y") - - modelMCMC <- buildMCMC(conf) - CmodelMCMC <- compileNimble(modelMCMC, project = model) - modrun <- runMCMC(CmodelMCMC, niter=100, setSeed = 0) - - onezro <- apply(modrun[,c(1,2,6,7)], 2, diff)!=0 - - expect_true(all(apply(onezro, 1, function(x) !any(x))), info= "Values of data do not change when MCMC sampling") -}) - test_that('asymptotic correct results from conjugate gamma - CAR_normal sampler', { num1 <- c(1, 2, 2, 1) adj1 <- c(2, 1, 3, 2, 4, 3) From 84474201e264f5530ed651e7b7b05d403d9914b9 Mon Sep 17 00:00:00 2001 From: Daniel Turek Date: Wed, 11 Feb 2026 16:17:35 -0500 Subject: [PATCH 10/16] updates per comments on PR --- packages/nimble/R/MCMC_samplers.R | 27 ++++++++++++++++----------- 1 file changed, 16 insertions(+), 11 deletions(-) diff --git a/packages/nimble/R/MCMC_samplers.R b/packages/nimble/R/MCMC_samplers.R index fef3c75ee..bb44d7527 100644 --- a/packages/nimble/R/MCMC_samplers.R +++ b/packages/nimble/R/MCMC_samplers.R @@ -3470,11 +3470,16 @@ sampler_partial_mvn <- nimbleFunction( multivariateNodesAsScalars <- extractControlElement(control, 'multivariateNodesAsScalars', default = getNimbleOption('MCMCmultivariateNodesAsScalars')) ## node list generation targetAsScalar <- model$expandNodeNames(target, returnScalarComponents = TRUE) - isdataBool <- sapply(targetAsScalar, function(targetAsScalar) eval(parse(text = targetAsScalar)[[1]], envir = model$isDataEnv)) - targetNonDataComponents <- targetAsScalar[!isdataBool] - predictiveBool <- sapply(targetNonDataComponents, function(n) length(model$getDependencies(n, self=FALSE, downstream=TRUE, dataOnly=TRUE)) == 0) - targetNonDataPP <- targetNonDataComponents[ predictiveBool] ## predictive - targetNonDataNP <- targetNonDataComponents[!predictiveBool] ## not predictive + isDataBool <- sapply(targetAsScalar, function(targetAsScalar) eval(parse(text = targetAsScalar)[[1]], envir = model$isDataEnv)) + targetNonDataComponents <- targetAsScalar[!isDataBool] + if(length(model$getDependencies(target, self = FALSE, downstream = TRUE, dataOnly = TRUE)) > 0) { + predictiveBool <- sapply(targetNonDataComponents, function(n) length(model$getDependencies(n, self = FALSE, downstream = TRUE, dataOnly = TRUE)) == 0) + targetNonDataPP <- targetNonDataComponents[ predictiveBool] ## predictive + targetNonDataNP <- targetNonDataComponents[!predictiveBool] ## not predictive + } else { + targetNonDataPP <- targetNonDataComponents ## entirely predictive + targetNonDataNP <- character() + } ## nested function and function list definitions samplerList <- nimbleFunctionList(sampler_BASE) if(length(targetNonDataNP) > 0) { @@ -3494,8 +3499,8 @@ sampler_partial_mvn <- nimbleFunction( samplerList[[ length(samplerList)+1 ]] <- sampler_partial_mvn_pp(model, mvSaved, targetNonDataPP) } ## checks - if(model$getDistribution(target) != 'dmnorm') stop(paste0('The node ', target, ' is parially observed. NIMBLE only handles this case for multivariate normal distibutions.')) - if(!model$isMixedData(target)) stop(paste0('The target node ', target, ' is not partially observed.')) + if(model$getDistribution(target) != 'dmnorm') stop('The node ', target, ' is parially observed. NIMBLE only handles this case for multivariate normal distibutions.') + if(!model$isMixedData(target)) stop('The target node ', target, ' is not partially observed.') }, run = function() { for(i in seq_along(samplerList)) { @@ -3535,10 +3540,10 @@ sampler_partial_mvn_pp <- nimbleFunction( ind1 <- match(target, mvNodeComponents) ind2 <- match(given, mvNodeComponents) ## checks - if(length(mvNode) != 1) stop('something went wrong in sampler_partial_mvn_pp setup') - if(model$getDistribution(mvNode) != 'dmnorm') stop('something went wrong in sampler_partial_mvn_pp setup') - if(n != length(mvNodeComponents)) stop('something went wrong in sampler_partial_mvn_pp setup') - if(n1*n2 == 0) stop('something went wrong in sampler_partial_mvn_pp setup') + if(length(mvNode) != 1) stop('unexpected error in sampler_partial_mvn_pp') + if(model$getDistribution(mvNode) != 'dmnorm') stop('unexpected error in sampler_partial_mvn_pp') + if(n != length(mvNodeComponents)) stop('unexpected error in sampler_partial_mvn_pp') + if(n1*n2 == 0) stop('unexpected error in sampler_partial_mvn_pp') }, run = function() { mu[,1] <<- model$getParam(mvNode, 'mean') From 2a4f5812efe8719ee6eda8cd6b898269e0e90554 Mon Sep 17 00:00:00 2001 From: Daniel Turek Date: Wed, 11 Feb 2026 16:36:00 -0500 Subject: [PATCH 11/16] partial_mvn sampler can assign barker sampler --- packages/nimble/R/MCMC_samplers.R | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/packages/nimble/R/MCMC_samplers.R b/packages/nimble/R/MCMC_samplers.R index bb44d7527..aea8e3e02 100644 --- a/packages/nimble/R/MCMC_samplers.R +++ b/packages/nimble/R/MCMC_samplers.R @@ -3491,7 +3491,11 @@ sampler_partial_mvn <- nimbleFunction( if(length(targetNonDataNP) == 1) { samplerList[[1]] <- sampler_RW(model, mvSaved, targetNonDataNP, control) } else { - samplerList[[1]] <- sampler_RW_block(model, mvSaved, targetNonDataNP, control) + if(getNimbleOption('MCMCuseBarkerAsDefaultMV')) { + samplerList[[1]] <- sampler_barker (model, mvSaved, targetNonDataNP, control) + } else { + samplerList[[1]] <- sampler_RW_block(model, mvSaved, targetNonDataNP, control) + } } } } From bc50d56fa1b6ee15c5e6b8d48278e50729dd0b2c Mon Sep 17 00:00:00 2001 From: Christopher Paciorek Date: Wed, 11 Feb 2026 13:43:08 -0800 Subject: [PATCH 12/16] Use Cholesky of Sigma22 in partial_mvn_pp. --- packages/nimble/R/MCMC_samplers.R | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/packages/nimble/R/MCMC_samplers.R b/packages/nimble/R/MCMC_samplers.R index bb44d7527..bdd4fc0fc 100644 --- a/packages/nimble/R/MCMC_samplers.R +++ b/packages/nimble/R/MCMC_samplers.R @@ -3537,6 +3537,8 @@ sampler_partial_mvn_pp <- nimbleFunction( Sigma12 <- array(0, c(n1, n2)) Sigma21 <- array(0, c(n2, n1)) Sigma22 <- array(0, c(n2, n2)) + tmp <- array(0, c(n1, 1)) + tmp2 <- array(0, c(n2, n1)) ind1 <- match(target, mvNodeComponents) ind2 <- match(given, mvNodeComponents) ## checks @@ -3554,9 +3556,14 @@ sampler_partial_mvn_pp <- nimbleFunction( Sigma12[1:n1,1:n2] <<- Sigma[ind1, ind2] Sigma21[1:n2,1:n1] <<- Sigma[ind2, ind1] Sigma22[1:n2,1:n2] <<- Sigma[ind2, ind2] - Sigma12 <<- Sigma12 %*% inverse(Sigma22) - mu1[,1] <<- mu1[,1] + (Sigma12 %*% (values(model,given) - mu2[,1]))[,1] - Sigma11 <<- Sigma11 - Sigma12 %*% Sigma21 + Sigma22 <<- chol(Sigma22) + ## Sigma12 <<- Sigma12 %*% inverse(Sigma22) + ## mu1[,1] <<- mu1[,1] + (Sigma12 %*% (values(model,given) - mu2[,1]))[,1] + tmp[,1] <<- backsolve(Sigma22, forwardsolve(t(Sigma22), values(model,given) - mu2[,1])) + mu1[,1] <<- mu1[,1] + (Sigma12 %*% tmp)[,1] + ## Sigma11 <<- Sigma11 - Sigma12 %*% Sigma21 + tmp2 <<- forwardsolve(t(Sigma22), Sigma21) + Sigma11 <<- Sigma11 - t(tmp2)%*%tmp2 Sigma11 <<- chol(Sigma11) values(model, target) <<- rmnorm_chol(1, mu1[,1], Sigma11, prec_param = 0) model$calculate(calcNodes) From b5e1fa89fd2ae733e481ba48b8a9a6a45b1e8e82 Mon Sep 17 00:00:00 2001 From: Daniel Turek Date: Wed, 11 Feb 2026 20:40:52 -0500 Subject: [PATCH 13/16] fixed dimension mistake --- packages/nimble/R/MCMC_samplers.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/packages/nimble/R/MCMC_samplers.R b/packages/nimble/R/MCMC_samplers.R index 5aa219ef5..6ada3f327 100644 --- a/packages/nimble/R/MCMC_samplers.R +++ b/packages/nimble/R/MCMC_samplers.R @@ -3541,7 +3541,7 @@ sampler_partial_mvn_pp <- nimbleFunction( Sigma12 <- array(0, c(n1, n2)) Sigma21 <- array(0, c(n2, n1)) Sigma22 <- array(0, c(n2, n2)) - tmp <- array(0, c(n1, 1)) + tmp <- array(0, c(n2, 1)) tmp2 <- array(0, c(n2, n1)) ind1 <- match(target, mvNodeComponents) ind2 <- match(given, mvNodeComponents) @@ -3567,7 +3567,7 @@ sampler_partial_mvn_pp <- nimbleFunction( mu1[,1] <<- mu1[,1] + (Sigma12 %*% tmp)[,1] ## Sigma11 <<- Sigma11 - Sigma12 %*% Sigma21 tmp2 <<- forwardsolve(t(Sigma22), Sigma21) - Sigma11 <<- Sigma11 - t(tmp2)%*%tmp2 + Sigma11 <<- Sigma11 - t(tmp2) %*% tmp2 Sigma11 <<- chol(Sigma11) values(model, target) <<- rmnorm_chol(1, mu1[,1], Sigma11, prec_param = 0) model$calculate(calcNodes) From 10c86032861ac2a35888366cba31cc41c829d83c Mon Sep 17 00:00:00 2001 From: Christopher Paciorek Date: Wed, 11 Feb 2026 17:57:16 -0800 Subject: [PATCH 14/16] Use Cholesky in sampler_partial_mvn_pp. --- packages/nimble/R/MCMC_samplers.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/packages/nimble/R/MCMC_samplers.R b/packages/nimble/R/MCMC_samplers.R index 5aa219ef5..1efc70ce4 100644 --- a/packages/nimble/R/MCMC_samplers.R +++ b/packages/nimble/R/MCMC_samplers.R @@ -3541,7 +3541,7 @@ sampler_partial_mvn_pp <- nimbleFunction( Sigma12 <- array(0, c(n1, n2)) Sigma21 <- array(0, c(n2, n1)) Sigma22 <- array(0, c(n2, n2)) - tmp <- array(0, c(n1, 1)) + tmp <- array(0, c(n2, 1)) tmp2 <- array(0, c(n2, n1)) ind1 <- match(target, mvNodeComponents) ind2 <- match(given, mvNodeComponents) From a2ffffb30cafa873fa47ae721aeaf7ad9b761405 Mon Sep 17 00:00:00 2001 From: Christopher Paciorek Date: Wed, 11 Feb 2026 17:57:31 -0800 Subject: [PATCH 15/16] Add more testing for sampler_partial_mvn. --- packages/nimble/tests/testthat/test-mcmc.R | 92 +++++++++++++++++++++- 1 file changed, 89 insertions(+), 3 deletions(-) diff --git a/packages/nimble/tests/testthat/test-mcmc.R b/packages/nimble/tests/testthat/test-mcmc.R index bfee35ba0..01988b416 100644 --- a/packages/nimble/tests/testthat/test-mcmc.R +++ b/packages/nimble/tests/testthat/test-mcmc.R @@ -3041,10 +3041,47 @@ test_that('partial_mvn sampler was given to dmnorm distribution when dependent o model <- nimbleModel(code = code, name = "model", constants = Consts, data = Data) conf <- configureMCMC(model, nodes = 'y[1:5]') - - expect_true(any(sapply(conf$getSamplers(), function(sc) sc$name)=='partial_mvn', info= "partial_mvn sampler not assigned to dmnorm distribution")) + mcmc <- buildMCMC(conf) + + + expect_true(any(sapply(conf$getSamplers(), function(sc) sc$name)=='partial_mvn', + info= "partial_mvn sampler not assigned to dmnorm distribution")) + expect_true(inherits(mcmc$samplerFunctions[[1]]$samplerList[[1]], "sampler_partial_mvn_pp")) + expect_identical(mcmc$samplerFunctions[[1]]$samplerList[[1]]$target, c("y[1]","y[2]","y[3]")) }) +test_that('partial_mvn sampler was given to dmnorm distribution when node is intermediate', { + code <- nimbleCode({ for (i in 1:N){ + theta[i] ~ dgamma(alpha,beta) + lambda[i] <- theta[i]*t[i] + x[i] ~ dpois(lambda[i]) + } + alpha ~ dexp(1.0) + beta ~ dgamma(0.1,1.0) + for(i in 1:10) { + mu[i] <- alpha+i + } + y[1:10]~dmnorm(mu[1:10], Sigma[1:10,1:10]) + for(i in 1:3) + z[i] ~ dnorm(y[i], 1) + }) + Consts <- list(N = 10, t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5)) + Data <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22), y=c(rep(NA,5),rep(3,5)), z = rep(1,3)) + + model <- nimbleModel(code = code, name = "model", constants = Consts, data = Data) + + conf <- configureMCMC(model, nodes = 'y[1:5]') + mcmc <- buildMCMC(conf) + + + expect_true(any(sapply(conf$getSamplers(), function(sc) sc$name)=='partial_mvn', + info= "partial_mvn sampler not assigned to dmnorm distribution")) + + expect_true(inherits(mcmc$samplerFunctions[[1]]$samplerList[[1]], "sampler_RW_block")) + expect_identical(mcmc$samplerFunctions[[1]]$samplerList[[1]]$target, c("y[1]","y[2]","y[3]")) + expect_true(inherits(mcmc$samplerFunctions[[1]]$samplerList[[2]], "sampler_partial_mvn_pp")) + expect_identical(mcmc$samplerFunctions[[1]]$samplerList[[2]]$target, c("y[4]", "y[5]")) +}) test_that('partial_mvn sampler was given to dmnorm distribution when not dependent on nodes in the model', { code <- nimbleCode({ for (i in 1:N){ @@ -3062,8 +3099,10 @@ test_that('partial_mvn sampler was given to dmnorm distribution when not depende model <- nimbleModel(code = code, name = "model", constants = consts, data = data) conf <- configureMCMC(model, nodes = 'y[1:10]') + mcmc <- buildMCMC(conf) expect_true(any(sapply(conf$getSamplers(), function(sc) sc$name)=='partial_mvn', info= "partial_mvn sampler not assigned to dmnorm distribution")) + expect_true(inherits(mcmc$samplerFunctions[[1]]$samplerList[[1]], "sampler_partial_mvn_pp")) }) test_that('partial_mvn sampler was given to dmnorm distribution when intermediate in a model', { @@ -3086,8 +3125,10 @@ test_that('partial_mvn sampler was given to dmnorm distribution when intermediat model <- nimbleModel(code = code, name = "model", constants = Consts, data = Data) conf <- configureMCMC(model, nodes = 'y[1:5]') + mcmc <- buildMCMC(conf) expect_true(any(sapply(conf$getSamplers(), function(sc) sc$name)=='partial_mvn', info= "partial_mvn sampler not assigned to dmnorm distribution")) + expect_true(inherits(mcmc$samplerFunctions[[1]]$samplerList[[1]], "sampler_partial_mvn_pp")) }) test_that('partial_mvn sampler was given to each dmnorm distribution', { @@ -3111,8 +3152,11 @@ test_that('partial_mvn sampler was given to each dmnorm distribution', { model <- nimbleModel(code = code, name = "model", constants = Consts, data = Data) conf <- configureMCMC(model, nodes = model$getNodeNames()[26:27]) + mcmc <- buildMCMC(conf) expect_true(sum(sapply(conf$getSamplers(), function(sc) sc$name)=='partial_mvn')==2, info= "partial_mvn sampler not assigned to each dmnorm distribution") + expect_true(inherits(mcmc$samplerFunctions[[1]]$samplerList[[1]], "sampler_partial_mvn_pp")) + expect_true(inherits(mcmc$samplerFunctions[[2]]$samplerList[[1]], "sampler_partial_mvn_pp")) }) test_that('partial_mvn sampler was given to each partially obs dmnorm distribution when both are at intermediate level in a model', { @@ -3139,8 +3183,11 @@ test_that('partial_mvn sampler was given to each partially obs dmnorm distributi model <- nimbleModel(code = code, name = "model", constants = Consts, data = Data) conf <- configureMCMC(model, nodes = model$getNodeNames()[26:27]) + mcmc <- buildMCMC(conf) expect_true(sum(sapply(conf$getSamplers(), function(sc) sc$name)=='partial_mvn')==2, info= "partial_mvn sampler not assigned to each dmnorm distribution") + expect_true(inherits(mcmc$samplerFunctions[[1]]$samplerList[[1]], "sampler_partial_mvn_pp")) + expect_true(inherits(mcmc$samplerFunctions[[2]]$samplerList[[1]], "sampler_partial_mvn_pp")) }) test_that('an error is given when trying to assign sampler to partially observed dmvt dist', { @@ -3184,7 +3231,7 @@ test_that('an error is given when trying to assign a sampler to partially observ expect_error(configureMCMC(model, print = FALSE), info = "model with dmvt given a sampler for every node" ) }) -test_that('trying to give a sampler to an observed node doesnt yield an error from configureMCMC', { +test_that('trying to give a sampler to an observed node does not yield an error from configureMCMC', { code <- nimbleCode({ for (i in 1:N){ theta[i] ~ dgamma(alpha,beta) lambda[i] <- theta[i]*t[i] @@ -3270,6 +3317,45 @@ test_that('Values change at same rows for unobserved parts of partially observed expect_true(all(apply(onezro, 1, function(x) all(x) || all(!x))), info= "Values of MCMC sampled nodes change at same time as other nodes") }) +test_that('partial_mvn_pp sampler with scalar cases', { + code <- nimbleCode({ + alpha ~ dexp(1.0) + for(i in 1:5) { + mu[i] <- alpha+i + } + y[1:5]~dmnorm(mu[1:5], Sigma[1:5,1:5]) + }) + ## Sample univariate + Data <- list(y = c(NA, rep(0,4))) + Const <- list(alpha = 1, Sigma = diag(5)) + model <- nimbleModel(code = code, name = "model", data = Data, constants = Const, inits = list(y=c(0, rep(1,4)))) + + conf <- configureMCMC(model, nodes = 'y[1:5]', monitors = 'y') + mcmc <- buildMCMC(conf) + + cmodel <- compileNimble(model) + cmcmc <- compileNimble(mcmc, project = model) + + expect_success(out <- runMCMC(mcmc, 5)) + expect_success(cout <- runMCMC(cmcmc, 5)) + + ## Condition on univariate + Data <- list(y = c(rep(NA, 4), 0)) + model <- nimbleModel(code = code, name = "model", data = Data, constants = Const, inits = list(y=rep(1,5))) + + ## Sample univariate + conf <- configureMCMC(model, nodes = 'y[1:5]', monitors = 'y') + mcmc <- buildMCMC(conf) + + cmodel <- compileNimble(model) + cmcmc <- compileNimble(mcmc, project = model) + + expect_success(out <- runMCMC(mcmc, 5)) + expect_success(cout <- runMCMC(cmcmc, 5)) + +}) + + test_that('asymptotic correct results from conjugate gamma - CAR_normal sampler', { num1 <- c(1, 2, 2, 1) From 74d0e65481277e69c3bb2f0474e6b15e92e29132 Mon Sep 17 00:00:00 2001 From: Daniel Turek Date: Wed, 11 Feb 2026 21:05:21 -0500 Subject: [PATCH 16/16] fix spacing again --- packages/nimble/R/MCMC_samplers.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/packages/nimble/R/MCMC_samplers.R b/packages/nimble/R/MCMC_samplers.R index df75616d6..6ada3f327 100644 --- a/packages/nimble/R/MCMC_samplers.R +++ b/packages/nimble/R/MCMC_samplers.R @@ -3541,7 +3541,7 @@ sampler_partial_mvn_pp <- nimbleFunction( Sigma12 <- array(0, c(n1, n2)) Sigma21 <- array(0, c(n2, n1)) Sigma22 <- array(0, c(n2, n2)) - tmp <- array(0, c(n2, 1)) + tmp <- array(0, c(n2, 1)) tmp2 <- array(0, c(n2, n1)) ind1 <- match(target, mvNodeComponents) ind2 <- match(given, mvNodeComponents)