diff --git a/packages/nimble/DESCRIPTION b/packages/nimble/DESCRIPTION index dd4de7382..7e8bc8b00 100644 --- a/packages/nimble/DESCRIPTION +++ b/packages/nimble/DESCRIPTION @@ -134,6 +134,7 @@ Collate: MCMC_build.R MCMC_run.R MCMC_samplers.R + MCMC_derived.R MCMC_conjugacy.R MCMC_autoBlock.R MCMC_RJ.R diff --git a/packages/nimble/NAMESPACE b/packages/nimble/NAMESPACE index 48c6b6736..943031b9d 100644 --- a/packages/nimble/NAMESPACE +++ b/packages/nimble/NAMESPACE @@ -85,6 +85,11 @@ export(dcar_proper) export(dcat) export(dconstraint) export(ddexp) +export(derived_BASE) +export(derived_logProb) +export(derived_mean) +export(derived_predictive) +export(derived_variance) export(ddirch) export(decide) export(decideAndJump) diff --git a/packages/nimble/R/MCMC_autoBlock.R b/packages/nimble/R/MCMC_autoBlock.R index cab187004..8755f59dd 100644 --- a/packages/nimble/R/MCMC_autoBlock.R +++ b/packages/nimble/R/MCMC_autoBlock.R @@ -416,7 +416,7 @@ autoBlockClass <- setRefClass( ## msg <- 'using \'posterior_predictive\' sampler may lead to results we don\'t want' ## cat(paste0('\nWARNING: ', msg, '\n\n')); warning(msg) ## } - if(grepl('^conjugate_', ss$name) && getNimbleOption('verifyConjugatePosteriors')) { + if(grepl('^conjugate_', ss$name) && getNimbleOption('MCMCverifyConjugatePosteriors')) { ##msg <- 'conjugate sampler running slow due to checking the posterior' ##cat(paste0('\nWARNING: ', msg, '\n\n')); warning(msg) warn <- TRUE diff --git a/packages/nimble/R/MCMC_build.R b/packages/nimble/R/MCMC_build.R index 6e4c673c2..78448eb55 100644 --- a/packages/nimble/R/MCMC_build.R +++ b/packages/nimble/R/MCMC_build.R @@ -30,6 +30,7 @@ #' #' \code{initializeModel}: Boolean specifying whether to run the initializeModel routine on the underlying model object, prior to beginning MCMC sampling (default = TRUE). #' +#' \code{chain}: Integer specifying the MCMC chain number. The chain number is passed to each MCMC sampler's before_chain method. The value for this argument is specified automatically from invocation via runMCMC, and need not be supplied when calling mcmc$run (default = 1). #' \code{chain}: Integer specifying the MCMC chain number. The chain number is passed to each MCMC sampler's before_chain and after_chain methods. The value for this argument is specified automatically from invocation via runMCMC, and genernally need not be supplied when calling mcmc$run (default = 1). #' #' \code{time}: Boolean specifying whether to record runtimes of the individual internal MCMC samplers. When \code{time = TRUE}, a vector of runtimes (measured in seconds) can be extracted from the MCMC using the method \code{mcmc$getTimes()} (default = FALSE). @@ -180,16 +181,47 @@ buildMCMC <- nimbleFunction( for(i in seq_along(conf$samplerConfs)) samplerFunctions[[i]] <- conf$samplerConfs[[i]]$buildSampler(model=model, mvSaved=mvSaved) } + + ## construct mvSamples and mvSamples2 + mvSamplesConf <- conf$getMvSamplesConf(1) + mvSamples2Conf <- conf$getMvSamplesConf(2) + mvSamples <- modelValues(mvSamplesConf) + mvSamples2 <- modelValues(mvSamples2Conf) + ## build derived quantity intervals + derivedFunctions <- nimbleFunctionList(derived_BASE) + numDerived <- length(conf$derivedConfs) + derivedIntervals <- numeric(max(numDerived, 2)) ## force to be a vector + for(i in seq_along(conf$derivedConfs)) { + derivedIntervals[i] <- conf$derivedConfs[[i]]$interval + } + + ## build derived quantity function list + derivedFunctions <- nimbleFunctionList(derived_BASE) + ## code below allows 'mcmc' to be a setup argument of derived quantity nimbleFunctions + on.exit({ + for(i in seq_along(conf$derivedConfs)) { + derivedFunctions[[i]] <- conf$derivedConfs[[i]]$buildDerived(model=model, mcmc=nfRefClassObject) + } + ## need to catch the case where buildMCMC errors out early, + ## prior to executing the final lines which actually create the nfRefClassObject object: + if(exists('nfRefClassObject', inherits = FALSE)) { + nfRefClassObject[['derivedFunctions']] <- nf_preProcessMemberDataObject(get('derivedFunctions')) + } + }, add = TRUE) + + ## for naming the derivedList return object from runMCMC + derivedTypes <- sapply(conf$derivedConfs, `[[`, 'name') + if(length(derivedTypes) == 0) derivedTypes <- character() + ## used for extracting names of derived quantities, + ## having as member data is necessary for compilation + derivedNames <- character(2) + samplerExecutionOrderFromConfPlusTwoZeros <- c(conf$samplerExecutionOrder, 0, 0) ## establish as a vector monitors <- mcmc_processMonitorNames(model, conf$monitors) monitors2 <- mcmc_processMonitorNames(model, conf$monitors2) thinFromConfVec <- c(conf$thin, conf$thin2) ## vector thinToUseVec <- c(0, 0) ## vector, needs to member data - mvSamplesConf <- conf$getMvSamplesConf(1) - mvSamples2Conf <- conf$getMvSamplesConf(2) - mvSamples <- modelValues(mvSamplesConf) - mvSamples2 <- modelValues(mvSamples2Conf) samplerTimes <- c(0,0) ## establish as a vector progressBarLength <- 52 ## multiples of 4 only progressBarDefaultSetting <- getNimbleOption('MCMCprogressBar') @@ -212,6 +244,8 @@ buildMCMC <- nimbleFunction( thinWAIC <- FALSE nburnin_extraWAIC <- 0 } + firstRun <- TRUE + setupOutputs(derivedTypes) }, run = function( @@ -231,23 +265,34 @@ buildMCMC <- nimbleFunction( if(niter < 0) stop('cannot specify niter < 0') if(nburnin < 0) stop('cannot specify nburnin < 0') if(nburnin > niter) stop('cannot specify nburnin > niter') - thinToUseVec <<- thinFromConfVec - if(thin != -1) thinToUseVec[1] <<- thin - if(thin2 != -1) thinToUseVec[2] <<- thin2 - for(iThin in 1:2) { - if(thinToUseVec[iThin] < 1) stop('cannot use thin < 1') - if(thinToUseVec[iThin] != floor(thinToUseVec[iThin])) stop('cannot use non-integer thin') - } - if(initializeModel) my_initializeModel$run() - nimCopy(from = model, to = mvSaved, row = 1, logProb = TRUE) + if(firstRun) reset <- TRUE ## compulsory reset on first run of MCMC + firstRun <<- FALSE if(reset) { - samplerTimes <<- numeric(length(samplerFunctions) + 1) ## default inititialization to zero + if(initializeModel) my_initializeModel$run() + thinToUseVec <<- thinFromConfVec + if(thin != -1) thinToUseVec[1] <<- thin + if(thin2 != -1) thinToUseVec[2] <<- thin2 + for(iThin in 1:2) { + if(thinToUseVec[iThin] < 1) stop('cannot use thin < 1') + if(thinToUseVec[iThin] != floor(thinToUseVec[iThin])) stop('cannot use non-integer thin') + } + for(i in seq_along(derivedFunctions)) { + if(derivedIntervals[i] == 0) { + derivedIntervals[i] <<- thinToUseVec[1] + derivedFunctions[[i]]$set_interval(thinToUseVec[1]) + } + } for(i in seq_along(samplerFunctions)) samplerFunctions[[i]]$reset() - for(i in seq_along(samplerFunctions)) samplerFunctions[[i]]$before_chain(niter, nburnin, chain) + for(i in seq_along(derivedFunctions)) derivedFunctions[[i]]$reset() + for(i in seq_along(samplerFunctions)) samplerFunctions[[i]]$before_chain(niter, nburnin, chain) + for(i in seq_along(derivedFunctions)) derivedFunctions[[i]]$before_chain(niter-nburnin, nburnin, thinToUseVec, chain) + samplerTimes <<- numeric(length(samplerFunctions) + 1) ## default inititialization to zero mvSamples_copyRow <- 0 mvSamples2_copyRow <- 0 } else { - if(nburnin != 0) stop('cannot specify nburnin when using reset = FALSE.') + if(nburnin != 0) stop('cannot specify nburnin when using reset = FALSE.') + if(thin != -1) stop('cannot specify thin when using reset = FALSE.') + if(thin2 != -1) stop('cannot specify thin2 when using reset = FALSE.') if(dim(samplerTimes)[1] != length(samplerFunctions) + 1) samplerTimes <<- numeric(length(samplerFunctions) + 1) ## first run: default inititialization to zero if (resetMV) { mvSamples_copyRow <- 0 @@ -257,8 +302,8 @@ buildMCMC <- nimbleFunction( mvSamples2_copyRow <- getsize(mvSamples2) } } - if(onlineWAIC & resetWAIC) - waicFun[[1]]$reset() + nimCopy(from = model, to = mvSaved, row = 1, logProb = TRUE) + if(onlineWAIC & resetWAIC) waicFun[[1]]$reset() resize(mvSamples, mvSamples_copyRow + floor((niter-nburnin) / thinToUseVec[1])) resize(mvSamples2, mvSamples2_copyRow + floor((niter-nburnin) / thinToUseVec[2])) ## reinstate samplerExecutionOrder as a runtime argument, once we support non-scalar default values for runtime arguments: @@ -278,6 +323,7 @@ buildMCMC <- nimbleFunction( if(niter < 1) return() for(iter in 1:niter) { checkInterrupt() + ## execute samplerFunctions if(time) { for(i in seq_along(samplerExecutionOrderToUse)) { ind <- samplerExecutionOrderToUse[i] @@ -289,26 +335,33 @@ buildMCMC <- nimbleFunction( samplerFunctions[[ind]]$run() } } - ## adding "accumulators" to MCMC - ## https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance if(iter > nburnin) { - sampleNumber <- iter - nburnin - if(sampleNumber %% thinToUseVec[1] == 0) { + ## save samples + iterPostBurnin <- iter - nburnin + if(iterPostBurnin %% thinToUseVec[1] == 0) { mvSamples_copyRow <- mvSamples_copyRow + 1 nimCopy(from = model, to = mvSamples, row = mvSamples_copyRow, nodes = monitors) } - if(sampleNumber %% thinToUseVec[2] == 0) { + if(iterPostBurnin %% thinToUseVec[2] == 0) { mvSamples2_copyRow <- mvSamples2_copyRow + 1 nimCopy(from = model, to = mvSamples2, row = mvSamples2_copyRow, nodes = monitors2) } + ## save WAIC if(enableWAIC & onlineWAIC & iter > nburnin + nburnin_extraWAIC) { if (!thinWAIC) { waicFun[[1]]$updateStats() - } else if (sampleNumber %% thinToUseVec[1] == 0){ + } else if (iterPostBurnin %% thinToUseVec[1] == 0) { waicFun[[1]]$updateStats() } } + ## execute derivedFunctions + for(i in seq_along(derivedFunctions)) { + if(iterPostBurnin %% derivedIntervals[i] == 0) { + derivedFunctions[[i]]$run( iterPostBurnin/derivedIntervals[i] ) + } + } } + ## progress bar if(progressBar & (iter == progressBarNextFloor)) { cat('-') progressBarNext <- progressBarNext + progressBarIncrement @@ -316,7 +369,9 @@ buildMCMC <- nimbleFunction( } } if(progressBar) print('|') + ## after_chain methods for(i in seq_along(samplerFunctions)) samplerFunctions[[i]]$after_chain() + for(i in seq_along(derivedFunctions)) derivedFunctions[[i]]$after_chain() returnType(void()) }, methods = list( @@ -324,7 +379,25 @@ buildMCMC <- nimbleFunction( returnType(double(1)) return(samplerTimes[1:(length(samplerTimes)-1)]) }, - ## Old-style post-sampling WAIC calculation. + getNumDerived = function() { + returnType(double()) + return(numDerived) + }, + getDerivedQuantityResults = function(ind = double()) { + if(ind > numDerived) { + print('there aren\'t that many derived functions') + return(array(0, c(0,0))) + } + returnType(double(2)) + return(derivedFunctions[[ind]]$get_results()) + }, + getDerivedQuantityNames = function(ind = double()) { + if(ind > numDerived) print('there aren\'t that many derived functions') + returnType(character(1)) + derivedNames <<- derivedFunctions[[ind]]$get_names() + return(derivedNames) + }, + ## old-style post-sampling WAIC calculation calculateWAIC = function(nburnin = integer(default = 0)) { if(!enableWAIC) { print('Error: One must set enableWAIC = TRUE in \'configureMCMC\' or \'buildMCMC\'. See \'help(configureMCMC)\' for additional information.') diff --git a/packages/nimble/R/MCMC_configuration.R b/packages/nimble/R/MCMC_configuration.R index 3f71bb949..ff3df2676 100644 --- a/packages/nimble/R/MCMC_configuration.R +++ b/packages/nimble/R/MCMC_configuration.R @@ -11,16 +11,16 @@ samplerConf <- setRefClass( ), methods = list( initialize = function(name, samplerFunction, target, control, model) { - baseClassName <<- environment(environment(samplerFunction)$contains)$className + baseClassName <<- getBaseClassName(samplerFunction) ## environment(environment(samplerFunction)$contains)$className if(is.null(baseClassName) || (baseClassName != 'sampler_BASE')) warning('MCMC sampler nimbleFunctions should inherit from (using "contains" argument) base class sampler_BASE.') setName(name) - setSamplerFunction(samplerFunction) + setFunction(samplerFunction) setTarget(target, model) setControl(control) if(name == 'crossLevel') control <<- c(control, list(dependent_nodes = model$getDependencies(target, self = FALSE, stochOnly = TRUE))) ## special case for printing dependents of crossLevel sampler (only) }, setName = function(name) name <<- name, - setSamplerFunction = function(fun) samplerFunction <<- fun, + setFunction = function(f) samplerFunction <<- f, setTarget = function(target, model) { target <<- target targetAsScalar <<- model$expandNodeNames(target, returnScalarComponents = TRUE, sort = TRUE) @@ -36,12 +36,56 @@ samplerConf <- setRefClass( mcmc_listContentsToStr(infoList, displayControlDefaults, displayNonScalars, displayConjugateDependencies) }, show = function() { - cat(toStr()) + cat(toStr()) ## later, change to: cat(toStr(), '\n') } ) ) + +derivedConf <- setRefClass( + Class = 'derivedConf', + fields = list( + name = 'ANY', + derivedFunction = 'ANY', + baseClassName = 'ANY', + interval = 'ANY', + control = 'ANY' + ), + methods = list( + initialize = function(name, derivedFunction, interval, control, model) { + baseClassName <<- getBaseClassName(derivedFunction) ## environment(environment(derivedFunction)$contains)$className + if(is.null(baseClassName) || (baseClassName != 'derived_BASE')) warning('MCMC derived quantity nimbleFunctions should inherit from (using "contains" argument) base class derived_BASE') + setName(name) + setFunction(derivedFunction) + setInterval(interval) + setControl(control) + }, + setName = function(name) name <<- name, + setFunction = function(f) derivedFunction <<- f, + setInterval = function(interval) interval <<- interval, + setControl = function(control) control <<- control, + buildDerived = function(model, mcmc) { + derivedFunction(model=model, mcmc=mcmc, interval=interval, control=control) + }, + toStr = function(displayNonScalars = FALSE) { + s <- paste0('derived quantity: ', name, ', ') + intervalString <- paste0('execution interval: ', if(interval == 0) 'thin' else interval) + s <- paste0(s, intervalString) + if(length(control)) { + controlString <- mcmc_listContentsToStr(control, displayNonScalars = displayNonScalars, removeCfunctions = FALSE, removeLengthZero = FALSE) + if(nchar(controlString) > 0) s <- paste0(s, ', ', controlString) + } + return(s) + }, + show = function() { + cat(toStr(), '\n') + } + ) +) + + + ## 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 @@ -83,6 +127,7 @@ MCMCconf <- setRefClass( enableWAIC = 'ANY', controlWAIC = 'ANY', samplerConfs = 'ANY', + derivedConfs = 'ANY', samplerExecutionOrder = 'ANY', controlDefaults = 'ANY', unsampledNodes = 'ANY', @@ -102,6 +147,8 @@ MCMCconf <- setRefClass( onlyRW = FALSE, onlySlice = FALSE, multivariateNodesAsScalars = getNimbleOption('MCMCmultivariateNodesAsScalars'), + samplePredictiveNodes = getNimbleOption('MCMCassignSamplersToPosteriorPredictiveNodes'), + mean = FALSE, variance = FALSE, logProb = FALSE, enableWAIC = getNimbleOption('MCMCenableWAIC'), controlWAIC = list(), print = TRUE, ...) { ' @@ -140,6 +187,14 @@ onlySlice: A logical argument, with default value FALSE. If specified as TRUE, multivariateNodesAsScalars: A logical argument, with default value FALSE. If specified as TRUE, then non-terminal multivariate stochastic nodes will have scalar samplers assigned to each of the scalar components of the multivariate node. The default value of FALSE results in a single block sampler assigned to the entire multivariate node. Note, multivariate nodes appearing in conjugate relationships will be assigned the corresponding conjugate sampler (provided useConjugacy == TRUE), regardless of the value of this argument. +samplePredictiveNodes: A logical argument. When TRUE, samplers will be assigned by default to update posterior predictive model nodes. When FALSE, no samplers will be assigned to predictive nodes. The default value of this argument is given by the nimble option \'MCMCassignSamplersToPosteriorPredictiveNodes\', which itself has a default value of TRUE. + +mean: Character or logical argument. When a vector of node names is provided, the \'mean\' derived quantity function will be used to calculate the running mean for all node names specified. If TRUE, the running mean will be calculated for nodes specified in monitors. + +variance: Character or logical argument. When a vector of node names is provided, the \'variance\' derived quantity function will be used to calculate the running variance for all node names specified. If TRUE, the running variance will be calculated for nodes specified in monitors. + +logProb: When TRUE, the summed log-density of all stochastic model nodes (including data nodes) will be calculated and returned, using the logProb derived quantity function. When provided as a character vector, the individual log density of each node in this vector will be recorded. When provided as a list, each list element may contain one or mode node names, and separately for the node(s) in each element of the list, the summed log-density list will be calculated. In addition, the keyword \'.all\' may also be provided in either the vector or list argument, which corresponds to the set of all stochastic model nodes (including data). + enableWAIC: A logical argument, specifying whether to enable WAIC calculations for the resulting MCMC algorithm. Defaults to the value of nimbleOptions(\'MCMCenableWAIC\'), which in turn defaults to FALSE. Setting nimbleOptions(\'MCMCenableWAIC\' = TRUE) will ensure that WAIC is enabled for all calls to \`configureMCMC\` and \`buildMCMC\`. controlWAIC A named list of inputs that control the behavior of the WAIC calculation, passed as the \'control\' input to \'buildWAIC\'. See \'help(waic)\`. @@ -162,6 +217,7 @@ print: A logical argument specifying whether to print the montiors and samplers. enableWAIC <<- enableWAIC controlWAIC <<- controlWAIC samplerConfs <<- list() + derivedConfs <<- list() samplerExecutionOrder <<- numeric() controlDefaults <<- list(...) ##namedSamplerLabelMaker <<- labelFunctionCreator('namedSampler') ## usage long since deprecated (Dec 2020) @@ -188,7 +244,7 @@ print: A logical argument specifying whether to print the montiors and samplers. } if(missing(nodes)) { - nodes <- model$getNodeNames(stochOnly = TRUE, includeData = FALSE) + nodes <- model$getNodeNames(stochOnly = TRUE, includeData = FALSE, includePredictive = samplePredictiveNodes) # Check of all(model$isStoch(nodes)) is not needed in this case } else if(is.null(nodes) || length(nodes)==0) { nodes <- character(0) @@ -201,7 +257,14 @@ print: A logical argument specifying whether to print the montiors and samplers. onlySlice = onlySlice, multivariateNodesAsScalars = multivariateNodesAsScalars, print = FALSE) - + + if(isTRUE (mean )) addDerivedQuantity('mean' , control = list(nodes = monitors)) + if(is.character(mean )) addDerivedQuantity('mean' , control = list(nodes = mean )) + if(isTRUE (variance)) addDerivedQuantity('variance', control = list(nodes = monitors)) + if(is.character(variance)) addDerivedQuantity('variance', control = list(nodes = variance)) + if(isTRUE (logProb )) addDerivedQuantity('logProb' , control = list(nodes = '.all' )) + if(!is.logical (logProb )) addDerivedQuantity('logProb' , control = list(nodes = logProb )) + if(print) show() ##printSamplers() }, @@ -243,7 +306,7 @@ For internal use. Adds default MCMC samplers to the specified nodes. ## convert to node IDs: nodeIDs <- model$expandNodeNames(nodes, returnType = 'ids') nodeIDsOrig <- nodeIDs - + ## determine which posterior predictive nodes should be sampled with posterior_predictive sampler. ## this requires some care, because it's only those nodes for which all ## downstream dependents are also slated for sampling. @@ -727,10 +790,10 @@ For internal use only if(all(model$isData(targetOne))) return() if(any(model$isData(targetOne))) targetOne <- filterOutDataNodes(targetOne) } - newSamplerInd <- length(samplerConfs) + 1 - samplerConfs[[newSamplerInd]] <<- samplerConf(name=thisSamplerName, samplerFunction=samplerFunction, target=targetOne, control=thisControlList, model=model) - samplerExecutionOrder <<- c(samplerExecutionOrder, newSamplerInd) - if(print) printSamplers(newSamplerInd) + newInd <- length(samplerConfs) + 1 + samplerConfs[[newInd]] <<- samplerConf(name=thisSamplerName, samplerFunction=samplerFunction, target=targetOne, control=thisControlList, model=model) + samplerExecutionOrder <<- c(samplerExecutionOrder, newInd) + if(print) printSamplers(newInd) }, filterOutDataNodes = function(nodes) { @@ -859,7 +922,7 @@ Arguments: ind: A numeric vector or character vector. A numeric vector may be used to specify the indices of the samplers to print, or a character vector may be used to indicate a set of target nodes and/or variables, for which all samplers acting on these nodes will be printed. For example, printSamplers(\'x\') will print all samplers whose target is model node \'x\', or whose targets are contained (entirely or in part) in the model variable \'x\'. If omitted, then all samplers are printed. -type: a character vector containing sampler type names. Only samplers with one of these specified types, as printed by this printSamplers method, will be displayed. Standard regular expression mathing using is also applied. +type: A character vector containing sampler type names. Only samplers with one of these specified types will be displayed. Regular expression matching is also used. displayConjugateDependencies: A logical argument, specifying whether to display the dependency lists of conjugate samplers (default FALSE). @@ -878,7 +941,7 @@ byType: A logical argument, specifying whether the nodes being sampled should be if(length(ind) > 0 && max(ind) > length(samplerConfs)) stop('MCMC configuration doesn\'t have that many samplers') if(!missing(type)) { if(!is.character(type)) stop('type argument must have type character') - ## find sampler indices with 'name' matching anything in 'type' argument: + ## find indices with 'name' matching anything in 'type' argument: typeInd <- unique(unname(unlist(lapply(type, grep, x = lapply(samplerConfs, `[[`, 'name'))))) ind <- intersect(ind, typeInd) } @@ -898,7 +961,7 @@ byType: A logical argument, specifying whether the nodes being sampled should be } else conjInfo <- "non-conjugate" info <- paste0(info, ", ", conjInfo) } - cat(paste0(info, "\n")) + cat(paste0(info, '\n')) } if(!executionOrder && !identical(as.numeric(samplerExecutionOrder), as.numeric(seq_along(samplerConfs)))) { messageIfVerbose('\n [Note] Samplers have a modified order of execution.') @@ -1048,7 +1111,198 @@ The indices of execution specified in this numeric vector correspond to the enum ' return(samplerExecutionOrder) }, - + + addDerivedQuantity = function(type, + interval = 0, + control = list(), + print = FALSE, + name, + ...) { + ' +Adds a derived quantity function to the MCMCconf object. + +Arguments: + +type: Character string, specifying the type of derived quantity function to add. This character string should correspond to the name of a derived quantity nimbleFunction. Alternatively, the type argument may be provided as a nimbleFunction itself rather than its name. In that case, the \'name\' argument may also be supplied to provide a meaningful name for this function. This argument has no default value, and must be provided. + +interval: A numeric value, specifying the number of MCMC iterations between each time this derived quantity function will execute. For example, if \'interval\' is 1, then this derived quantity function will execute at the end of every MCMC sampling iteration, and if \'interval\' is 10, then this derived quantity function will execute at the end of MCMC sampling iterations 10, 20, 30, etc. If \'interval\' is omitted, then the default behavior (the default frequency of execution) will match the thinning interval (\'thin\') on which samples are saved. + +control: An optional list of control arguments to derived quantity function. + +print: Logical argument, specifying whether to print the details of newly added function. + +name: Optional character string name for the derived quantity function, which is used by the printDerivedQuantities method. If \'name\' is not provided, the \'type\' argument is used to generate the name. + +...: Additional named arguments passed through ... will be used as additional control list elements. + +Details: + +Derived quantity functions are added to the end of the list for this MCMCconf object, and do not replace any existing functions. Derived quantity functions can be removed using the removeDerivedQuantities method. + +Invisibly returns a list of the current derived quantity function configurations, which are derivedConf reference class objects. +' + nameProvided <- !missing(name) + if(is.character(type)) { + thisDerivedName <- if(nameProvided) name else gsub('^derived_', '', type) ## removes 'derived_' from beginning of name, if present + if(exists(type, inherits = TRUE) && is.nfGenerator(eval(as.name(type)))) { ## try to find derived quantity function 'type' + derivedFunction <- eval(as.name(type)) + } else { + derived_type <- paste0('derived_', type) ## next, try to find derived quantity function 'derived_type' + if(exists(derived_type) && is.nfGenerator(eval(as.name(derived_type)))) { ## try to find derived quantity function 'derived_type' + derivedFunction <- eval(as.name(derived_type)) + } else stop(paste0('cannot find derived quantity function \'', type, '\'')) + } + } else if(is.function(type)) { + if(nameProvided) { + thisDerivedName <- name + } else { + typeArg <- substitute(type) + if(is.name(typeArg) || is.call(typeArg)) { + thisDerivedName <- deparse(typeArg) + if(grepl("::", thisDerivedName)) + thisDerivedName <- gsub('.*:', '', thisDerivedName) + thisDerivedName <- gsub('^derived_', '', thisDerivedName) + } else { + thisDerivedName <- 'custom_function' + } + } + derivedFunction <- type + } else stop('derived quantity function type must be character name or a function') + if(!is.character(thisDerivedName)) stop('derived quantity function name should be a character string') + if(!is.function(derivedFunction)) stop('derived quantity function type does not specify a function') + + if(!is.numeric(interval)) stop('derived quantity interval should be numeric type') + if(length(interval) != 1) stop('derived quantity interval should be a single number') + if(interval < 0) stop('derived quantity interval must be at least 1') ## interval = 0 corresponds to matching thin interval + if(floor(interval) != interval) stop('derived quantity interval must be an integer') + + thisControlList <- c(control, list(...)) + + addOneDerivedQuantity(thisDerivedName, derivedFunction, interval, thisControlList, print) + + return(invisible(derivedConfs)) + }, + + addOneDerivedQuantity = function(thisDerivedName, derivedFunction, interval, thisControlList, print) { + ' +For internal use only +' + newInd <- length(derivedConfs) + 1 + derivedConfs[[newInd]] <<- derivedConf(name=thisDerivedName, derivedFunction=derivedFunction, interval=interval, control=thisControlList, model=model) + if(print) printDerivedQuantities(newInd) + }, + + removeDerivedQuantities = function(..., ind, print = FALSE) { + ' +Removes one or more derived quantity functions from an MCMCconf object. + +Arguments: + +...: Numeric indices, used to specify the indices of the derived quantity functions to remove. + +ind: A numeric vector specifying the indices of the derived quantity functions to remove. If omitted, and no indices are provided via the ... argument, then all derived quantity functions are removed. + +print: A logical argument specifying whether to print the current list of derived quantity functions once the removal has been done (default FALSE). + +' + if(missing(ind)) { + ind <- list(...) + ind <- unname(unlist(ind)) + if(is.null(ind)) ind <- seq_along(derivedConfs) + } + if(length(ind) > 0 && max(ind) > length(derivedConfs)) stop('MCMC configuration doesn\'t have that many derived quantity functions') + derivedConfs[ind] <<- NULL + if(print) printDerivedQuantities() + return(invisible(NULL)) + }, + + removeDerivedQuantity = function(...) { + ' +Alias for removeDerivedQuantities method +' + removeDerivedQuantities(...) + }, + + printDerivedQuantities = function(ind, type, displayNonScalars = FALSE, byType = FALSE) { + ' +Prints details of derived quantity functions. + +Arguments: + +ind: A numeric vector, used to specify the indices of the derived quantity functions to print. If omitted, then all derived quantity functions are printed. + +type: A character vector containing derived quantity function types. Only derived quantity functions with one of these specified types will be displayed. Regular expression matching is also used. + +displayNonScalars: A logical argument, specifying whether to display the values of non-scalar control list elements (default FALSE). + +byType: A logical argument, specifying whether a summary of the derived quantity functions should be printed, instead of the details of each individual function (default FALSE). +' + if(missing(ind)) { + ind <- seq_along(derivedConfs) + } + if(length(ind) > 0 && max(ind) > length(derivedConfs)) stop('MCMC configuration doesn\'t have that many derived quantity functions') + if(!missing(type)) { + if(!is.character(type)) stop('type argument must have type character') + ## find indices with 'name' matching anything in 'type' argument: + typeInd <- unique(unname(unlist(lapply(type, grep, x = lapply(derivedConfs, `[[`, 'name'))))) + ind <- intersect(ind, typeInd) + } + if(byType) { + printDerivedQuantitiesByType(ind) + return(invisible(NULL)) + } + makeSpaces <- if(length(ind) > 0) newSpacesFunction(max(ind)) else NULL + for(i in ind) { + info <- paste0('[', i, '] ', makeSpaces(i), derivedConfs[[i]]$toStr(displayNonScalars)) + cat(paste0(info, '\n')) + } + return(invisible(NULL)) + }, + + printDerivedQuantitiesByType = function(ind) { + if(length(ind) == 0) return(invisible(NULL)) + derivedTypes <- sapply(ind, function(i) derivedConfs[[i]]$name) + derivedTypesTable <- table(derivedTypes) + for(i in seq_along(derivedTypesTable)) { + name <- names(derivedTypesTable)[i] + num <- as.numeric(derivedTypesTable)[i] + cat(paste0('- ', name, ' (', num, ')\n')) + } + }, + + getDerivedQuantities = function(ind) { + ' +Returns a list of derivedConf objects. + +Arguments: + +ind: A numeric vector used to specify the indices of the derivedConf objects to return. If omitted, then all derivedConf objects in this MCMC configuration object are returned. +' + if(missing(ind)) ind <- seq_along(derivedConfs) + if(length(ind) > 0 && max(ind) > length(derivedConfs)) stop('MCMC configuration doesn\'t have that many derived quantity functions') + return(derivedConfs[ind]) + }, + + getDerivedQuantityDefinition = function(ind, print = FALSE) { + ' +Returns the nimbleFunction definition of aa derived quantity function. + +Arguments: + +ind: A numeric index used to specify the index of the derived quantity function definition to return. If more than one derived quantity function is specified, only the first is returned. + +Returns a list object, containing the setup function, run function, and additional member methods for the specified derived quantity nimbleFunction. +' + if(length(ind) > 1) { + messageIfVerbose(' [Note] More than one ved quantity functio specified, only returning the first.') + ind <- ind[1] + } + if((ind <= 0) || (ind > length(derivedConfs))) stop('Invalid derived quantity function specified') + if(print) printDerivedQuantities(ind) + def <- getDefinition(derivedConfs[[ind]]$derivedFunction) + return(def) + }, + addMonitors = function(..., ind = 1, print = TRUE) { ' Adds variables to the list of monitors. @@ -1065,7 +1319,7 @@ See the initialize() function ' if(isMvSamplesReady(ind)){ - messageIfVerbose(' [Note] Changing monitors, even though an MCMC has been built already. When compiling the MCMC, use resetFunctions = TRUE option.') + messageIfVerbose(' [Note] Changing monitors, even though an MCMC has been built already. When compiling the MCMC, use resetFunctions = TRUE option.') if(ind == 1) mvSamples1Conf <<- NULL if(ind == 2) @@ -1299,7 +1553,7 @@ See the initialize() function msg <- paste0(' [Warning] No samplers assigned for ', numUnsampled, ' node', sTag) if(includeConfGetUnsampledNodes) msg <- paste0(msg, ', use conf$getUnsampledNodes() for node name', sTag) msg <- paste0(msg, '.') - messageIfVerbose(msg) + message(msg) } }, @@ -1320,6 +1574,10 @@ See the initialize() function printMonitors() cat('===== Samplers =====\n') if(length(samplerConfs)) printSamplers(byType = TRUE) else cat('(no samplers assigned)\n') + if(length(derivedConfs)) { + cat('===== DerivedQ =====\n') + printDerivedQuantities(byType = TRUE) + } printComments(...) } ) @@ -1356,6 +1614,10 @@ See the initialize() function #'@param onlyRW A logical argument, with default value FALSE. If specified as TRUE, then Metropolis-Hastings random walk samplers (\link{sampler_RW}) will be assigned for all non-terminal continuous-valued nodes nodes. Discrete-valued nodes are assigned a slice sampler (\link{sampler_slice}), and terminal nodes are assigned a posterior_predictive sampler (\link{sampler_posterior_predictive}). #'@param onlySlice A logical argument, with default value FALSE. If specified as TRUE, then a slice sampler is assigned for all non-terminal nodes. Terminal nodes are still assigned a posterior_predictive sampler. #'@param multivariateNodesAsScalars A logical argument, with default value FALSE. If specified as TRUE, then non-terminal multivariate stochastic nodes will have scalar samplers assigned to each of the scalar components of the multivariate node. The default value of FALSE results in a single block sampler assigned to the entire multivariate node. Note, multivariate nodes appearing in conjugate relationships will be assigned the corresponding conjugate sampler (provided \code{useConjugacy == TRUE}), regardless of the value of this argument. +#'@param samplePredictiveNodes A logical argument. When TRUE, samplers will be assigned by default to update posterior predictive model nodes. When FALSE, no samplers will be assigned to predictive nodes. The default value of this argument is given by the nimble option \code{MCMCassignSamplersToPosteriorPredictiveNodes}, which itself has a default value of TRUE. +#'@param mean Character or logical argument. When a vector of node names is provided, the \code{mean} derived quantity function will be used to calculate the running mean for all node names specified. If \code{TRUE}, the running mean will be calculated for nodes specified in \code{monitors}. +#'@param variance Character or logical argument. When a vector of node names is provided, the \code{variance} derived quantity function will be used to calculate the running variance for all node names specified. If \code{TRUE}, the running variance will be calculated for nodes specified in \code{monitors}. +#'@param logProb When \code{TRUE}, the summed log-density of all stochastic model nodes (including data nodes) will be calculated and returned, using the \code{logProb} derived quantity function. When provided as a character vector, the individual log density of each node in this vector will be recorded. When provided as a list, each list element may contain one or mode node names, and separately for the node(s) in each element of the list, the summed log-density list will be calculated. In addition, the keyword \code{".all"} may also be provided in either the vector or list argument, which corresponds to the set of all stochastic model nodes (including data). #'@param enableWAIC A logical argument, specifying whether to enable WAIC calculations for the resulting MCMC algorithm. Defaults to the value of \code{nimbleOptions('MCMCenableWAIC')}, which in turn defaults to FALSE. Setting \code{nimbleOptions('enableWAIC' = TRUE)} will ensure that WAIC is enabled for all calls to \code{\link{configureMCMC}} and \code{\link{buildMCMC}}. #'@param controlWAIC A named list of inputs that control the behavior of the WAIC calculation. See \code{help(waic)}. #'@param print A logical argument, specifying whether to print the ordered list of default samplers. @@ -1371,6 +1633,8 @@ configureMCMC <- function(model, nodes, control = list(), useConjugacy = getNimbleOption('MCMCuseConjugacy'), onlyRW = FALSE, onlySlice = FALSE, multivariateNodesAsScalars = getNimbleOption('MCMCmultivariateNodesAsScalars'), + samplePredictiveNodes = getNimbleOption('MCMCassignSamplersToPosteriorPredictiveNodes'), + mean = FALSE, variance = FALSE, logProb = FALSE, enableWAIC = getNimbleOption('MCMCenableWAIC'), controlWAIC = list(), print = getNimbleOption('verbose'), autoBlock = FALSE, oldConf, @@ -1397,6 +1661,8 @@ configureMCMC <- function(model, nodes, control = list(), useConjugacy = useConjugacy, onlyRW = onlyRW, onlySlice = onlySlice, multivariateNodesAsScalars = multivariateNodesAsScalars, + samplePredictiveNodes = samplePredictiveNodes, + mean = mean, variance = variance, logProb = logProb, enableWAIC = enableWAIC, controlWAIC = controlWAIC, print = print, ...) return(invisible(thisConf)) @@ -1404,28 +1670,3 @@ configureMCMC <- function(model, nodes, control = list(), -# This is function which builds a new MCMCconf from an old MCMCconf -# This is required to be able to a new C-based MCMC without recompiling -makeNewConfFromOldConf <- function(oldMCMCconf){ - newMCMCconf <- configureMCMC(oldMCMCconf$model, nodes = NULL, print = FALSE) - newMCMCconf$monitors <- oldMCMCconf$monitors - newMCMCconf$monitors2 <- oldMCMCconf$monitors2 - newMCMCconf$thin <- oldMCMCconf$thin - newMCMCconf$thin2 <- oldMCMCconf$thin2 - newMCMCconf$samplerConfs <- oldMCMCconf$samplerConfs - newMCMCconf$samplerExecutionOrder <- oldMCMCconf$samplerExecutionOrder - newMCMCconf$controlDefaults <- oldMCMCconf$controlDefaults - ##newMCMCconf$namedSamplerLabelMaker <- oldMCMCconf$namedSamplerLabelMaker ## usage long since deprecated (Dec 2020) - newMCMCconf$mvSamples1Conf <- oldMCMCconf$mvSamples1Conf - newMCMCconf$mvSamples2Conf <- oldMCMCconf$mvSamples2Conf - return(newMCMCconf) -} - - -newSpacesFunction <- function(m) { - log10max <- floor(log10(m)) - function(i) paste0(rep(' ', log10max-floor(log10(i))), collapse = '') -} - - - diff --git a/packages/nimble/R/MCMC_conjugacy.R b/packages/nimble/R/MCMC_conjugacy.R index 541299b64..a4a22753d 100644 --- a/packages/nimble/R/MCMC_conjugacy.R +++ b/packages/nimble/R/MCMC_conjugacy.R @@ -653,7 +653,7 @@ conjugacyClass <- setRefClass( ## only if we're verifying conjugate posterior distributions: get initial targetValue, and modelLogProb -- model$getLogProb(calcNodes) - if(getNimbleOption('verifyConjugatePosteriors')) { + if(getNimbleOption('MCMCverifyConjugatePosteriors')) { functionBody$addCode({ modelLogProb0 <- model$getLogProb(calcNodes) origTargetValue <- model[[target]] @@ -672,7 +672,7 @@ conjugacyClass <- setRefClass( nimCopy(from = model, to = mvSaved, row = 1, nodes = calcNodes, logProb = TRUE) }, list(RPOSTERIORCALL = posteriorObject$rCallExpr)) ## only if we're verifying conjugate posterior distributions: figure out if conjugate posterior distribution is correct - if(getNimbleOption('verifyConjugatePosteriors')) { + if(getNimbleOption('MCMCverifyConjugatePosteriors')) { functionBody$addCode({ modelLogProb1 <- model$getLogProb(calcNodes) posteriorLogDensity0 <- DPOSTERIORCALL_ORIG diff --git a/packages/nimble/R/MCMC_derived.R b/packages/nimble/R/MCMC_derived.R new file mode 100644 index 000000000..55011e4d9 --- /dev/null +++ b/packages/nimble/R/MCMC_derived.R @@ -0,0 +1,457 @@ + +#################################################################### +### virtual nimbleFunction template for all derived quantities ##### +#################################################################### + +#' @rdname derived +#' @export +derived_BASE <- nimbleFunctionVirtual( + name = 'derived_BASE', + run = function(timesRan = double()) { }, + methods = list( + set_interval = function(newInterval = double()) { }, + before_chain = function(niter = double(), nburnin = double(), thin = double(1), chain = double()) { }, + after_chain = function() { }, + get_results = function() { returnType(double(2)) }, + get_names = function() { returnType(character(1)) }, + reset = function() { } + ), + methodControl = list( + before_chain = list(required = FALSE), + after_chain = list(required = FALSE), + reset = list(required = FALSE) + ) +) + + + +#################################################################### +### derived quantity: mean ######################################### +#################################################################### + +#' @rdname derived +#' @export +derived_mean <- nimbleFunction( + name = 'derived_mean', + contains = derived_BASE, + setup = function(model, mcmc, interval, control) { + ## control list extraction + nodes <- extractControlElement(control, 'nodes', defaultValue = character()) + recordingFrequency <- extractControlElement(control, 'recordingFrequency', defaultValue = 0) + ## node list generation + if(is.list(nodes)) nodes <- unlist(nodes) + nodes <- model$expandNodeNames(nodes) + ## names generation + names <- if(length(nodes) < 2) c(nodes,'','') else nodes ## vector + ## numeric value generation + nSamples <- 0 + nResults <- length(nodes) + vals <- numeric(max(nResults, 2)) ## vector + onlineMean <- numeric(max(nResults, 2)) ## vector + results <- array(0, c(1, nResults)) + }, + run = function(timesRan = double()) { + if(nResults == 0) return() + nSamples <<- nSamples + 1 + vals <<- values(model, nodes) + if(nSamples == 1) { + onlineMean <<- vals + } else { + onlineMean <<- onlineMean + (vals - onlineMean) / nSamples + } + if(recordingFrequency != 0 & timesRan %% recordingFrequency == 0) { + results[timesRan/recordingFrequency,] <<- onlineMean + } + }, + methods = list( + set_interval = function(newInterval = double()) { + interval <<- newInterval + }, + before_chain = function(niter = double(), nburnin = double(), thin = double(1), chain = double()) { + if(recordingFrequency == 0) { + nKeep <- 1 + } else { + nKeep <- floor(niter / (interval * recordingFrequency)) + } + results <<- nimArray(NA, c(nKeep, nResults)) + }, + after_chain = function() { + if(recordingFrequency == 0) { + results[1,] <<- onlineMean + } + }, + get_results = function() { + returnType(double(2)) + return(results) + }, + get_names = function() { + returnType(character(1)) + return(names) + }, + reset = function() { + nSamples <<- 0 + vals <<- nimNumeric(nResults, value = NA) + onlineMean <<- nimNumeric(nResults, value = NA) + } + ) +) + + + +#################################################################### +### derived quantity: variance ##################################### +#################################################################### + +#' @rdname derived +#' @export +derived_variance <- nimbleFunction( + name = 'derived_variance', + contains = derived_BASE, + setup = function(model, mcmc, interval, control) { + ## control list extraction + nodes <- extractControlElement(control, 'nodes', defaultValue = character()) + recordingFrequency <- extractControlElement(control, 'recordingFrequency', defaultValue = 0) + ## node list generation + if(is.list(nodes)) nodes <- unlist(nodes) + nodes <- model$expandNodeNames(nodes) + ## names generation + names <- if(length(nodes) < 2) c(nodes,'','') else nodes ## vector + ## numeric value generation + nSamples <- 0 + nResults <- length(nodes) + vals <- numeric(max(nResults, 2)) ## vector + prvMean <- numeric(max(nResults, 2)) ## vector + newMean <- numeric(max(nResults, 2)) ## vector + sumSqur <- numeric(max(nResults, 2)) ## vector + results <- array(0, c(1, nResults)) + }, + run = function(timesRan = double()) { + if(nResults == 0) return() + nSamples <<- nSamples + 1 + vals <<- values(model, nodes) + ## Welford's algorithm for stable online variance + if(nSamples == 1) { + newMean <<- vals + sumSqur <<- nimNumeric(nResults) + } else { + prvMean <<- newMean + newMean <<- prvMean + (vals - prvMean) / nSamples + sumSqur <<- sumSqur + (vals - prvMean) * (vals - newMean) + } + if(recordingFrequency != 0 & timesRan %% recordingFrequency == 0) { + if(nSamples == 1) { + results[timesRan/recordingFrequency,] <<- rep(NA, nResults) + } else { + results[timesRan/recordingFrequency,] <<- sumSqur / (nSamples-1) + } + } + }, + methods = list( + set_interval = function(newInterval = double()) { + interval <<- newInterval + }, + before_chain = function(niter = double(), nburnin = double(), thin = double(1), chain = double()) { + if(recordingFrequency == 0) { + nKeep <- 1 + } else { + nKeep <- floor(niter / (interval * recordingFrequency)) + } + results <<- nimArray(NA, c(nKeep, nResults)) + }, + after_chain = function() { + if(recordingFrequency == 0) { + if(nSamples == 1) { + results[1,] <<- rep(NA, nResults) + } else { + results[1,] <<- sumSqur / (nSamples-1) + } + } + }, + get_results = function() { + returnType(double(2)) + return(results) + }, + get_names = function() { + returnType(character(1)) + return(names) + }, + reset = function() { + nSamples <<- 0 + vals <<- nimNumeric(nResults, value = NA) + prvMean <<- nimNumeric(nResults, value = NA) + newMean <<- nimNumeric(nResults, value = NA) + sumSqur <<- nimNumeric(nResults, value = NA) + } + ) +) + + + +#################################################################### +### derived quantity: logProb ###################################### +#################################################################### + +getLogProb_virtual <- nimbleFunctionVirtual( + run = function() { returnType(double()) } +) + +getLogProbNF <- nimbleFunction( + contains = getLogProb_virtual, + setup = function(model, nodes) {}, + run = function() { + returnType(double()) + return(model$getLogProb(nodes)) + } +) + +#' @rdname derived +#' @export +derived_logProb <- nimbleFunction( + name = 'derived_logProb', + contains = derived_BASE, + setup = function(model, mcmc, interval, control) { + ## control list extraction + nodes <- extractControlElement(control, 'nodes', defaultValue = '.all') + silent <- extractControlElement(control, 'silent', defaultValue = FALSE) + ## node list generation + nodeList <- if(is.character(nodes)) { + as.list(unlist(lapply(nodes, function(x) if(identical(x,'.all')) '.all' else model$expandNodeNames(x)))) + } else nodes + allBool <- sapply(nodeList, function(x) identical(x, '.all')) + nodeList <- lapply(nodeList, function(x) if(identical(x,'.all')) model$getNodeNames(stochOnly=TRUE) else x) + ## names generation + if(is.list(nodes) && !is.null(names(nodes))) { + names <- names(nodes) + } else { + names <- character(length(nodeList)) + sumIndex <- 0 + for(i in seq_along(nodeList)) { + if(allBool[i]) { + names[i] <- '_all_nodes_' + next } + if(length(nodeList[[i]]) > 1) { + sumIndex <- sumIndex + 1 + names[i] <- paste0('sum', sumIndex) + next } + names[i] <- nodeList[[i]] + } + } + if(length(names) < 2) names <- c(names, '', '') ## vector + ## numeric value generation + nResults <- length(nodeList) + results <- array(0, c(1, nResults)) + ## nested function and function list definitions + getLogProbNFL <- nimbleFunctionList(getLogProb_virtual) + for(i in seq_along(nodeList)) getLogProbNFL[[i]] <- getLogProbNF(model, nodeList[[i]]) + ## checks + uniqueNodes <- unique(unlist(nodeList)) + missingInd <- sapply(uniqueNodes, function(x) length(model$expandNodeNames(x)) == 0) + missingNodes <- uniqueNodes[missingInd] + if(length(missingNodes) && !silent) + warning('logProb derived quantity function is using node names which are not in the model: ', + paste0(missingNodes, collapse=', '), call. = FALSE) + }, + run = function(timesRan = double()) { + if(nResults == 0) return() + for(i in 1:nResults) { + results[timesRan, i] <<- getLogProbNFL[[i]]$run() + } + }, + methods = list( + set_interval = function(newInterval = double()) { + interval <<- newInterval + }, + before_chain = function(niter = double(), nburnin = double(), thin = double(1), chain = double()) { + nKeep <- floor(niter / interval) + results <<- nimArray(NA, c(nKeep, nResults)) + }, + get_results = function() { + returnType(double(2)) + return(results) + }, + get_names = function() { + returnType(character(1)) + return(names) + }, + reset = function() { + results <<- nimArray(NA, c(1, nResults)) + } + ) +) + + + +#################################################################### +### derived quantity: predictive ################################### +#################################################################### + +#' @rdname derived +#' @export +derived_predictive <- nimbleFunction( + name = 'derived_predictive', + contains = derived_BASE, + setup = function(model, mcmc, interval, control) { + ## control list extraction + nodes <- extractControlElement(control, 'nodes', defaultValue = '.missing') + simNodes <- extractControlElement(control, 'simNodes', defaultValue = '.missing') + sort <- extractControlElement(control, 'sort', defaultValue = TRUE) + silent <- extractControlElement(control, 'silent', defaultValue = FALSE) + ## node list generation + ## all predictive stochastic nodes, and their deterministic dependencies + ppNodesAndDeps <- model$getDependencies(model$getNodeNames(predictiveOnly = TRUE), downstream = TRUE) + if(length(nodes) == 1 && nodes == '.missing') { + saveNodes <- model$getNodeNames(endOnly = TRUE, includeData = FALSE) + isEnd <- sapply(saveNodes, function(x) length(model$getDependencies(x, self = FALSE)) == 0) + saveNodes <- saveNodes[isEnd] + } else if(length(nodes) == 1 && nodes == '.all') { + ## this approach missed predicted deterministic nodes, which also have a stochastic dependency: + ##saveNodes <- unique(c( + ## ppNodesAndDeps, ## predictive stochastic nodes and deterministic dependencies + ## model$getNodeNames(endOnly = TRUE, includeData = FALSE) ## deterministic derived quantities + ##)) + allModelNodes <- model$getNodeNames() + allNodeDownstreamDeps <- lapply(allModelNodes, function(x) model$getDependencies(x, downstream = TRUE)) + haveDataDepsBool <- sapply(allNodeDownstreamDeps, function(x) any(model$isData(x))) + saveNodes <- allModelNodes[!haveDataDepsBool] + saveNodes <- model$topologicallySortNodes(saveNodes) + } else { + saveNodes <- model$expandNodeNames(nodes) + } + sampledNodesList <- lapply( + mcmc$samplerFunctions$contentsList, + function(x) { + nodes <- character() + if('target' %in% ls(x)) nodes <- c(nodes, x$target) + if('targetAsScalar' %in% ls(x)) nodes <- c(nodes, x$targetAsScalar) + if('simNodes' %in% ls(x)) nodes <- c(nodes, x$simNodes) + return(nodes) + }) + sampledNodes <- model$expandNodeNames(unlist(sampledNodesList)) + if(simNodes == '.missing') { + ## figuring this case out was not easy -DT May 2025 + upToDateNodes <- setdiff( ## nodes which should be kept up-to-date by the MCMC + model$getNodeNames(includePredictive = FALSE), ## all model nodes, excluding predictive stochastic nodes + ppNodesAndDeps ## predictive stochastic nodes and deterministic dependencies + ) + ## now we add any predictive nodes (and their deterministic dependencies), which might have samplers assigned + sampledNodeDeps <- model$getDependencies(sampledNodes) + sampledNodeDetermDeps <- sampledNodeDeps[model$isDeterm(sampledNodeDeps)] + upToDateNodes <- unique(c( + upToDateNodes, + sampledNodes, ## sampled stochastic nodes + sampledNodeDetermDeps ## deterministic dependencies of sampled nodes + )) + saveNodesParents <- model$getParents(saveNodes, self = TRUE, upstream = TRUE) ## everything upstream from (and including) saveNodes + simNodes <- setdiff(saveNodesParents, upToDateNodes) + } + if(sort) simNodes <- model$topologicallySortNodes(simNodes) + calcNodes <- model$getDependencies(simNodes) + mvSaved <- mcmc$mvSaved + ## names generation + names <- if(length(saveNodes) < 2) c(saveNodes,'','') else saveNodes ## vector + ## numeric value generation + nResults <- length(saveNodes) + results <- array(0, c(1, nResults)) + ## checks + otherwiseSampledSimNodes <- intersect(simNodes, sampledNodes) + if(length(otherwiseSampledSimNodes) & !silent) { + message(' [Warning] predictive derived quantity function is simulating nodes being updated by other MCMC samplers: ', + paste0(otherwiseSampledSimNodes, collapse=', ')) + } + stochSaveNodes <- saveNodes[model$isStoch(saveNodes)] + nonPPsaveNodes <- setdiff(stochSaveNodes, ppNodesAndDeps) + if(length(nonPPsaveNodes) & !silent) { + message(' [Warning] predictive derived quantity function is operating on non-posterior-predictive nodes: ', + paste0(nonPPsaveNodes, collapse=', ')) + } + }, + run = function(timesRan = double()) { + if(nResults == 0) return() + model$simulate(simNodes) + model$calculate(calcNodes) + results[timesRan,] <<- values(model, saveNodes) + nimCopy(from = model, to = mvSaved, row = 1, nodes = calcNodes, logProb = TRUE) + }, + methods = list( + set_interval = function(newInterval = double()) { + interval <<- newInterval + }, + before_chain = function(niter = double(), nburnin = double(), thin = double(1), chain = double()) { + nKeep <- floor(niter / interval) + results <<- nimArray(NA, c(nKeep, nResults)) + }, + get_results = function() { + returnType(double(2)) + return(results) + }, + get_names = function() { + returnType(character(1)) + return(names) + }, + reset = function() { + results <<- nimArray(0, c(1, nResults)) + } + ) +) + + + +#' MCMC Derived Quantities +#' +#' Details of the NIMBLE MCMC engine handles derived quantities, which are deterministic functions that can be calaculated and recorded after each MCMC sampling iteration. +#' +#' @section Mean and Variance +#' +#' The \code{mean} and \code{variance} derived quantity functions calculate the running mean and variance, respectively, for each node specified in the \code{nodes} argument. If added to an MCMC configuration object using the \code{addDerivedQuantity} method, then a value of the \code{interval} argument may also be provided to \code{addDerivedQuantity}. In that case, the value of \code{interval} specifies the number of MCMC iterations between calculations of the statistic. When the statistic is calculated, only the current value of each node is used to update the statistic. For example, if \code{interval} is 2, then every other MCMC iteration is used to calculate an updated value of the statistic. If no value of \code{interval} is provided as an argument to \code{addDerivedQuantity}, then the default value is the thinning interval \code{thin} of the MCMC. +#' +#' The \code{mean} and \code{variance} derived quantity functions both accept the following control list elements: +#' \itemize{ +#' \item nodes. The set of model nodes used for tracking the statistic. +#' \item recordingFrequency. The frequency (number of calculations of the statistic) afer which the value of the statistic is saved. For example, if \code{recordingFrequency} is 1, then the value of the statistic is saved after every update of its value. But if \code{recordingFrequency} is 10, then the value of the statistic is only saved after every tenth update of its value. The dafault value of \code{recordingFrequency} is 0, which corresponds to a special case: the value of the statistic is only recorded a single time, which is on the final iteration of the MCMC chain. +#' } +#' +#' @section Model Log-Densities +#' +#' The \code{logProb} derived quantity function calculates and records values of the log-density of individual nodes or (summed) groups of nodes. If added to an MCMC configuration object using the \code{addDerivedQuantity} method, then a value of the \code{interval} argument may also be provided to \code{addDerivedQuantity}. In that case, the value of \code{interval} specifies the number of MCMC iterations between recordings of the log-density values. For example, if \code{interval} is 2, then log-density values will be recorded upon every other MCMC iteration. If no value of \code{interval} is provided as an argument to \code{addDerivedQuantity}, then the default value is the thinning interval \code{thin} of the MCMC. +#' +#' The \code{logProb} derived quantity function accepts the following control list elements: +#' \itemize{ +#' \item nodes. The \code{nodes} argument determines the individual nodes, or (summed) groups of nodes, for recording log-density values. When provided as a character vector, the individual log-density of each node in this vector will be recorded. When provided as a list, each list element may contain one or mode node names, and separately for the node(s) in each element of the list, the summed log-density list will be calculated. In addition, the keyword \code{".all"} may also be provided in either the vector or list argument, which corresponds to the set of all stochastic model nodes (including data). +#' \item silent. By default, the \code{logProb} derived quantity function will issue a warning when the \code{nodes} argument includes node names which are not present in the model. This warning may be suppressed by setting \code{silent} to \code{TRUE}. +#' } +#' +#' @section Posterior Predictive Nodes and Derived Quantities +#' +#' The \code{predictive} derived quantity function simulates the values of posterior predictive nodes in the model and stores these simulated values. This may be useful when a model structure includes posterior predictive nodes (or deterministically defined posterior derived quantities), but for reasons of efficiency, these nodes may not undergo MCMC sampling. In such cases, the \code{predictive} derived quantity function may be assigned to these nodes, and when executed it will simulate new values for these nodes and record the simulated values. +#' +#' When added to an MCMC configuration object using the \code{addDerivedQuantity} method, a value of the \code{interval} argument may also be provided to \code{addDerivedQuantity}. In that case, the value of \code{interval} specifies the number of MCMC iterations between operations of the \code{predictive} function. For example, if \code{interval} is 2, then prediction and storing values takes place every other MCMC iteration. If no value of \code{interval} is provided as an argument to \code{addDerivedQuantity}, then the default value is the thinning interval \code{thin} of the MCMC. +#' +#' The \code{predictive} derived quantity function accepts the following control list elements: +#' \itemize{ +#' \item nodes. The \code{nodes} argument defines the predictive nodes which will have their values recorded and returned. The \code{predictive} function will automatically determine which nodes upstream of \code{nodes} need to be simulated, prior to storing the values of \code{nodes}. If omitted, the default set of \code{nodes} will be all (non-data) terminal model nodes. In addition, specifying the value \code{nodes = ".all"} will take \code{nodes} to be all predictive nodes in the model (all deterministic and stochastic nodes which have no downstream data dependencies). +#' \item simNodes. The \code{simNodes} specifies the set of model nodes which will be simulated, prior to recording the values of \code{nodes}. If ommited, which would be the usual case, \code{simNodes} will automatically represent the set of all model nodes which must be simulated in order to reach the predictive \code{nodes}. By providing \code{simNodes}, more fine-grained control of the simulation process is possible, including omitting simulation of some nodes. +#' \item sort. The \code{sort} argument determines whether the simulation of \code{simNodes} takes place in topological order. This argument has a default value of \code{TRUE}. When specified as \code{FALSE} the simulation of \code{simNodes} will take place in the order in which they were specified in the \code{simNodes} argument, which may not be in their natural order of dependency. +#' \item silent. By default, the \code{predictive} derived quantity function will issue a warning when the \code{simNodes} argument includes node names which are being updated by some sampler function in the MCMC. This warning may be suppressed by setting \code{silent} to \code{TRUE}. +#' } +#' +#' @name derived +#' +#' @aliases derived_mean derived_variance derived_logProb +#' +#' @examples +#' conf$addDerivedQuantity("mean", nodes = c("a", "b")) +#' +#' conf$addDerivedQuantity("mean", nodes = "theta", interval = 5) +#' +#' conf$addDerivedQuantity("variance", nodes = "x[1:4]", control = list(recordingFrequency = 10)) +#' +#' conf$addDerivedQuantity("logProb", nodes = c('alpha', 'beta')) +#' +#' conf <- configureMCMC(model, mean = 'a', variance = 'b', logProb = TRUE) +#' +#' @seealso \code{\link{configureMCMC}} \code{\link{addDerivedQuantity}} \code{\link{buildMCMC}} \code{\link{runMCMC}} \code{\link{nimbleMCMC}} +#' +#' @author Daniel Turek +#' +NULL + diff --git a/packages/nimble/R/MCMC_run.R b/packages/nimble/R/MCMC_run.R index 6411aeea0..9277eca30 100644 --- a/packages/nimble/R/MCMC_run.R +++ b/packages/nimble/R/MCMC_run.R @@ -28,6 +28,8 @@ #' #' @param summary Logical argument. When \code{TRUE}, summary statistics for the posterior samples of each parameter are also returned, for each MCMC chain. This may be returned in addition to the posterior samples themselves. Default value is \code{FALSE}. See details. #' +#' @param derivedQuantities Logical argument. When \code{TRUE}, a list containing the derivedQuantity results will be returned. This list has length equal to the number of derived quantity functions, and each element is the results matrix corresponding to one derived quantity function. In the case of multiple MCMC chains, a list (of length nchains) is returned consisting of one such list for each MCMC chain. +#' #' @param WAIC Logical argument. When \code{TRUE}, the WAIC (Watanabe, 2010) of the model is calculated and returned. Note that in order for the WAIC to be calculated, the \code{mcmc} object must have also been created with the argument `enableWAIC = TRUE`. If multiple chains are run, then a single WAIC value is calculated using the posterior samples from all chains. Default value is \code{FALSE}. See \code{help(waic)}. #' #' @param perChainWAIC Logical argument. When \code{TRUE} and multiple chains are run, the WAIC for each chain is returned as a means of helping assess the stability of the WAIC estimate. Default value is \code{FALSE}. @@ -93,12 +95,13 @@ runMCMC <- function(mcmc, samples = TRUE, samplesAsCodaMCMC = FALSE, summary = FALSE, + derivedQuantities = (mcmc$getNumDerived()>0) && getNimbleOption('MCMCreturnDerivedQuantities'), WAIC = FALSE, perChainWAIC = FALSE) { if(missing(mcmc)) stop('must provide a NIMBLE MCMC algorithm') if(!identical(nfGetDefVar(mcmc, 'name'), 'MCMC')) stop('mcmc argument must be a NIMBLE MCMC algorithm') if(!is.Cnf(mcmc)) messageIfVerbose(' [Warning] Running an uncompiled MCMC algorithm. Use compileNimble() for faster execution.') - if(!samples && !summary && !WAIC) stop('no output specified, use samples = TRUE, summary = TRUE, or WAIC = TRUE') + if(!samples && !summary && !derivedQuantities && !WAIC) stop('no output specified, use samples = TRUE, summary = TRUE, derivedQuantities = TRUE, or WAIC = TRUE') if(nchains < 1) stop('must have nchains > 0') if(!missing(inits)) { if(!is.function(inits) && !is.list(inits)) stop('inits must be a function, a list of initial values, or a list (of length nchains) of lists of inital values') @@ -114,6 +117,7 @@ runMCMC <- function(mcmc, hasMonitors2 <- length(if(is.Cnf(mcmc)) mcmc$Robject$monitors2 else mcmc$monitors2) > 0 samplesList <- vector('list', nchains); names(samplesList) <- paste0('chain', 1:nchains) samplesList2 <- vector('list', nchains); names(samplesList2) <- paste0('chain', 1:nchains) + derivedList <- vector('list', nchains); names(derivedList) <- paste0('chain', 1:nchains) thinToUseVec <- c(0, 0) thinToUseVec[1] <- if(!missing(thin)) thin else mcmc$thinFromConfVec[1] thinToUseVec[2] <- if(!missing(thin2)) thin2 else mcmc$thinFromConfVec[2] @@ -124,29 +128,48 @@ runMCMC <- function(mcmc, ## if(thinToUseVec[1] > 1 && nburnin > 0) message("runMCMC's handling of nburnin changed in nimble version 0.6-11. Previously, nburnin samples were discarded *post-thinning*. Now nburnin samples are discarded *pre-thinning*. The number of samples returned will be floor((niter-nburnin)/thin).") ## reinstate samplerExecutionOrder as a runtime argument, once we support non-scalar default values for runtime arguments: ##samplerExecutionOrderToUse <- if(!missing(samplerExecutionOrder)) samplerExecutionOrder else mcmc$samplerExecutionOrderFromConfPlusTwoZeros[mcmc$samplerExecutionOrderFromConfPlusTwoZeros>0] - for(i in 1:nchains) { - messageIfVerbose('running chain ', i, '...') + numDerived <- mcmc$getNumDerived() + for(chain in 1:nchains) { + messageIfVerbose('running chain ', chain, '...') if(is.numeric(setSeed)) { if(length(setSeed) != nchains) stop('setSeed argument has different length from nchains.') - set.seed(setSeed[i]) - } else if(setSeed) set.seed(i) + set.seed(setSeed[chain]) + } else if(setSeed) set.seed(chain) if(!missing(inits)) { if(is.function(inits)) { theseInits <- inits() } else if(is.list(inits) && (length(inits) > 0) && is.list(inits[[1]])) { - theseInits <- inits[[i]] + theseInits <- inits[[chain]] } else theseInits <- inits model$setInits(theseInits) } ##model$calculate() # shouldn't be necessary, since mcmc$run() includes call to my_initializeModel$run() - mcmc$run(niter, nburnin = nburnin, thin = thinToUseVec[1], thin2 = thinToUseVec[2], progressBar = progressBar, resetWAIC = ifelse(i == 1, TRUE, FALSE), chain = i) #, samplerExecutionOrder = samplerExecutionOrderToUse) + mcmc$run(niter, + nburnin = nburnin, + thin = thinToUseVec[1], + thin2 = thinToUseVec[2], + progressBar = progressBar, + resetWAIC = ifelse(chain == 1, TRUE, FALSE), + chain = chain + ## samplerExecutionOrder = samplerExecutionOrderToUse + ) tmp <- as.matrix(mcmc$mvSamples) if(!is.null(tmp)) - samplesList[[i]] <- tmp + samplesList[[chain]] <- tmp if(hasMonitors2) { tmp <- as.matrix(mcmc$mvSamples2) if(!is.null(tmp)) - samplesList2[[i]] <- tmp + samplesList2[[chain]] <- tmp + } + if(derivedQuantities && (numDerived>0)) { + tempList <- lapply(1:numDerived, function(i) mcmc$getDerivedQuantityResults(i)) + for(i in 1:numDerived) { + theseNames <- mcmc$getDerivedQuantityNames(i) + if(ncol(tempList[[i]]) > 0) + colnames(tempList[[i]]) <- theseNames[1:ncol(tempList[[i]])] + } + names(tempList) <- mcmc$derivedTypes + derivedList[[chain]] <- tempList } } if(WAIC) { @@ -184,6 +207,8 @@ runMCMC <- function(mcmc, if(hasMonitors2) if(length(samplesList2)) samplesList2 <- samplesList2[[1]] else samplesList2 <- NULL ## returns matrix when nchains = 1 + if(derivedQuantities) + derivedList <- derivedList[[1]] else derivedList <- NULL ## returns matrix when nchains = 1 } if(summary) { if(nchains == 1) { @@ -206,6 +231,7 @@ runMCMC <- function(mcmc, if(summary) retList$summary <- summaryObject if(WAIC) retList$WAIC <- WAICvalue if(perChainWAIC) retList$perChainWAIC <- perChainWAICvalue + if(derivedQuantities) retList$derived <- derivedList if(length(retList) == 1) retList <- retList[[1]] return(retList) } @@ -250,7 +276,13 @@ runMCMC <- function(mcmc, #' @param samplesAsCodaMCMC Logical argument. If \code{TRUE}, then a \code{coda} \code{mcmc} object is returned instead of an R matrix of samples, or when \code{nchains > 1} a \code{coda} \code{mcmc.list} object is returned containing \code{nchains} \code{mcmc} objects. This argument is only used when \code{samples} is \code{TRUE}. Default value is \code{FALSE}. See details. #' #' @param summary Logical argument. When \code{TRUE}, summary statistics for the posterior samples of each parameter are also returned, for each MCMC chain. This may be returned in addition to the posterior samples themselves. Default value is \code{FALSE}. See details. -#'z +#' +#' @param mean Character or logical argument. When a vector of node names is provided, the \code{mean} derived quantity function will be used to calculate the running mean for all node names specified. If \code{TRUE}, the running mean will be calculated for nodes specified in \code{monitors}. +#' +#' @param variance Character or logical argument. When a vector of node names is provided, the \code{variance} derived quantity function will be used to calculate the running variance for all node names specified. If \code{TRUE}, the running variance will be calculated for nodes specified in \code{monitors}. +#' +#' @param logProb When \code{TRUE}, the summed log-density of all stochastic model nodes (including data nodes) will be calculated and returned, using the \code{logProb} derived quantity function. When provided as a character vector, the individual log density of each node in this vector will be recorded. When provided as a list, each list element may contain one or mode node names, and separately for the node(s) in each element of the list, the summed log-density list will be calculated. In addition, the keyword \code{".all"} may also be provided in either the vector or list argument, which corresponds to the set of all stochastic model nodes (including data). +#' #' @param WAIC Logical argument. When \code{TRUE}, the WAIC (Watanabe, 2010) of the model is calculated and returned. If multiple chains are run, then a single WAIC value is calculated using the posterior samples from all chains. Default value is \code{FALSE}. Note that the version of WAIC used is the default WAIC conditional on random effects/latent states and without any grouping of data nodes. See \code{help(waic)} for more details. If a different version of WAIC is desired, do not use \code{nimbleMCMC}. Instead, specify the \code{controlWAIC} argument to \code{configureMCMC} or \code{buildMCMC}, and then use \code{runMCMC}. #' #' @return A list is returned with named elements depending on the arguments passed to \code{nimbleMCMC}, unless only one among samples, summary, and WAIC are requested, in which case only that element is returned. These elements may include \code{samples}, \code{summary}, and \code{WAIC}. When \code{nchains = 1}, posterior samples are returned as a single matrix, and summary statistics as a single matrix. When \code{nchains > 1}, posterior samples are returned as a list of matrices, one matrix for each chain, and summary statistics are returned as a list containing \code{nchains+1} matrices: one matrix corresponding to each chain, and the final element providing a summary of all chains, combined. If \code{samplesAsCodaMCMC} is \code{TRUE}, then posterior samples are provided as \code{coda} \code{mcmc} and \code{mcmc.list} objects. When \code{WAIC} is \code{TRUE}, a WAIC summary object is returned. @@ -313,6 +345,9 @@ nimbleMCMC <- function(code, samples = TRUE, samplesAsCodaMCMC = FALSE, summary = FALSE, + mean = FALSE, + variance = FALSE, + logProb = FALSE, WAIC = FALSE, userEnv = parent.frame()) { #### process 'code' argument, to accept a filename, or a function @@ -329,7 +364,7 @@ nimbleMCMC <- function(code, if(!samples && !summary && !WAIC) stop('no output specified, use samples = TRUE, summary = TRUE, or WAIC = TRUE') if(!missing(code) && inherits(code, 'modelBaseClass')) model <- code ## let's handle it, if model object is provided as un-named first argument to nimbleMCMC Rmodel <- mcmc_createModelObject(model, inits, nchains, setSeed, code, constants, data, dimensions, check, userEnv = userEnv) - conf <- configureMCMC(Rmodel, monitors = monitors, thin = thin, enableWAIC = WAIC, print = FALSE) + conf <- configureMCMC(Rmodel, monitors = monitors, thin = thin, mean = mean, variance = variance, logProb = logProb, enableWAIC = WAIC, print = FALSE) Rmcmc <- buildMCMC(conf) compiledList <- compileNimble(Rmodel, Rmcmc) ## only one compileNimble() call Cmcmc <- compiledList$Rmcmc diff --git a/packages/nimble/R/MCMC_utils.R b/packages/nimble/R/MCMC_utils.R index 07cf58c96..90d4e758d 100644 --- a/packages/nimble/R/MCMC_utils.R +++ b/packages/nimble/R/MCMC_utils.R @@ -266,7 +266,8 @@ mcmc_generateControlListArgument <- function(control, controlDefaults) { -mcmc_listContentsToStr <- function(ls, displayControlDefaults=FALSE, displayNonScalars=FALSE, displayConjugateDependencies=FALSE) { +mcmc_listContentsToStr <- function(ls, displayControlDefaults=FALSE, displayNonScalars=FALSE, displayConjugateDependencies=FALSE, + removeCfunctions = TRUE, removeLengthZero = TRUE) { ##if(any(unlist(lapply(ls, is.function)))) warning('probably provided wrong type of function argument') if(!displayConjugateDependencies) { if(grepl('^conjugate_d', names(ls)[1])) ls <- ls[1] ## for conjugate samplers, remove all 'dep_dnorm', etc, control elements (don't print them!) @@ -287,7 +288,7 @@ mcmc_listContentsToStr <- function(ls, displayControlDefaults=FALSE, displayNonS for(i in seq_along(ls)) { controlName <- names(ls)[i] controlValue <- ls[[i]] - if(length(controlValue) == 0) next ## remove length 0 + if(length(controlValue) == 0 && removeLengthZero) next ## remove length 0 ##if(!displayControlDefaults) ## if(controlName %in% names(defaultOptions)) ## skip default control values ## if(identical(controlValue, defaultOptions[[controlName]])) next @@ -305,7 +306,7 @@ mcmc_listContentsToStr <- function(ls, displayControlDefaults=FALSE, displayNonS ##if(length(ls2) == 1) ## str <- paste0(str, ', default') str <- gsub('\"', '', str) - str <- gsub('c\\((.*?)\\)', '\\1', str) + if(removeCfunctions) str <- gsub('c\\((.*?)\\)', '\\1', str) return(str) } @@ -551,6 +552,37 @@ mcmc_checkTargetAD <- function(model, targetNodes, samplerType) { } +# This is function which builds a new MCMCconf from an old MCMCconf +# This is required to be able to a new C-based MCMC without recompiling +makeNewConfFromOldConf <- function(oldMCMCconf){ + newMCMCconf <- configureMCMC(oldMCMCconf$model, nodes = NULL, print = FALSE) + newMCMCconf$monitors <- oldMCMCconf$monitors + newMCMCconf$monitors2 <- oldMCMCconf$monitors2 + newMCMCconf$thin <- oldMCMCconf$thin + newMCMCconf$thin2 <- oldMCMCconf$thin2 + newMCMCconf$samplerConfs <- oldMCMCconf$samplerConfs + newMCMCconf$samplerExecutionOrder <- oldMCMCconf$samplerExecutionOrder + newMCMCconf$controlDefaults <- oldMCMCconf$controlDefaults + ##newMCMCconf$namedSamplerLabelMaker <- oldMCMCconf$namedSamplerLabelMaker ## usage long since deprecated (Dec 2020) + newMCMCconf$mvSamples1Conf <- oldMCMCconf$mvSamples1Conf + newMCMCconf$mvSamples2Conf <- oldMCMCconf$mvSamples2Conf + return(newMCMCconf) +} + + +newSpacesFunction <- function(m) { + log10max <- floor(log10(m)) + function(i) paste0(rep(' ', log10max-floor(log10(i))), collapse = '') +} + + +getBaseClassName <- function(nf) { + baseName <- environment(environment(nf)$contains)$className + if(is.null(baseName)) warning('cannot find base class name') + return(baseName) +} + + diff --git a/packages/nimble/R/cppDefs_nimbleFunction.R b/packages/nimble/R/cppDefs_nimbleFunction.R index 15fe92136..73b228eaf 100644 --- a/packages/nimble/R/cppDefs_nimbleFunction.R +++ b/packages/nimble/R/cppDefs_nimbleFunction.R @@ -771,7 +771,7 @@ modifyForAD_indexingBracket <- function(code, symTab, workEnv) { } ## modifyForAD_issuePowWarning <- function(code, symTab, workEnv) { -## message(" [Note] Operator `pow` may cause derivative problems with negative arguments. If the exponent is guaranteed to be an integer, use `pow_int` instead.") +## message(" [Note] Operator `pow` may cause derivative problems with negative arguments. If the exponent is guaranteed to be an integer, use `pow_int` instead.") ## invisible(NULL) ## } @@ -935,12 +935,12 @@ modifyForAD_getDerivs_wrapper <- function(code, symTab, workEnv) { modifyForAD_nfMethod <- function(code, symTab, workEnv) { if(code$args[[1]]$name != "cppPointerDereference") - message(" [Note] In modifyForAD_nfMethod, was expecting cppPointerDereference. There must be another case that needs implementation.") + message(" [Note] In modifyForAD_nfMethod, was expecting cppPointerDereference. There must be another case that needs implementation.") objName <- code$args[[1]]$args[[1]]$name NFsymObj <- workEnv$RsymTab$getSymbolObject(objName, TRUE) methodName <- code$args[[2]] if(!is.character(methodName)) - message(" [Note] In modifyForAD_nfMethod, was expecting method name to be character. There must be another case that needs implementation.") + message(" [Note] In modifyForAD_nfMethod, was expecting method name to be character. There must be another case that needs implementation.") methodSymObj <- NFsymObj$nfProc$compileInfos[[methodName]]$newLocalSymTab$getSymbolObject(methodName, TRUE) buildDerivs <- methodSymObj$nfMethodRCobj$buildDerivs if(!is.null(buildDerivs)) diff --git a/packages/nimble/R/options.R b/packages/nimble/R/options.R index c8405fe9f..c017a653b 100644 --- a/packages/nimble/R/options.R +++ b/packages/nimble/R/options.R @@ -196,13 +196,12 @@ nimOptimMethod("bobyqa", verbose = TRUE, verboseErrors = FALSE, + showCompilerOutput = FALSE, + ## verifies the correct posterior is created for any conjugate samplers, at run-time. ## if this option is changed, then congugate sampler functions can be rebuilt using: ## buildConjugateSamplerFunctions() - verifyConjugatePosteriors = FALSE, - - showCompilerOutput = FALSE, - + MCMCverifyConjugatePosteriors = FALSE, MCMCprogressBar = TRUE, MCMCsaveHistory = FALSE, MCMCmultivariateNodesAsScalars = FALSE, @@ -211,11 +210,14 @@ nimOptimMethod("bobyqa", MCMCorderPriorSamplesSamplersFirst = TRUE, MCMCorderPosteriorPredictiveSamplersLast = TRUE, MCMCusePredictiveDependenciesInCalculations = FALSE, - MCMCusePosteriorPredictiveSampler = TRUE, + MCMCassignSamplersToPosteriorPredictiveNodes = TRUE, ## whether any samplers are assigned (by default) to PP nodes + MCMCusePosteriorPredictiveSampler = TRUE, ## for PP nodes being sampled, use post_pred (or otherwise RW, etc) MCMCwarnUnsampledStochasticNodes = TRUE, MCMCRJcheckHyperparam = TRUE, MCMCenableWAIC = FALSE, MCMCuseBarkerAsDefaultMV = FALSE, + MCMCreturnDerivedQuantities = TRUE, + parameterTransformWarnUserDists = TRUE, useClearCompiledInADTesting = TRUE, unsupportedDerivativeHandling = 'error', # default is error, other options are 'warn' and 'ignore'. Handled in updateADproxyModelMethods in cppDefs_nimbleFunction.R @@ -244,7 +246,7 @@ setNimbleOption <- function(name, value) { #' @export #' @return The value of the option. #' @examples -#' getNimbleOption('verifyConjugatePosteriors') +#' getNimbleOption('MCMCverifyConjugatePosteriors') getNimbleOption <- function(x) { option <- try(get(x, envir = .nimbleOptions), silent = TRUE) if(inherits(option, 'try-error')) @@ -274,7 +276,7 @@ getNimbleOption <- function(x) { #' #' @examples #' # Set one option: -#' nimbleOptions(verifyConjugatePosteriors = FALSE) +#' nimbleOptions(MCMCverifyConjugatePosteriors = FALSE) #' #' # Compactly print all options: #' str(nimbleOptions(), max.level = 1) diff --git a/packages/nimble/tests/testthat/test-derived-quantities.R b/packages/nimble/tests/testthat/test-derived-quantities.R new file mode 100644 index 000000000..27230f367 --- /dev/null +++ b/packages/nimble/tests/testthat/test-derived-quantities.R @@ -0,0 +1,624 @@ +source(system.file(file.path('tests', 'testthat', 'test_utils.R'), package = 'nimble')) + +code <- nimbleCode({ + b[1] ~ dnorm(0, sd = 2) + b[2] ~ dnorm(0, sd = 2) + sigma ~ dunif(0, 2) + for (i in 1:n){ + mu[i] <- b[1] + b[2]*x[i] + y[i] ~ dnorm(mu[i], sd = sigma) + } +}) + +test_that("mean derived parameter", { + + # Single derived parameter (mean) + set.seed(1) + mod <- nimbleModel(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10))) + conf <- configureMCMC(mod, print=FALSE) + conf$addDerivedQuantity("mean", nodes="b[1]") + expect_is(conf$derivedConfs[[1]], "derivedConf") + + mcmc <- buildMCMC(conf) + out <- runMCMC(mcmc, niter=10) + + expect_is(out, "list") + expect_equal(names(out), c("samples", "derived")) + expect_equivalent(out$derived$mean[1,1], + mean(out$samples[,"b[1]"])) + + # Means for two parameters + conf <- configureMCMC(mod, print=FALSE) + conf$addDerivedQuantity("mean", nodes=c("b[1]", "b[2]")) + mcmc <- buildMCMC(conf) + set.seed(1) + out <- runMCMC(mcmc, niter=10) + expect_equal(out$derived$mean[1,], + colMeans(out$samples[,1:2])) + + # Using range + conf <- configureMCMC(mod, print=FALSE) + conf$addDerivedQuantity("mean", nodes=c("b[1:2]")) + mcmc <- buildMCMC(conf) + set.seed(1) + out <- runMCMC(mcmc, niter=10) + expect_equal(out$derived$mean[1,], + colMeans(out$samples[,1:2])) + + # Combination of two separate parameters + conf <- configureMCMC(mod, print=FALSE) + conf$addDerivedQuantity("mean", nodes=list(c("b[1]", "sigma"))) + mcmc <- buildMCMC(conf) + set.seed(1) + out <- runMCMC(mcmc, niter=10) + expect_equivalent(round(out$derived$mean[1,], 4), round(apply(out$samples[,c("b[1]", "sigma")], 2, mean),4)) + + # What if the parameter doesn't exist in the model? + # result is empty matrix + conf <- configureMCMC(mod, print=FALSE) + conf$addDerivedQuantity("mean", nodes=c("alpha"), interval=11) + mcmc <- buildMCMC(conf) + set.seed(1) + out <- runMCMC(mcmc, niter=10) + expect_equal(out$derived$mean, matrix(0, 1, 0)) + + # If parameter index is too big + conf <- configureMCMC(mod, print=FALSE) + conf$addDerivedQuantity("mean", nodes=c("b[3]"), interval=11) + mcmc <- buildMCMC(conf) + set.seed(1) + out <- runMCMC(mcmc, niter=10) + expect_equal(out$derived$mean, matrix(0, 1, 0)) +}) + +test_that("interval and thinning rate for derived parameter", { + + # Changing interval + # every other mcmc iteration will be used to calculate mean (starting with 2) + mod <- nimbleModel(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10))) + conf <- configureMCMC(mod, print=FALSE) + conf$addDerivedQuantity("mean", nodes=c("b[1]"), interval=2) + mcmc <- buildMCMC(conf) + set.seed(1) + out <- runMCMC(mcmc, niter=10) + expect_equivalent(mean(out$samples[c(2,4,6,8,10),1]), + out$derived$mean[1,1]) + + # Try interval of 3 + conf <- configureMCMC(mod, print=FALSE) + conf$addDerivedQuantity("mean", nodes=c("b[1]"), interval=3) + mcmc <- buildMCMC(conf) + set.seed(1) + out <- runMCMC(mcmc, niter=10) + expect_equivalent(mean(out$samples[c(3,6,9),1]), + out$derived$mean[1,1]) + + # What happens if interval is larger than sample? + # output is NA + conf <- configureMCMC(mod, print=FALSE) + conf$addDerivedQuantity("mean", nodes=c("b[1]"), interval=11) + mcmc <- buildMCMC(conf) + set.seed(1) + out <- runMCMC(mcmc, niter=10) + expect_true(is.na(out$derived$mean[1,1])) + + # How does interval interact with thinning rate? + # First create unthinned samples for reference + set.seed(1) + mod <- nimbleModel(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10))) + conf <- configureMCMC(mod, print=FALSE) + mcmc <- buildMCMC(conf) + out_unthinned <- runMCMC(mcmc, niter=10) + + set.seed(1) + mod <- nimbleModel(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10))) + conf <- configureMCMC(mod, thin = 2, print=FALSE) + # Without setting interval; this should set interval to the same as thin automatically + conf$addDerivedQuantity("mean", nodes=c("b[1]")) + mcmc <- buildMCMC(conf) + out <- runMCMC(mcmc, niter=10) + expect_equivalent(mean(out_unthinned[c(2,4,6,8,10), 1]), + out$derived$mean[1,1]) + + # What if we manually set interval to 1? + # then all samples (thinned and unthinned) are used in mean calculation + set.seed(1) + mod <- nimbleModel(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10))) + conf <- configureMCMC(mod, thin = 2, print=FALSE) + # Without setting interval; this should set interval to the same as thin automatically + conf$addDerivedQuantity("mean", nodes=c("b[1]"), interval=1) + mcmc <- buildMCMC(conf) + out <- runMCMC(mcmc, niter=10) + expect_equivalent(mean(out_unthinned[, 1]), + out$derived$mean[1,1]) + + # Setting different thin and interval rates + set.seed(1) + mod <- nimbleModel(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10))) + conf <- configureMCMC(mod, thin = 2, print=FALSE) + # Without setting interval; this should set interval to the same as thin automatically + conf$addDerivedQuantity("mean", nodes=c("b[1]"), interval=3) + mcmc <- buildMCMC(conf) + out <- runMCMC(mcmc, niter=10) + expect_equivalent(mean(out_unthinned[c(3,6,9), 1]), + out$derived$mean[1,1]) +}) + +test_that("recording frequency for derived parameters", { + # Testing recording frequency + # The default 0 (calculated once at the end) + set.seed(1) + mod <- nimbleModel(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10))) + conf <- configureMCMC(mod, thin = 1, print=FALSE) + conf$addDerivedQuantity("mean", nodes=c("b[1]"), interval=1, recordingFrequency=0) + mcmc <- buildMCMC(conf) + out_ref <- runMCMC(mcmc, niter=10) + expect_equivalent(mean(out_ref$samples[,1]), out_ref$derived$mean[1,1]) + + # frequency 1 (cumulative mean) + set.seed(1) + mod <- nimbleModel(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10))) + conf <- configureMCMC(mod, thin = 1, print=FALSE) + conf$addDerivedQuantity("mean", nodes=c("b[1]"), interval=1, recordingFrequency=1) + mcmc <- buildMCMC(conf) + out <- runMCMC(mcmc, niter=10) + expect_equivalent(cumsum(out$samples[,1]) / seq_along(out$samples[,1]), + out$derived$mean[,1]) + + # Vary the thinning rate + # This should not affect derived quantity results when we set interval and frequency to 1 + set.seed(1) + mod <- nimbleModel(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10))) + conf <- configureMCMC(mod, thin = 2, print=FALSE) + conf$addDerivedQuantity("mean", nodes=c("b[1]"), interval=1, recordingFrequency=1) + mcmc <- buildMCMC(conf) + out <- runMCMC(mcmc, niter=10) + # Comparing to the reference output above with no thinning + expect_equivalent(cumsum(out_ref$samples[,1]) / seq_along(out_ref$samples[,1]), + out$derived$mean[,1]) + + # frequency = 2 + # cumulative mean but calculate it half as often + set.seed(1) + mod <- nimbleModel(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10))) + conf <- configureMCMC(mod, thin = 1, print=FALSE) + conf$addDerivedQuantity("mean", nodes=c("b[1]"), interval=1, recordingFrequency=2) + mcmc <- buildMCMC(conf) + out <- runMCMC(mcmc, niter=10) + expect_equivalent(c(mean(out$samples[1:2]), mean(out$samples[1:4]), mean(out$samples[1:6]), + mean(out$samples[1:8]), mean(out$samples[1:10])), + out$derived$mean[,1]) + + # interval = 2 and frequency = 2 + # only use every other sample and only calculate after 2 recorded samples + # result should be length 2 + set.seed(1) + mod <- nimbleModel(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10))) + conf <- configureMCMC(mod, thin = 1, print=FALSE) + conf$addDerivedQuantity("mean", nodes=c("b[1]"), interval=2, recordingFrequency=2) + mcmc <- buildMCMC(conf) + out <- runMCMC(mcmc, niter=10) + expect_equivalent(out$derived$mean, + c(mean(out$samples[c(2,4),1]), mean(out$samples[c(2,4,6,8)]))) + + # What if recording frequency is greater than number of iterations? + # Output is numeric(0) + conf <- configureMCMC(mod, thin = 1, print=FALSE) + conf$addDerivedQuantity("mean", nodes=c("b[1]"), interval=1, recordingFrequency=11) + mcmc <- buildMCMC(conf) + out <- runMCMC(mcmc, niter=10) + expect_equivalent(out$derived$mean, numeric(0)) + + # Same thing in this case when combo of interval + recording frequency means + # there aren't any derived values + conf <- configureMCMC(mod, thin = 1, print=FALSE) + conf$addDerivedQuantity("mean", nodes=c("b[1]"), interval=6, recordingFrequency=2) + mcmc <- buildMCMC(conf) + out <- runMCMC(mcmc, niter=10) + expect_equivalent(out$derived$mean, numeric(0)) +}) + +test_that("mean derived parameter with compiled mcmc", { + # Check compiled version works as expected + set.seed(1) + mod <- nimbleModel(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10))) + conf <- configureMCMC(mod, print=FALSE) + conf$addDerivedQuantity("mean", nodes="b[1]") + mcmc <- buildMCMC(conf) + Cmod <- compileNimble(mod) + Cmcmc <- compileNimble(mcmc, project=mod) + out <- runMCMC(Cmcmc, niter=10) + expect_equivalent(mean(out$samples[,1]), out$derived$mean[1,1]) +}) + +test_that("variance derived parameter", { + # Variance + # fewer tests here since it should have same behavior as mean + set.seed(1) + mod <- nimbleModel(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10))) + conf <- configureMCMC(mod, print=FALSE) + conf$addDerivedQuantity("variance", nodes="b[1]") + mcmc <- buildMCMC(conf) + out <- runMCMC(mcmc, niter=10) + expect_equivalent(var(out$samples[,1]), out$derived$variance[1,1]) + + # Check compiled version + set.seed(1) + mod <- nimbleModel(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10))) + conf <- configureMCMC(mod, print=FALSE) + conf$addDerivedQuantity("variance", nodes="b[1]") + mcmc <- buildMCMC(conf) + Cmod <- compileNimble(mod) + Cmcmc <- compileNimble(mcmc, project=mod) + out <- runMCMC(Cmcmc, niter=10) + expect_equivalent(var(out$samples[,1]), out$derived$variance[1,1]) + + # Combination of both mean and variance + set.seed(1) + mod <- nimbleModel(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10))) + conf <- configureMCMC(mod, print=FALSE) + conf$addDerivedQuantity("mean", nodes="b[1]") + conf$addDerivedQuantity("variance", nodes="b[1]") + mcmc <- buildMCMC(conf) + out <- runMCMC(mcmc, niter=10) + expect_equivalent(out$derived$mean, mean(out$samples[,1])) + expect_equivalent(out$derived$variance, var(out$samples[,1])) +}) + +test_that("logProb derived quantity", { + + set.seed(1) + inits <- list(b=c(0, 0), sigma=1) + mod <- nimbleModel(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10)), + inits=inits) + conf <- configureMCMC(mod, print=FALSE) + + # Single node + conf$addDerivedQuantity("logProb", nodes="b[1]") + mcmc <- buildMCMC(conf) + + out <- runMCMC(mcmc, niter=10) + expect_equal(dim(out$derived$logProb), c(10,1)) + expect_equivalent(out$derived$logProb[10,"b[1]"], mod$logProb_b[1]) + + # non-scalar parameter + set.seed(1) + mod <- nimbleModel(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10)), + inits=list(b=c(0,0), sigma=1)) + conf <- configureMCMC(mod, print=FALSE) + conf$addDerivedQuantity("logProb", nodes="b") + mcmc <- buildMCMC(conf) + out <- runMCMC(mcmc, niter=10) + expect_equal(colnames(out$derived$logProb), c("b[1]", "b[2]")) + expect_equivalent(out$derived$logProb[10,], mod$logProb_b) + + # Multiple parameters + set.seed(1) + mod <- nimbleModel(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10)), + inits=list(b=c(0,0), sigma=1)) + conf <- configureMCMC(mod, print=FALSE) + conf$addDerivedQuantity("logProb", nodes=c("b[1]", "sigma")) + mcmc <- buildMCMC(conf) + out <- runMCMC(mcmc, niter=10) + expect_equal(colnames(out$derived$logProb), c("b[1]", "sigma")) + expect_equivalent(out$derived$logProb[10,], c(mod$logProb_b[1], mod$logProb_sigma)) + + # Multiple parameters including non-scalars + set.seed(1) + mod <- nimbleModel(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10)), + inits=list(b=c(0,0), sigma=1)) + conf <- configureMCMC(mod, print=FALSE) + conf$addDerivedQuantity("logProb", nodes=c("b", "sigma")) + mcmc <- buildMCMC(conf) + out <- runMCMC(mcmc, niter=10) + expect_equal(colnames(out$derived$logProb), c("b[1]", "b[2]", "sigma")) + expect_equivalent(out$derived$logProb[10,], c(mod$logProb_b, mod$logProb_sigma)) + + # All parameters + set.seed(1) + mod <- nimbleModel(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10)), + inits=list(b=c(0,0), sigma=1)) + conf <- configureMCMC(mod, print=FALSE) + # Here .all means the sum of logprobs of all parameters + conf$addDerivedQuantity("logProb", nodes=".all") + mcmc <- buildMCMC(conf) + out <- runMCMC(mcmc, niter=10) + expect_equal(colnames(out$derived$logProb), c("_all_nodes_")) + expect_equivalent(out$derived$logProb[10,1], sum(c(mod$logProb_b, mod$logProb_sigma, mod$logProb_y))) + + # Functions of parameters + set.seed(1) + mod <- nimbleModel(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10)), + inits=list(b=c(0,0), sigma=1)) + conf <- configureMCMC(mod, print=FALSE) + # Here taking sum logprobs of b + conf$addDerivedQuantity("logProb", nodes=list("b", "sigma")) + mcmc <- buildMCMC(conf) + out <- runMCMC(mcmc, niter=10) + expect_equal(colnames(out$derived$logProb), c("b", "sigma")) + expect_equivalent(out$derived$logProb[10,], c(sum(mod$logProb_b), mod$logProb_sigma)) + + set.seed(1) + mod <- nimbleModel(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10)), + inits=list(b=c(0,0), sigma=1)) + conf <- configureMCMC(mod, print=FALSE) + # Here taking sum of logprobs of b[1] and sigma + conf$addDerivedQuantity("logProb", nodes=list(c("b[1]", "sigma"))) + mcmc <- buildMCMC(conf) + out <- runMCMC(mcmc, niter=10) + expect_equal(colnames(out$derived$logProb), "sum1") + expect_equivalent(out$derived$logProb[10,1], sum(c(mod$logProb_b[1], mod$logProb_sigma))) + + # Create name for sum + set.seed(1) + mod <- nimbleModel(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10)), + inits=list(b=c(0,0), sigma=1)) + conf <- configureMCMC(mod, print=FALSE) + # Here taking sum of logprobs of b[1] and sigma + conf$addDerivedQuantity("logProb", nodes=list(b_and_sigma=c("b[1]", "sigma"))) + mcmc <- buildMCMC(conf) + out <- runMCMC(mcmc, niter=10) + expect_equal(colnames(out$derived$logProb), "b_and_sigma") + expect_equivalent(out$derived$logProb[10,1], sum(c(mod$logProb_b[1], mod$logProb_sigma))) + + # What happens if parameter doesn't exist? An empty matrix is returned + set.seed(1) + mod <- nimbleModel(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10)), + inits=list(b=c(0,0), sigma=1)) + conf <- configureMCMC(mod, print=FALSE) + # alpha not in model + conf$addDerivedQuantity("logProb", nodes="alpha") + mcmc <- buildMCMC(conf) + out <- runMCMC(mcmc, niter=10) + expect_equal(out$derived$logProb, matrix(0, 10, 0)) + + # What if an impossible index range is requested? + set.seed(1) + mod <- nimbleModel(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10)), + inits=list(b=c(0,0), sigma=1)) + conf <- configureMCMC(mod, print=FALSE) + # alpha not in model + conf$addDerivedQuantity("logProb", nodes="b[1:3]") + mcmc <- buildMCMC(conf) + out <- runMCMC(mcmc, niter=10) + # The extra element of b is ignored + expect_equal(colnames(out$derived$logProb), c("b[1]", "b[2]")) + + # Request only the non-existing index - result is empty matrix + set.seed(1) + mod <- nimbleModel(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10)), + inits=list(b=c(0,0), sigma=1)) + conf <- configureMCMC(mod, print=FALSE) + # alpha not in model + conf$addDerivedQuantity("logProb", nodes="b[3]") + mcmc <- buildMCMC(conf) + out <- runMCMC(mcmc, niter=10) + expect_equal(out$derived$logProb, matrix(0, 10, 0)) + + # Change interval + set.seed(1) + mod <- nimbleModel(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10)), + inits=list(b=c(0,0), sigma=1)) + conf <- configureMCMC(mod, print=FALSE) + # Here taking sum of logprobs of b[1] and sigma + conf$addDerivedQuantity("logProb", nodes="b[1]", interval=5) + mcmc <- buildMCMC(conf) + out <- runMCMC(mcmc, niter=10) + expect_equal(dim(out$derived$logProb), c(2,1)) +}) + +test_that("predictive derived quantity", { + + code2 <- nimbleCode({ + b[1] ~ dnorm(0, sd = 2) + b[2] ~ dnorm(0, sd = 2) + sigma ~ dunif(0, 2) + for (i in 1:n){ + mu[i] <- b[1] + b[2]*x[i] + y[i] ~ dnorm(mu[i], sd = sigma) + } + + bsq <- b[1]^2 # derived quantity + z ~ dnorm(bsq, sd = sigma) # posterior predictive node + }) + + # Simple example + set.seed(1) + mod <- nimbleModel(code2, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10)), + inits=list(b=c(0,0), sigma=1)) + conf <- configureMCMC(mod, print=FALSE, samplePredictiveNodes=FALSE) + + conf$addDerivedQuantity("predictive", nodes=c("bsq", "z")) + mcmc <- buildMCMC(conf) + out <- runMCMC(mcmc, niter=10) + expect_equal(dim(out$derived$predictive), c(10,2)) + expect_equal(colnames(out$derived$predictive), c("bsq", "z")) + + # Only the predictive node + set.seed(1) + mod <- nimbleModel(code2, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10)), + inits=list(b=c(0,0), sigma=1)) + conf <- configureMCMC(mod, print=FALSE, samplePredictiveNodes=FALSE) + + conf$addDerivedQuantity("predictive", nodes="z") + mcmc <- buildMCMC(conf) + out <- runMCMC(mcmc, niter=10) + expect_equal(dim(out$derived$predictive), c(10,1)) + expect_equal(colnames(out$derived$predictive), "z") + + # Change the interval + set.seed(1) + mod <- nimbleModel(code2, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10)), + inits=list(b=c(0,0), sigma=1)) + conf <- configureMCMC(mod, print=FALSE, samplePredictiveNodes=FALSE) + + conf$addDerivedQuantity("predictive", nodes="z", interval=5) + mcmc <- buildMCMC(conf) + out <- runMCMC(mcmc, niter=10) + expect_equal(dim(out$derived$predictive), c(2,1)) + expect_equal(colnames(out$derived$predictive), "z") + + # Set a non-predictive node + set.seed(1) + mod <- nimbleModel(code2, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10)), + inits=list(b=c(0,0), sigma=1)) + conf <- configureMCMC(mod, print=FALSE, samplePredictiveNodes=FALSE) + + conf$addDerivedQuantity("predictive", nodes=c("b[1]", "z")) + # warning message here + expect_message(mcmc <- buildMCMC(conf)) + out <- runMCMC(mcmc, niter=10) + expect_equal(dim(out$derived$predictive), c(10, 2)) + expect_equal(colnames(out$derived$predictive), c("b[1]", "z")) + + # Set .all + set.seed(1) + mod <- nimbleModel(code2, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10)), + inits=list(b=c(0,0), sigma=1)) + conf <- configureMCMC(mod, print=FALSE, samplePredictiveNodes=FALSE) + + # Check that samples of both bsq and z are returned when setting .all + conf$addDerivedQuantity("predictive", nodes=".all") + mcmc <- buildMCMC(conf) + out <- runMCMC(mcmc, niter=10) + expect_equal(dim(out$derived$predictive), c(10,2)) + expect_equal(colnames(out$derived$predictive), c("bsq", "z")) + + # Non-existing parameter - empty matrix is returned + set.seed(1) + mod <- nimbleModel(code2, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10)), + inits=list(b=c(0,0), sigma=1)) + conf <- configureMCMC(mod, print=FALSE, samplePredictiveNodes=FALSE) + conf$addDerivedQuantity("predictive", nodes="alpha") + mcmc <- buildMCMC(conf) + out <- runMCMC(mcmc, niter=10) + expect_equal(out$derived$predictive, matrix(0, 10, 0)) + +}) + +test_that("user-facing mcmc configuration methods", { + + set.seed(1) + inits <- list(b=c(0, 0), sigma=1) + mod <- nimbleModel(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10)), + inits=inits) + conf <- configureMCMC(mod, print=FALSE) + + # Add node + conf$addDerivedQuantity("logProb", nodes="b[1]") + mcmc <- buildMCMC(conf) + + # Print + co <- capture.output(conf$printDerivedQuantities()) + expect_equal(co, "[1] derived quantity: logProb, execution interval: thin, nodes: b[1]") + + # Get + conf$addDerivedQuantity("logProb", nodes="sigma") + dq <- conf$getDerivedQuantities() + expect_equal(length(dq), 2) + expect_is(dq[[1]], "derivedConf") + + dqf <- conf$getDerivedQuantityDefinition(1) + expect_is(dqf, "list") + + # Remove nodes + conf2 <- conf + conf2$removeDerivedQuantities() + co <- capture.output(conf2$printDerivedQuantities()) + expect_equal(co, character(0)) + + # Remove nodes + conf2 <- conf + conf2$removeDerivedQuantity() + co <- capture.output(conf2$printDerivedQuantities()) + expect_equal(co, character(0)) +}) + +test_that("configureMCMC and nimbleMCMC derived quantity interface", { + + set.seed(1) + inits <- list(b=c(0, 0), sigma=1) + mod <- nimbleModel(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10)), + inits=inits) + conf <- configureMCMC(mod, print=FALSE, mean="b[1]", logProb="sigma") + co <- capture.output(conf$printDerivedQuantities()) + expect_equal(co, c("[1] derived quantity: mean, execution interval: thin, nodes: b[1]", + "[2] derived quantity: logProb, execution interval: thin, nodes: sigma")) + + # logprob=TRUE + set.seed(1) + inits <- list(b=c(0, 0), sigma=1) + mod <- nimbleModel(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10)), + inits=inits) + conf <- configureMCMC(mod, print=FALSE, logProb=TRUE) + co <- capture.output(conf$printDerivedQuantities()) + expect_equal(co, "[1] derived quantity: logProb, execution interval: thin, nodes: .all") + + # nimbleMCMC + mod <- nimbleMCMC(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10)), + mean="b[1]", logProb="sigma", nchain=1, niter=10, nburnin=0) + expect_equal(names(mod$derived), c("mean", "logProb")) + + mod <- nimbleMCMC(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10)), + logProb=TRUE, nchain=1, niter=10, nburnin=0) + expect_equal(dim(mod$derived$logProb), c(10,1)) +}) + +test_that("custom derived quantity nimbleFunction", { + # multiplies the specified nodes by some value multValue + test_derived <- nimbleFunction( + name = 'test_derived', + contains = derived_BASE, ## all derived functions must contain (inherit from) the derived_BASE class + setup = function(model, mcmc, interval, control) { + nodes <- extractControlElement(control, 'nodes', defaultValue = character()) + multValue <- extractControlElement(control, 'multValue', defaultValue = 2) + ## node list generation + nodes <- model$expandNodeNames(nodes) + ## names generation + names <- if(length(nodes) < 2) c(nodes,'','') else nodes ## vector + ## numeric value generation + nResults <- length(nodes) + vals <- numeric(max(nResults, 2)) + results <- array(0, c(1, nResults)) + }, + run = function(timesRan = double()) { + if(nResults == 0) return() + vals <<- values(model, nodes) + results[1,] <<- vals * multValue + }, + methods = list( + set_interval = function(newInterval = double()) { + interval <<- newInterval + }, + get_results = function() { + returnType(double(2)) + return(results) + }, + get_names = function() { + returnType(character(1)) + return(names) + } + ) + ) + temporarilyAssignInGlobalEnv(test_derived) + + set.seed(1) + inits <- list(b=c(0, 0), sigma=1) + mod <- nimbleModel(code, constants=list(n=10, x=rnorm(10)), data = list(y=rnorm(10)), + inits=inits) + conf <- configureMCMC(mod, print=FALSE) + conf$addDerivedQuantity("test_derived", nodes=c("b[1]", "sigma"), multValue=2) + mcmc <- buildMCMC(conf) + out <- runMCMC(mcmc, niter=10) + + expect_equal(out$samples[10,c(1,3)]*2, out$derived$test_derived[1,]) + + # Check in C++ + Cmod <- compileNimble(mod) + Cmcmc <- compileNimble(mcmc, project=mod) + out <- runMCMC(Cmcmc, niter=10) + expect_equal(out$samples[10,c(1,3)]*2, out$derived$test_derived[1,]) + + # check error when custom derived quantity function not defined + expect_error(conf$addDerivedQuantity("fake", nodes=c("b[1]", "sigma"))) +})