diff --git a/packages/RELEASE_INSTRUCTIONS b/packages/RELEASE_INSTRUCTIONS index dd26a136b..8ad20f6b5 100644 --- a/packages/RELEASE_INSTRUCTIONS +++ b/packages/RELEASE_INSTRUCTIONS @@ -85,11 +85,13 @@ R CMD check --as-cran nimble_${VERSION}.tar.gz Officially this should be done on a development build of R, so this is just a first pass and we check with r-devel via winbuilder or r-hub below. -Recently I've been getting this warning on my Ubuntu machine: +Recently I've been getting warnings like this on my Ubuntu machine: ``` Compilation used the following non-portable flag(s): ‘-Wdate-time’ ‘-Werror=format-security’ ‘-Wformat’ +Compilation used the following non-portable flag(s): + ‘-mno-omit-leaf-frame-pointer’ ``` Seems like it might just be my local system as it doesn't show up on winbuilder or r-hub checking. diff --git a/packages/nimble/R/BUGS_model.R b/packages/nimble/R/BUGS_model.R index b4269207b..2cf6e609b 100644 --- a/packages/nimble/R/BUGS_model.R +++ b/packages/nimble/R/BUGS_model.R @@ -499,6 +499,22 @@ 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))) + 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}, 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]) + 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)) { @@ -787,15 +803,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/R/MCMC_configuration.R b/packages/nimble/R/MCMC_configuration.R index 3a91dc2d6..e663e7b49 100644 --- a/packages/nimble/R/MCMC_configuration.R +++ b/packages/nimble/R/MCMC_configuration.R @@ -245,6 +245,11 @@ print: A logical argument specifying whether to print the monitors and samplers. if(missing(nodes)) { nodes <- model$getNodeNames(stochOnly = TRUE, includeData = FALSE, includePredictive = samplePredictiveNodes) # 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 @@ -431,6 +436,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)) { @@ -807,7 +819,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 827ac9126..6ada3f327 100644 --- a/packages/nimble/R/MCMC_samplers.R +++ b/packages/nimble/R/MCMC_samplers.R @@ -3456,7 +3456,129 @@ sampler_barker <- 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', 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] + 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) { + if(multivariateNodesAsScalars) { + for(i in seq_along(targetNonDataNP)) { + samplerList[[i]] <- sampler_RW(model, mvSaved, targetNonDataNP[i], control) + } + } else { + if(length(targetNonDataNP) == 1) { + samplerList[[1]] <- sampler_RW(model, mvSaved, targetNonDataNP, control) + } else { + if(getNimbleOption('MCMCuseBarkerAsDefaultMV')) { + samplerList[[1]] <- sampler_barker (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('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)) { + samplerList[[i]]$run() + } + }, + methods = list( + reset = function() { + 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) { + ## 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 + 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)) + tmp <- array(0, c(n2, 1)) + tmp2 <- array(0, c(n2, n1)) + ind1 <- match(target, mvNodeComponents) + ind2 <- match(given, mvNodeComponents) + ## checks + 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') + 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] + 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) + 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. @@ -3881,6 +4003,15 @@ sampler_barker <- 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: #' @@ -3896,7 +4027,7 @@ sampler_barker <- 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 sampler_barker barker +#' @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 sampler_barker barker sampler_partial_mvn 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 46bd3658a..3fc68d378 100644 --- a/packages/nimble/man/samplers.Rd +++ b/packages/nimble/man/samplers.Rd @@ -24,6 +24,7 @@ \alias{sampler_CAR_proper} \alias{sampler_polyagamma} \alias{sampler_barker} +\alias{sampler_partial_mvn} \alias{samplers} \alias{sampler} \alias{binary} @@ -51,6 +52,7 @@ \alias{RW_lkj_corr_cholesky} \alias{RW_block_lkj_corr_cholesky} \alias{barker} +\alias{partial_mvn} \alias{sampler_RJ_fixed_prior} \alias{sampler_RJ_indicator} \alias{sampler_RJ_toggled} @@ -105,6 +107,8 @@ sampler_polyagamma(model, mvSaved, target, control) sampler_barker(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) @@ -595,6 +599,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}{ diff --git a/packages/nimble/tests/testthat/test-mcmc.R b/packages/nimble/tests/testthat/test-mcmc.R index 1f4cf8088..01988b416 100644 --- a/packages/nimble/tests/testthat/test-mcmc.R +++ b/packages/nimble/tests/testthat/test-mcmc.R @@ -79,32 +79,32 @@ test_mcmc('dugongs', numItsC = 1000, resampleData = TRUE) test_mcmc('epil', model = 'epil2.bug', inits = 'epil-inits.R', - data = 'epil-data.R', numItsC = 1000, resampleData = TRUE) + data = 'epil-data.R', numItsC = 1000, resampleData = TRUE) # looks ok test_mcmc('epil', model = 'epil3.bug', inits = 'epil-inits.R', - data = 'epil-data.R', numItsC = 1000, resampleData = TRUE) + data = 'epil-data.R', numItsC = 1000, resampleData = TRUE) # looks ok test_mcmc('seeds', model = 'seedsuni.bug', inits = 'seeds-init.R', - data = 'seeds-data.R', numItsC = 1000, resampleData = TRUE) + data = 'seeds-data.R', numItsC = 1000, resampleData = TRUE) # looks fine - intervals for b's seem a bit large but probably ok # particularly since default seeds.bug seems fine # results compared to JAGS look fine test_mcmc('seeds', model = 'seedssig.bug', inits = 'seeds-init.R', - data = 'seeds-data.R', numItsC = 1000, resampleData = TRUE) + data = 'seeds-data.R', numItsC = 1000, resampleData = TRUE) # looks fine - intervals for b's seem a bit large but probably ok test_mcmc('birats', model = 'birats1.bug', inits = 'birats-inits.R', - data = 'birats-data.R', numItsC = 1000, resampleData = TRUE) + data = 'birats-data.R', numItsC = 1000, resampleData = TRUE) # seems fine test_mcmc('birats', model = 'birats3.bug', inits = 'birats-inits.R', - data = 'birats-data.R', numItsC = 1000, resampleData = TRUE) + data = 'birats-data.R', numItsC = 1000, resampleData = TRUE) # seems fine test_mcmc('birats', model = 'birats2.bug', inits = 'birats-inits.R', - data = 'birats-data.R', numItsC = 1000, resampleData = TRUE) + data = 'birats-data.R', numItsC = 1000, resampleData = TRUE) # looks fine now that values() returns in order # result changes as of v0.4 because in v0.3-1 'omega.beta' was found # as both topNode and nontopNode and was being simulated into @@ -118,22 +118,22 @@ test_mcmc('ice', model = 'icear.bug', inits = 'ice-inits.R', # are simulated from their priors to have large magnitude values test_that('ice example reworked', { - # rework ice example so that beta[1] and beta[2] will be top nodes + # rework ice example so that beta[1] and beta[2] will be top nodes system.in.dir(paste("sed 's/tau\\*1.0E-6/1.0E-6/g' icear.bug > ", file.path(tempdir(), "icear.bug")), dir = system.file('classic-bugs','vol2','ice', package = 'nimble')) test_mcmc(model = file.path(tempdir(), "icear.bug"), inits = system.file('classic-bugs', 'vol2', 'ice','ice-inits.R', package = 'nimble'), data = system.file('classic-bugs', 'vol2', 'ice','ice-data.R', package = 'nimble'), numItsC = 1000, resampleData = TRUE, avoidNestedTest = TRUE) - # looks fine, but alpha and beta values shifted a bit (systematically) relative to JAGS results - on further inspection this is because mixing for this model is poor in both NIMBLE and JAGS - with longer runs they seem to agree (as best as one can tell given the mixing without doing a super long run) + # looks fine, but alpha and beta values shifted a bit (systematically) relative to JAGS results - on further inspection this is because mixing for this model is poor in both NIMBLE and JAGS - with longer runs they seem to agree (as best as one can tell given the mixing without doing a super long run) }) test_mcmc('beetles', model = 'beetles-logit.bug', inits = 'beetles-inits.R', data = 'beetles-data.R', numItsC = 1000, resampleData = TRUE) - # getting warning; deterministic model node is NA or NaN in model initialization - # weirdness with llike.sat[8] being NaN on init (actually that makes sense), and with weird lifting of RHS of llike.sat +# getting warning; deterministic model node is NA or NaN in model initialization +# weirdness with llike.sat[8] being NaN on init (actually that makes sense), and with weird lifting of RHS of llike.sat test_that('leuk example setup', { writeLines(c("var","Y[N,T],","dN[N,T];"), con = file.path(tempdir(), "leuk.bug")) ## echo doesn't seem to work on Windows - # need nimStep in data block as we no longer have step + # need nimStep in data block as we no longer have step system.in.dir(paste("cat leuk.bug >> ", file.path(tempdir(), "leuk.bug")), dir = system.file('classic-bugs','vol1','leuk',package = 'nimble')) - # need nimStep in data block as we no longer have step + # need nimStep in data block as we no longer have step system.in.dir(paste("sed -i -e 's/step/nimStep/g'", file.path(tempdir(), "leuk.bug"))) test_mcmc(model = file.path(tempdir(), "leuk.bug"), name = 'leuk', inits = system.file('classic-bugs', 'vol1', 'leuk','leuk-init.R', package = 'nimble'), data = system.file('classic-bugs', 'vol1', 'leuk','leuk-data.R', package = 'nimble'), numItsC = 1000, @@ -145,14 +145,14 @@ test_that('salm example setup', { writeLines(paste("var","logx[doses];"), con = file.path(tempdir(), "salm.bug")) system.in.dir(paste("cat salm.bug >>", file.path(tempdir(), "salm.bug")), dir = system.file('classic-bugs','vol1','salm', package = 'nimble')) test_mcmc(model = file.path(tempdir(), "salm.bug"), name = 'salm', inits = system.file('classic-bugs', 'vol1', 'salm','salm-init.R', package = 'nimble'), data = system.file('classic-bugs', 'vol1', 'salm','salm-data.R', package = 'nimble'), numItsC = 1000, avoidNestedTest = TRUE) - # looks good compared to JAGS + # looks good compared to JAGS }) test_that('air example setup', { file.copy(system.file('classic-bugs','vol2','air','air.bug', package = 'nimble'), file.path(tempdir(), "air.bug"), overwrite=TRUE) system.in.dir(paste("sed -i -e 's/mean(X)/mean(X\\[\\])/g'", file.path(tempdir(), "air.bug"))) test_mcmc(model = file.path(tempdir(), "air.bug"), name = 'air', inits = system.file('classic-bugs', 'vol2', 'air','air-inits.R', package = 'nimble'), data = system.file('classic-bugs', 'vol2', 'air','air-data.R', package = 'nimble'), numItsC = 1000, avoidNestedTest = TRUE) - # theta[2] posterior is a bit off from JAGS - would be worth more investigation + # theta[2] posterior is a bit off from JAGS - would be worth more investigation }) test_that('jaw-linear setup', { @@ -160,7 +160,7 @@ test_that('jaw-linear setup', { test_mcmc(model = file.path(tempdir(), "jaw-linear.bug"), name = 'jaw-linear', inits = system.file('classic-bugs', 'vol2', 'jaw','jaw-inits.R', package = 'nimble'), data = system.file('classic-bugs', 'vol2', 'jaw','jaw-data.R', package = 'nimble'), numItsC = 1000, avoidNestedTest = TRUE) # , knownFailures = list('R MCMC' = 'Cholesky of NA matrix fails in R 3.4.2 in calculate(model) of initializeModel() but not in R 3.4.1')) }) ## note R MCMC used to fail when tried to do Cholesky of 0 matrix in 2-point method, but no longer doing multiplicative link for Wishart targets - + test_mcmc('pump', resampleData = TRUE, results = list(mean = list( @@ -229,53 +229,53 @@ test_that('very simple example setup', { ### linear Gaussian state-space model (of length 5) test_that('linear Gaussian state-space model MCMC works', { - set.seed(0) - n <- 5 - a <- 6 - b <- 0.8 - sigmaPN <- 2 - sigmaOE <- 4 - set.seed(0) - x <- numeric(n) - y <- numeric(n) - x[1] <- 1 - y[1] <- rnorm(1, x[1], sigmaOE) - for(i in 2:n) { - x[i] <- rnorm(1, a+b*x[i-1], sigmaPN) - y[i] <- rnorm(1, x[i], sigmaOE) - } - code <- nimbleCode({ - a ~ dnorm(0, sd = 10000) - b ~ dnorm(0, sd = 10000) - sigmaPN ~ dunif(0, 10000) - sigmaOE ~ dunif(0, 10000) - x[1] ~ dnorm(0, sd = 10000) - y[1] ~ dnorm(x[1], sd = sigmaOE) - for(t in 2:N) { - x[t] ~ dnorm(a + b*x[t-1], sd = sigmaPN) - y[t] ~ dnorm(x[t], sd = sigmaOE) + set.seed(0) + n <- 5 + a <- 6 + b <- 0.8 + sigmaPN <- 2 + sigmaOE <- 4 + set.seed(0) + x <- numeric(n) + y <- numeric(n) + x[1] <- 1 + y[1] <- rnorm(1, x[1], sigmaOE) + for(i in 2:n) { + x[i] <- rnorm(1, a+b*x[i-1], sigmaPN) + y[i] <- rnorm(1, x[i], sigmaOE) } - }) - constants <- list(N = length(y)) - data <- list(y = y) - inits <- list(a=6, b=0.8, sigmaOE=4, sigmaPN=2, x=y+rnorm(length(y))) - - test_mcmc(model = code, - name = 'linear Gaussian state-space model', - data = c(data, constants), - inits = inits, - seed = 123, - resampleData = FALSE, - results = list( - mean = list(a = 37.36, b = -2.26, sigmaPN = 5.27, sigmaOE = 5.08), - sd = list(a = 30.14, b = 2.60, sigmaPN = 6.54, sigmaOE = 2.96)), - ## The expected results come from exactly these run conditions, - ## so they are essentially exact MCMC chain replication results. - ## For that reason, tolerances are all set to the precision with - ## which the numbers were entered, to nearest 0.01. - resultsTolerance = list(mean = list(a = .01, b = .01, sigmaPN=0.01, sigmaOE=0.01), - sd = list(a = .01, b = .01, sigmaPN=0.01, sigmaOE=0.01)), - avoidNestedTest = TRUE) + code <- nimbleCode({ + a ~ dnorm(0, sd = 10000) + b ~ dnorm(0, sd = 10000) + sigmaPN ~ dunif(0, 10000) + sigmaOE ~ dunif(0, 10000) + x[1] ~ dnorm(0, sd = 10000) + y[1] ~ dnorm(x[1], sd = sigmaOE) + for(t in 2:N) { + x[t] ~ dnorm(a + b*x[t-1], sd = sigmaPN) + y[t] ~ dnorm(x[t], sd = sigmaOE) + } + }) + constants <- list(N = length(y)) + data <- list(y = y) + inits <- list(a=6, b=0.8, sigmaOE=4, sigmaPN=2, x=y+rnorm(length(y))) + + test_mcmc(model = code, + name = 'linear Gaussian state-space model', + data = c(data, constants), + inits = inits, + seed = 123, + resampleData = FALSE, + results = list( + mean = list(a = 37.36, b = -2.26, sigmaPN = 5.27, sigmaOE = 5.08), + sd = list(a = 30.14, b = 2.60, sigmaPN = 6.54, sigmaOE = 2.96)), + ## The expected results come from exactly these run conditions, + ## so they are essentially exact MCMC chain replication results. + ## For that reason, tolerances are all set to the precision with + ## which the numbers were entered, to nearest 0.01. + resultsTolerance = list(mean = list(a = .01, b = .01, sigmaPN=0.01, sigmaOE=0.01), + sd = list(a = .01, b = .01, sigmaPN=0.01, sigmaOE=0.01)), + avoidNestedTest = TRUE) }) @@ -368,7 +368,7 @@ test_that('slice sampler example setup', { avoidNestedTest = TRUE) nimbleOptions(MCMCusePredictiveDependenciesInCalculations = nimbleUsePredictiveDependenciesSetting) - }) +}) ### elliptical slice sampler 'ess' @@ -403,7 +403,7 @@ test_that('elliptical slice sampler setup', { numItsC = 100000, samplers = list(list(type = 'ess', target = 'x')), avoidNestedTest = TRUE) - }) +}) ### demo2 of check conjugacy @@ -463,7 +463,7 @@ test_that('Weibull-gamma conjugacy setup', { conf <- configureMCMC(m) samplers <- conf$getSamplers() expect_identical(samplers[[1]]$name, 'conjugate_dgamma_dweib_multiplicative', - info = "dweibull-dgamma conjugacy with dependency using lambda not detected") + info = "dweibull-dgamma conjugacy with dependency using lambda not detected") mcmc <- buildMCMC(conf) comp <- compileNimble(m, mcmc) set.seed(0) @@ -481,18 +481,18 @@ test_that('Weibull-gamma conjugacy setup', { smpMan <- manualSampler(10, y, depShape, c, shape, rate) expect_identical(smp[,1], smpMan, - info = "NIMBLE gamma-Weibull conjugate sampler and manual sampler results differ") + info = "NIMBLE gamma-Weibull conjugate sampler and manual sampler results differ") }) test_that('Dirichlet-multinomial conjugacy setup', { -### Dirichlet-multinomial conjugacy + ### Dirichlet-multinomial conjugacy -# as of v0.4, exact numerical results here have changed because -# ddirch now sometimes returns NaN rather than -Inf (when an -# alpha is proposed to be negative) -- this changes the RNG -# sequence because NaN values result in no runif() call in decide() + # as of v0.4, exact numerical results here have changed because + # ddirch now sometimes returns NaN rather than -Inf (when an + # alpha is proposed to be negative) -- this changes the RNG + # sequence because NaN values result in no runif() call in decide() -# single multinomial + # single multinomial set.seed(0) n <- 100 alpha <- c(10, 30, 15, 60, 1) @@ -517,7 +517,7 @@ test_that('Dirichlet-multinomial conjugacy setup', { resultsTolerance = list(mean = list(p = rep(.06, K))), avoidNestedTest = TRUE) }) ## bad mixing for alphas; probably explains why posterior estimates for alphas changed so much as of v 0.4 - + ## with replication test_that('Dirichlet-multinomial with replication setup', { set.seed(0) @@ -552,7 +552,7 @@ test_that('Dirichlet-multinomial with replication setup', { resultsTolerance = list(mean = list(p = matrix(.05, m, K), alpha = c(5,10,10,20,.5))), knownFailures = list('MCMC match to known posterior: p mean 39' = 'KNOWN ISSUE: two samples outside resultsTolerance', - 'MCMC match to known posterior: p mean 76' = 'KNOWN ISSUE: two samples outside resultsTolerance'), avoidNestedTest = TRUE) + 'MCMC match to known posterior: p mean 76' = 'KNOWN ISSUE: two samples outside resultsTolerance'), avoidNestedTest = TRUE) }) # note alphas mix poorly (and are highly correlated), # presumably because of cross-level dependence between @@ -560,9 +560,9 @@ test_that('Dirichlet-multinomial with replication setup', { # or, of course, integrating over the p's test_that('Dirichlet-categorical conjugacy setup', { -### Dirichlet-categorical conjugacy + ### Dirichlet-categorical conjugacy -# single multinomial represented as categorical + # single multinomial represented as categorical set.seed(0) n <- 100 alpha <- c(10, 30, 15, 60, 1) @@ -573,7 +573,7 @@ test_that('Dirichlet-categorical conjugacy setup', { code <- function() { for(i in 1:n) - y[i] ~ dcat(p[1:K]) + y[i] ~ dcat(p[1:K]) p[1:K] ~ ddirch(alpha[1:K]) for(i in 1:K) { alpha[i] ~ dgamma(.001, .001); @@ -613,7 +613,7 @@ test_that('block sampler on MVN node setup', { var = list(x = c(.1, .03, .01))), samplers = list( list(type = 'RW_block', target = 'x[1:3]')), avoidNestedTest = TRUE) - # caution: setting targetNodes='x' works but the initial end sampler is not removed because x[1:3] in targetNode in default sampler != 'x' in targetNodes passed in + # caution: setting targetNodes='x' works but the initial end sampler is not removed because x[1:3] in targetNode in default sampler != 'x' in targetNodes passed in if(FALSE) { Rmodel <- nimbleModel(code, constants = list(Q=Q)) mcmcspec <- MCMCspec(Rmodel, nodes = NULL) @@ -643,7 +643,7 @@ test_that('block sampler on MVN node setup', { }) test_that('second block sampler on multivariate node', { -### DT's model + ### DT's model mu <- c(1,2,3) corr <- matrix(c(1,.8,0.3,.8,1,0,0.3,0,1), nrow=3) varr <- c(1,2,3) @@ -663,10 +663,10 @@ test_that('second block sampler on multivariate node', { var = list(x = c(.1,.1,.1))), samplers = list( list(type = 'RW_block', target = 'x[1:3]')), avoidNestedTest = TRUE) - }) +}) test_that('MVN conjugate setup', { -### MVN conjugate update + ### MVN conjugate update set.seed(0) mu0 = 1:3 @@ -675,7 +675,7 @@ test_that('MVN conjugate setup', { a = c(-2, .5, 1) B = matrix(rnorm(9), 3) -##### not currently working - see Perry's email of ~ 10/6/14 + ##### not currently working - see Perry's email of ~ 10/6/14 ## code <- nimbleCode({ ## mu[1:3] ~ dmnorm(mu0[1:3], Q0[1:3, 1:3]) ## y[1:3] ~ dmnorm(asCol(a[1:3]) + B[1:3, 1:3] %*% asCol(mu[1:3]), Q[1:3, 1:3]) @@ -689,7 +689,7 @@ test_that('MVN conjugate setup', { mu <- mu0 + chol(solve(Q0)) %*% rnorm(3) - # make sure y is a vec not a 1-col matrix or get a dimensionality error + # make sure y is a vec not a 1-col matrix or get a dimensionality error y <- c(a + B%*%mu + chol(solve(Q)) %*% rnorm(3)) data = list(mu0 = mu0, Q0 = Q0, Q = Q, a = a, B = B, y = y) @@ -703,7 +703,7 @@ test_that('MVN conjugate setup', { cov = list(mu = matrix(.01, 3, 3))), avoidNestedTest = TRUE) -### scalar RW updates in place of conjugate mv update + ### scalar RW updates in place of conjugate mv update test_mcmc(model = code, name = 'two-level multivariate normal with scalar updaters', data = data, seed = 0, numItsC = 100000, results = list(mean = list(mu = muMeanTrue), @@ -839,7 +839,7 @@ test_that('conjugate Wishart setup', { sd = list(Omega = OmegaSimTrueSDs)), resultsTolerance = list(mean = list(Omega = matrix(.05, M,M)), sd = list(Omega = matrix(0.06, M, M))), avoidNestedTest = TRUE) - # issue with Chol in R MCMC - probably same issue as in jaw-linear + # issue with Chol in R MCMC - probably same issue as in jaw-linear }) @@ -889,7 +889,7 @@ test_that('conjugate Wishart setup with scaling', { sd = list(Omega = OmegaSimTrueSDs)), resultsTolerance = list(mean = list(Omega = matrix(.05, M,M)), sd = list(Omega = matrix(0.06, M, M))), avoidNestedTest = TRUE) - # issue with Chol in R MCMC - probably same issue as in jaw-linear + # issue with Chol in R MCMC - probably same issue as in jaw-linear }) @@ -1146,7 +1146,7 @@ test_that('using LKJ randomw walk samplers', { outSigma <- matrix(0, nrow(out), p*p) for(i in 1:nrow(outSigma)) outSigma[i,] <- t(matrix(out[i,], p, p)) %*% matrix(out[i,],p,p) - + conf <- configureMCMC(m, nodes = NULL, thin = 10) conf$addSampler('Ustar', 'RW_lkj_corr_cholesky', control = list(scale = .1)) mcmc <- buildMCMC(conf) @@ -1374,7 +1374,7 @@ test_that('conjugate MVN with ragged dependencies', { expect_true(all(abs(pmean - obsmean) / pmean < 0.01), info = 'ragged dmnorm conjugate posterior mean') expect_true(all(abs(pprec - obsprec) / pprec < 0.005), info = 'ragged dmnorm conjugate posterior precision') }) - + ## testing binary sampler test_that('binary sampler setup', { cat('===== Starting MCMC test for binary sampler. =====') @@ -1441,8 +1441,8 @@ test_that('binary sampler setup', { expect_lt(abs(means[['h']] - 0.5), tol) }) - - ## testing the binary sampler handles 'out of bounds' ok + +## testing the binary sampler handles 'out of bounds' ok test_that('binary sampler handles out of bounds', { cat('===== Starting MCMC test for binary sampler handles out of bounds. =====') code <- nimbleCode({ @@ -1596,9 +1596,9 @@ test_that('RW_dirichlet sampler consistent with conjugate multinomial sampler', means <- apply(Csamples, 2, mean) expect_lt(max(abs(means[c('p[1]','p[2]','p[3]')] - means[c('p2[1]','p2[2]','p2[3]')])), 0.001, - label = 'agreement between RW_dirichlet and conjugate dirichlet sampling' ) + label = 'agreement between RW_dirichlet and conjugate dirichlet sampling' ) }) - + ## testing RW_dirichlet sampler ## more complicated -- intermediate deterministic nodes, and non-conjugate ## test agreement between RW_dirichlet sampler, and writing model with component gammas @@ -1663,8 +1663,8 @@ test_that('dnorm-dmnorm conjugacies NIMBLE fails to detect', { m = nimbleModel(code, inits = list(pr0 = diag(3), pr = diag(3))) conf <- configureMCMC(m) expect_failure(expect_match(conf$getSamplers()[[1]]$name, "conjugate_dmnorm_dnorm", - info = "failed to detect dmnorm-dnorm conjugacy"), - info = "EXPECTED FAILURE NOT FAILING: this known failure should occur because of limitations in conjugacy detection with dmnorm dependents of dnorm target") + info = "failed to detect dmnorm-dnorm conjugacy"), + info = "EXPECTED FAILURE NOT FAILING: this known failure should occur because of limitations in conjugacy detection with dmnorm dependents of dnorm target") code = nimbleCode({ for(i in 1:3) @@ -1675,8 +1675,8 @@ test_that('dnorm-dmnorm conjugacies NIMBLE fails to detect', { m = nimbleModel(code, inits = list(z = rep(0,3), pr = diag(3))) conf <- configureMCMC(m) expect_failure(expect_match(conf$getSamplers()[[2]]$name, "conjugate_dnorm_dmnorm", - info = "failed to detect dmnorm-dnorm conjugacy"), - info = "EXPECTED FAILURE NOT FAILING: this known failure should occur because of limitations in conjugacy detection with dnorm dependents of dmnorm target") + info = "failed to detect dmnorm-dnorm conjugacy"), + info = "EXPECTED FAILURE NOT FAILING: this known failure should occur because of limitations in conjugacy detection with dnorm dependents of dmnorm target") }) ## dnorm prior in vectorized regression mean (inprod, matrix multiplication) @@ -1697,7 +1697,7 @@ test_that('NIMBLE detects dnorm-dnorm conjugacy via inprod() or %*%', { m <- nimbleModel(code, data = data, constants = constants) conf <- configureMCMC(m) expect_identical(conf$getSamplers()[[1]]$name, 'conjugate_dnorm_dnorm_linear', - info = "conjugacy with inprod not detected") + info = "conjugacy with inprod not detected") code <- nimbleCode({ for(i in 1:n) @@ -1712,7 +1712,7 @@ test_that('NIMBLE detects dnorm-dnorm conjugacy via inprod() or %*%', { m <- nimbleModel(code, data = data, constants = constants) conf <- configureMCMC(m) expect_identical(conf$getSamplers()[[1]]$name, 'conjugate_dnorm_dnorm_linear', - info = "conjugacy with inprod not detected") + info = "conjugacy with inprod not detected") ## compare to conjugate sampler using summed contributions @@ -1731,7 +1731,7 @@ test_that('NIMBLE detects dnorm-dnorm conjugacy via inprod() or %*%', { m <- nimbleModel(code, data = data, constants = constants) conf <- configureMCMC(m) expect_identical(conf$getSamplers()[[1]]$name, 'conjugate_dnorm_dnorm_linear', - info = "conjugacy with sum not detected") + info = "conjugacy with sum not detected") ## compare to conjugate sampler using summed contributions @@ -1766,7 +1766,7 @@ test_that('NIMBLE detects dnorm-dnorm conjugacy via inprod() or %*%', { m <- nimbleModel(code, data = data, constants = constants) conf <- configureMCMC(m) expect_identical(conf$getSamplers()[[1]]$name, 'RW', - info = "conjugacy with inprod improperly detected") + info = "conjugacy with inprod improperly detected") code <- nimbleCode({ for(i in 1:n) @@ -1778,7 +1778,7 @@ test_that('NIMBLE detects dnorm-dnorm conjugacy via inprod() or %*%', { m <- nimbleModel(code, data = data, constants = constants) conf <- configureMCMC(m) expect_identical(conf$getSamplers()[[1]]$name, 'RW', - info = "conjugacy with inprod improperly detected") + info = "conjugacy with inprod improperly detected") code <- nimbleCode({ for(i in 1:n) @@ -1790,7 +1790,7 @@ test_that('NIMBLE detects dnorm-dnorm conjugacy via inprod() or %*%', { m <- nimbleModel(code, data = data, constants = constants) conf <- configureMCMC(m) expect_identical(conf$getSamplers()[[1]]$name, 'RW', - info = "conjugacy with sum improperly detected") + info = "conjugacy with sum improperly detected") code <- nimbleCode({ for(i in 1:n) @@ -1802,7 +1802,7 @@ test_that('NIMBLE detects dnorm-dnorm conjugacy via inprod() or %*%', { m <- nimbleModel(code, data = data, constants = constants) conf <- configureMCMC(m) expect_identical(conf$getSamplers()[[1]]$name, 'conjugate_dnorm_dnorm_linear', - info = "conjugacy with matrix multiplication not detected") + info = "conjugacy with matrix multiplication not detected") mcmc <- buildMCMC(conf) cm <- compileNimble(m) cmcmc <- compileNimble(mcmc, project = m) @@ -1835,8 +1835,8 @@ test_that('NIMBLE detects dnorm-dnorm conjugacy via inprod() or %*%', { m <- nimbleModel(code, data = data, constants = constants) conf <- configureMCMC(m) expect_identical(conf$getSamplers()[[1]]$name, 'conjugate_dnorm_dnorm_linear', - info = "conjugacy with inprod not detected") - + info = "conjugacy with inprod not detected") + code <- nimbleCode({ for(i in 1:n) y[i] ~ dnorm(b0 + inprod(zbeta[1:p], X[i, 1:p]), 1) @@ -1854,7 +1854,7 @@ test_that('NIMBLE detects dnorm-dnorm conjugacy via inprod() or %*%', { m <- nimbleModel(code, data = data, constants = constants) conf <- configureMCMC(m) expect_identical(conf$getSamplers()[[1]]$name, 'conjugate_dnorm_dnorm_linear', - info = "conjugacy with inprod not detected") + info = "conjugacy with inprod not detected") code <- nimbleCode({ for(i in 1:n) @@ -1873,7 +1873,7 @@ test_that('NIMBLE detects dnorm-dnorm conjugacy via inprod() or %*%', { expect_identical(conf$getSamplers()[[1]]$name, 'conjugate_dnorm_dnorm_linear', info = "conjugacy with inprod not detected") - code <- nimbleCode({ + code <- nimbleCode({ for(i in 1:n) y[i] ~ dnorm(b0 + inprod(zbeta[1:p], X[i, 1:p]), 1) for(i in 1:p) { @@ -1893,11 +1893,11 @@ test_that('NIMBLE detects dnorm-dnorm conjugacy via inprod() or %*%', { info = "conjugacy with inprod mistakenly detected") expect_identical(nimble:::cc_checkLinearity( - quote(structureExpr(z[1] * exp(w * beta[1]), z[2] * exp(w * beta[2]))), - 'beta[2]'), NULL) + quote(structureExpr(z[1] * exp(w * beta[1]), z[2] * exp(w * beta[2]))), + 'beta[2]'), NULL) output <- nimble:::cc_checkLinearity( - quote(structureExpr(z[1] * exp(w * beta[1]), a + z[2] * (d + w * beta[2]))), - 'beta[2]') + quote(structureExpr(z[1] * exp(w * beta[1]), a + z[2] * (d + w * beta[2]))), + 'beta[2]') expect_identical(is.list(output), TRUE) ## should be a list with scale/offset }) @@ -1987,15 +1987,15 @@ test_that('cc_checkScalar operates correctly', { test_that('cc_stripExpr operates correctly', { expr <- 'coeff * (log(value) - offset) * taulog' expect_identical(deparse(nimble:::cc_stripExpr(parse(text = 'coeff^2 * tau')[[1]], TRUE, TRUE)), - '1 * tau') + '1 * tau') expect_identical(deparse(nimble:::cc_stripExpr(parse(text = expr)[[1]], TRUE, TRUE)), - '(log(value)) * taulog') + '(log(value)) * taulog') expect_identical(deparse(nimble:::cc_stripExpr(parse(text = expr)[[1]], TRUE, FALSE)), - 'coeff * (log(value)) * taulog') + 'coeff * (log(value)) * taulog') expect_identical(deparse(nimble:::cc_stripExpr(parse(text = expr)[[1]], FALSE, TRUE)), - '(log(value) - offset) * taulog') + '(log(value) - offset) * taulog') expect_identical(deparse(nimble:::cc_stripExpr(parse(text = expr)[[1]], FALSE, FALSE)), - expr) + expr) }) @@ -2016,7 +2016,7 @@ test_that("realized conjugacy links are working", { mu[2] ~ dnorm(0,1) }) m <- nimbleModel(code, data = list (y1 = rnorm(2), - y2=rnorm(2), y3=rnorm(2), y4= rnorm(2), y5 = rnorm(2)), + y2=rnorm(2), y3=rnorm(2), y4= rnorm(2), y5 = rnorm(2)), inits = list(mu = rnorm(2), x = matrix(rnorm(4), 2))) conf <- configureMCMC(m) mcmc <- buildMCMC(conf) @@ -2028,13 +2028,13 @@ test_that("realized conjugacy links are working", { expect_identical(mcmc$samplerFunctions[[1]]$N_dep_dnorm_linear, 4L) expect_identical(c('dep_dnorm_identity_coeff', 'dep_dnorm_additive_coeff') %in% - ls(mcmc$samplerFunctions[[1]]), rep(FALSE, 2)) + ls(mcmc$samplerFunctions[[1]]), rep(FALSE, 2)) expect_identical(c('dep_dnorm_multiplicative_coeff', 'dep_dnorm_linear_coeff') %in% - ls(mcmc$samplerFunctions[[1]]), rep(TRUE, 2)) + ls(mcmc$samplerFunctions[[1]]), rep(TRUE, 2)) expect_identical(c('dep_dnorm_identity_offset', 'dep_dnorm_multiplicative_offset') %in% - ls(mcmc$samplerFunctions[[1]]), rep(FALSE, 2)) + ls(mcmc$samplerFunctions[[1]]), rep(FALSE, 2)) expect_identical(c('dep_dnorm_additive_offset', 'dep_dnorm_linear_offset') %in% - ls(mcmc$samplerFunctions[[1]]), rep(TRUE, 2)) + ls(mcmc$samplerFunctions[[1]]), rep(TRUE, 2)) expect_identical(mcmc$samplerFunctions[[1]]$dep_dnorm_identity_nodeNames, c('y1[1]', 'y1[2]')) expect_identical(mcmc$samplerFunctions[[1]]$dep_dnorm_additive_nodeNames, c('y2[1]', 'y2[2]')) @@ -2054,7 +2054,7 @@ test_that("realized conjugacy links are working", { tau ~ dgamma(1, 1) }) m <- nimbleModel(code, data = list (y1 = rnorm(2), - y2=rnorm(2), y3= rpois(2, 1), y4 = rpois(2, 1)), + y2=rnorm(2), y3= rpois(2, 1), y4 = rpois(2, 1)), inits = list(mu = rnorm(1), tau = 1)) conf <- configureMCMC(m) mcmc <- buildMCMC(conf) @@ -2066,13 +2066,13 @@ test_that("realized conjugacy links are working", { expect_identical(mcmc$samplerFunctions[[1]]$N_dep_dpois_multiplicative, 2L) expect_identical(c('dep_dnorm_identity_offset', 'dep_dnorm_multiplicative_offset', 'dep_dnorm_additive_offset', 'dep_dnorm_linear_offset') %in% - ls(mcmc$samplerFunctions[[1]]), rep(FALSE, 4)) + ls(mcmc$samplerFunctions[[1]]), rep(FALSE, 4)) expect_identical(c('dep_dnorm_identity_coeff', 'dep_dnorm_multiplicative_coeff', 'dep_dnorm_additive_coeff', 'dep_dnorm_linear_coeff') %in% - ls(mcmc$samplerFunctions[[1]]), c(FALSE, TRUE, FALSE, FALSE)) + ls(mcmc$samplerFunctions[[1]]), c(FALSE, TRUE, FALSE, FALSE)) expect_identical(c('dep_dpois_identity_offset', 'dep_dpois_multiplicative_offset', 'dep_dpois_additive_offset', 'dep_dpois_linear_offset') %in% - ls(mcmc$samplerFunctions[[1]]), rep(FALSE, 4)) + ls(mcmc$samplerFunctions[[1]]), rep(FALSE, 4)) expect_identical(c('dep_dpois_identity_coeff', 'dep_dpois_multiplicative_coeff', 'dep_dpois_additive_coeff', 'dep_dpois_linear_coeff') %in% - ls(mcmc$samplerFunctions[[1]]), c(FALSE, TRUE, FALSE, FALSE)) + ls(mcmc$samplerFunctions[[1]]), c(FALSE, TRUE, FALSE, FALSE)) expect_identical(mcmc$samplerFunctions[[1]]$dep_dnorm_identity_nodeNames, c('y1[1]', 'y1[2]')) expect_identical(mcmc$samplerFunctions[[1]]$dep_dnorm_multiplicative_nodeNames, c('y2[1]', 'y2[2]')) @@ -2098,9 +2098,9 @@ test_that("realized conjugacy links are working", { } }) m <- nimbleModel(code, data = list (y1 = matrix(rnorm(6),2), - y2 = matrix(rnorm(6),2), - y3 = matrix(rnorm(6),2), - y4 = matrix(rnorm(6),2)), + y2 = matrix(rnorm(6),2), + y3 = matrix(rnorm(6),2), + y4 = matrix(rnorm(6),2)), inits = list(b0 = rnorm(3), A=matrix(1:9, 3), pr = diag(3))) conf <- configureMCMC(m) mcmc <- buildMCMC(conf) @@ -2112,13 +2112,13 @@ test_that("realized conjugacy links are working", { expect_identical(mcmc$samplerFunctions[[1]]$N_dep_dmnorm_linear, 2L) expect_identical(c('dep_dmnorm_identity_coeff', 'dep_dmnorm_additive_coeff') %in% - ls(mcmc$samplerFunctions[[1]]), rep(FALSE, 2)) + ls(mcmc$samplerFunctions[[1]]), rep(FALSE, 2)) expect_identical(c('dep_dmnorm_multiplicative_coeff', 'dep_dmnorm_linear_coeff') %in% - ls(mcmc$samplerFunctions[[1]]), rep(TRUE, 2)) + ls(mcmc$samplerFunctions[[1]]), rep(TRUE, 2)) expect_identical(c('dep_dmnorm_identity_offset', 'dep_dmnorm_multiplicative_offset') %in% - ls(mcmc$samplerFunctions[[1]]), rep(FALSE, 2)) + ls(mcmc$samplerFunctions[[1]]), rep(FALSE, 2)) expect_identical(c('dep_dmnorm_additive_offset', 'dep_dmnorm_linear_offset') %in% - ls(mcmc$samplerFunctions[[1]]), rep(TRUE, 2)) + ls(mcmc$samplerFunctions[[1]]), rep(TRUE, 2)) expect_identical(mcmc$samplerFunctions[[1]]$dep_dmnorm_identity_nodeNames, c('y1[1, 1:3]', 'y1[2, 1:3]')) expect_identical(mcmc$samplerFunctions[[1]]$dep_dmnorm_additive_nodeNames, c('y2[1, 1:3]', 'y2[2, 1:3]')) @@ -2137,7 +2137,7 @@ test_that("realized conjugacy links are working", { pr[1:3,1:3] ~ dwish(R[1:3,1:3], 8) }) m <- nimbleModel(code, data = list (y1 = matrix(rnorm(6),2), - y2 = matrix(rnorm(6),2)), + y2 = matrix(rnorm(6),2)), inits = list(pr = diag(3), R = diag(3))) conf <- configureMCMC(m) mcmc <- buildMCMC(conf) @@ -2147,11 +2147,11 @@ test_that("realized conjugacy links are working", { expect_identical(mcmc$samplerFunctions[[1]]$N_dep_dmnorm_multiplicativeScalar, 2L) expect_identical('dep_dmnorm_identity_coeff' %in% - ls(mcmc$samplerFunctions[[1]]), FALSE) + ls(mcmc$samplerFunctions[[1]]), FALSE) expect_identical('dep_dmnorm_multiplicativeScalar_coeff' %in% - ls(mcmc$samplerFunctions[[1]]), TRUE) + ls(mcmc$samplerFunctions[[1]]), TRUE) expect_identical(c('dep_dmnorm_identity_offset', 'dep_dmnorm_multiplicativeScalar_offset') %in% - ls(mcmc$samplerFunctions[[1]]), rep(FALSE, 2)) + ls(mcmc$samplerFunctions[[1]]), rep(FALSE, 2)) expect_identical(mcmc$samplerFunctions[[1]]$dep_dmnorm_identity_nodeNames, c('y1[1, 1:3]', 'y1[2, 1:3]')) expect_identical(mcmc$samplerFunctions[[1]]$dep_dmnorm_multiplicativeScalar_nodeNames, c('y2[1, 1:3]', 'y2[2, 1:3]')) @@ -3021,6 +3021,342 @@ 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]') + 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){ + 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]') + 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', { + 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]') + 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', { + 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]) + 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', { + 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]) + 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', { + 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 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] + 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 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('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) adj1 <- c(2, 1, 3, 2, 4, 3)