diff --git a/packages/nimble/DESCRIPTION b/packages/nimble/DESCRIPTION index 7e8bc8b00..e0c55e57a 100644 --- a/packages/nimble/DESCRIPTION +++ b/packages/nimble/DESCRIPTION @@ -151,4 +151,4 @@ Collate: nimble-package.r QuadratureGrids.R zzz.R -RoxygenNote: 7.3.1 +RoxygenNote: 7.3.2 diff --git a/packages/nimble/R/miscAlgorithms.r b/packages/nimble/R/miscAlgorithms.r new file mode 100644 index 000000000..6999c5d19 --- /dev/null +++ b/packages/nimble/R/miscAlgorithms.r @@ -0,0 +1,253 @@ + +## Not exported function to get approximate number of iterations. +## Based on Sherlock 2021 Appendix A and code from: +## https://github.com/ChrisGSherlock/expQ/blob/master/rexpQ.cpp +## get_m(). Switching back to qpois for now. Will turn use when +## buildDerivs = TRUE in next iteration. +# getNPrec = nimbleFunction( + # run = function(rho = double(), prec = double(0, default= 1e-8)) { + + # logprec <- log(prec) + # pi <- 3.14159265 + # logroot2pi <- 0.5*log(2*pi) + + # mhi <- rho + (-logprec)/3.0*(1.0+sqrt(1.0+18.0*rho/(-logprec))) + # A <- 2*rho*(1.0-mhi/rho+mhi/rho*log(mhi/rho)) + # mlo <- rho+sqrt(2*rho)*sqrt(-logroot2pi-logprec-1.5*log(A)+log(A-1)) + # returnType(integer()) + # return(ceiling((mlo+mhi)/2)) + # } #, buildDerivs = TRUE +# ) + +## Not exported function to bound the spectrum of a square matrix: +## Based on eigenDisc from RTMB, which allows for the +## generalization of the Sherlock 2021 algorithm for more than just +## generator matrices. +spectrumBound <- nimbleFunction( + run = function(A = double(2)){ + q <- dim(A) + Min <- Inf + Max <- -Inf + for( i in 1:q[1] ){ + centre <- A[i,i] + radius <- sum(abs(A[i,])) - abs(centre) + minC <- centre - radius + maxC <- centre + radius + if(minC < Min) Min <- minC + if(maxC > Max) Max <- maxC + } + ans <- nimNumeric(value = 0, length = 2) + returnType(double(1)) + ans[1] <- (Max + Min)/2 + ans[2] <- (Max - Min)/2 + return(ans) + }#, buildDerivs = TRUE +) + +#' Matrix Exponential times a vector +#' +#' Compute the combined term expm(A) %*% v +#' to avoid a full matrix exponentiation. +#' +#' @name expAv +#' +#' @param A Square matrix. +#' @param v vector to multiply by the matrix exponential exp(A) %*% v. +#' @param tol level of accuracy required (default = 1e-8). +#' @param rescaleFreq How frequently should the terms be scaled to avoid underflow/overflow (default = 10). +#' @param Nmax Maximum number of iterations to compute (default = 10000). +#' @param sparse (logical) specify if the matrix may be sparse and to do sparse computation (default = TRUE). +#' @author Paul van Dam-Bates +#' @details For large matrix exponentials it is much more efficient to compute exp(A) %*% v, than to actually compute the entire matrix exponential. +#' +#' This function follows the function `expAv` from the R package RTMB (Kristensen, 2025), and theory outlined in Sherlock (2021). It is developed for working with continuous times +#' Markov chains. If using the matrix exponential to create a transition probability matrix in a HMM context just once, +#' this function may be slower than the one time call to compute the full matrix exponentiation. If a full matrix exponentiation is required, refer to +#' `expm` to compute. Choosing sparse = TRUE will check which matrix A values are non-zero and do sparse linear algebra. +#' Note that for computation efficiency matrix uniformization is done by A* = A + rho I, where rho = max(abs(diag(A))), see Algorithm 2' in Sherlock (2021). +#' When the row sums of the matrix are not zero, then uniformization is not done, and the number of iterations to reach tolerance are approximated based on the +#' a bound of the spectrum, similar to `RTMB' (Kristensen, 2025). +#' +#' @return \code{expAv} gives a vector that is ans = exp(A) %*% v. +#' +#' @references +#' Sherlock, C. (2021). Direct statistical inference for finite Markov jump processes via the matrix exponential. Computational Statistics, 36(4), 2863-2887. +#' +#' Kristensen K (2025). _RTMB: 'R' Bindings for 'TMB'_. R package version 1.7, commit 6bd7a16403ccb4d3fc13ff7526827540bf27b352, +#' . +#' +#' @examples +#' A <- rbind(c(-1, 0.25, 0.75), c(0, -2, 2), c(0.25, 0.25, -0.5)) +#' v <- c(0.35, 0.25, 0.1) +#' expAv(A, v) +NULL + +#' @rdname expAv +#' @export +expAv <- nimbleFunction( + run = function(A = double(2), v = double(1), + tol = double(0, default = 1e-8), + rescaleFreq = double(0, default = 10), + Nmax = integer(0, default = 10000), + sparse = logical(0, default = TRUE)){ + q <- dim(A) + if(q[2] != q[1]) stop("A must be a square matrix.") + if(length(v) != q[1]) stop("Length of v must be equal to the number of columns of A.") + + ## Only centre, perform uniformization, if no diag > 0. + CR <- spectrumBound(A) + diag(A) <- diag(A) - CR[1] + ans <- v + term <- v + log_scale <- 0 + m <- 1L + Niter <- 1L + # R <- ADbreak(CR[2]) + # Niter <- getNPrec(CR[2], tol) + Niter <- qpois(tol, CR[2], lower.tail = FALSE) + if(Niter > Nmax) Niter <- Nmax + + ## Check and build sparse matrix if sparsity exists. + ## Compressed sparse col format + if(sparse){ + m <- prod(q) + colptr <- nimNumeric(value = 0, length = q[2]+1) + rowindex <- nimNumeric(value = 0, length = m) + values <- nimNumeric(value = 0, length = m) + i <- 1L; j <- 1L; k <- 1L + colptr[1] <- 1 + for( j in 1:q[2] ){ + for( i in 1:q[1] ){ + if(A[i,j] != 0){ + rowindex[k] <- i + values[k] <- A[i,j] + k <- k+1 + } + } + colptr[j+1] <- k + } + m <- k-1 + rowindex <- rowindex[1:m] + values <- values[1:m] + } + check <- rescaleFreq + for(n in 1:Niter) { + if(sparse){ + tmp <- nimNumeric(0, length = q[1]) ## Row Vector + ## Sparse column matrix vector multiplication: term <- A %*% term/n + ## for each column, add a[rowindices,j] * term[j] + for( j in 1:q[2] ){ + if(colptr[j+1] > colptr[j]) + tmp[rowindex[colptr[j]:(colptr[j+1]-1)]] <- tmp[rowindex[colptr[j]:(colptr[j+1]-1)]] + term[j]*values[colptr[j]:(colptr[j+1]-1)] + } + term <- tmp / n + }else{ + term <- (A %*% asCol(term/n))[,1] + } + ans <- ans + term + if( check == n ) { + s <- sum(abs(ans)) + term <- term / s + ans <- ans / s + log_scale <- log_scale + log(s) + check <- check + rescaleFreq + } + } + if(Niter == Nmax) + cat(" Warning: Nmax in `expAv` is less than N required for the specified tolerance. Consider increasing Nmax.\n") + + ans <- exp(CR[1] + log_scale) * ans + returnType(double(1)) + return(ans) + },# buildDerivs = TRUE +) + +#' Matrix Exponential +#' +#' Compute the the matrix exponential expm(A) +#' by scaling and squaring. +#' +#' @name expm +#' +#' @param A Square matrix. +#' @param tol level of accuracy required (default = 1e-8). +#' @author Paul van Dam-Bates +#' +#' @details +#' This function follows the scaling and squaring algorithm from Sherlock (2021), except that we compute the full +#' matrix exponential. It differs from the standard Taylor scaling and squaring algorithm reviewed by Ruiz et al (2016) and found in common texts, +#' by doing uniformization if the matrix is a generator matrix from a continuous time Markov chain. If using the matrix exponential to create a +#' transition probability matrix in a HMM context just once, this function may be efficient. If dimension is large, we recommend avoiding the +#' matrix exponential and using `expAv` instead. Note that for computation efficiency, when the columns are non-positive, matrix uniformization +#' is done by A* = A + rho I, where rho = max(abs(diag(A))). When the row sums of the matrix are not zero, then uniformization is not done, +#' and the number of iterations to reach tolerance are approximated based on the a bound of the spectrum, similar to `RTMB' (Kristensen, 2025). +#' +#' @return \code{expm} gives a matrix that is ans = exp(A). +#' +#' @references +#' Sherlock, C. (2021). Direct statistical inference for finite Markov jump processes via the matrix exponential. Computational Statistics, 36(4), 2863-2887. +#' +#' Ruiz, P., Sastre, J., Ibáñez, J., & Defez, E. (2016). High performance computing of the matrix exponential. +#' Journal of computational and applied mathematics, 291, 370-379. +#' +#' Kristensen K (2025). _RTMB: 'R' Bindings for 'TMB'_. R package version 1.7, commit 6bd7a16403ccb4d3fc13ff7526827540bf27b352, +#' . +#' +#' @examples +#' A <- rbind(c(-1, 0.25, 0.75), c(0, -2, 2), c(0.25, 0.25, -0.5)) +#' Lambda <- diag(c(0.25, 0.1, 0)) +#' expm((A-Lambda)*2.5) +NULL + + +#' @rdname expm +#' @export +expm <- nimbleFunction( + run = function(A = double(2), tol = double(0, default = 1e-8)){ + returnType(double(2)) + + if(dim(A)[1] != dim(A)[2]) stop("`A' must be a square matrix") + # if(dim(tol)[1] > 1) stop("`tol' must be a scalar.") + + ## Only centre, perform uniformization, if no diag > 0. + CR <- spectrumBound(A) + Niter <- 1L + C <- CR[1] + R <- CR[2] + # Niter <- getNPrec(R, tol) + + n <- dim(A)[1] + log2 <- log(2) + shat <- ceiling((log(R)+log(log2))/log2) + smin <- max(c(0, shat)) + smax <- smin+6; + rhosmall <- R/exp(smin*log(2)) + opt <- Inf + for( k in smin:smax ){ + rhosmall <- rhosmall/2 + # test <- getNPrec(rhosmall, tol) + k + test <- qpois(tol, rhosmall, lower.tail = FALSE) + k + if(test < opt){ + opt <- test + s <- k + } + } + twos <- 2^s + Msmall <- (A - C*diag(n))/twos + # m <- getNPrec(R/twos, tol) + m <- qpois(tol, R/twos, lower.tail = FALSE) + m <- max(2, m) + + expMsmall <- diag(n) + Msmall + expMpro <- Msmall + for( i in 2:m ) { + expMpro <- expMpro %*% Msmall / i + expMsmall <- expMsmall + expMpro + } + + expMsmall <- expMsmall * exp(C/twos) + + for( i in 1:s ) expMsmall <- expMsmall %*% expMsmall + return(expMsmall) + }#, buildDerivs = TRUE +) diff --git a/packages/nimble/man/MCMCconf-class.Rd b/packages/nimble/man/MCMCconf-class.Rd index 2cfa9fa82..4c00c67dc 100644 --- a/packages/nimble/man/MCMCconf-class.Rd +++ b/packages/nimble/man/MCMCconf-class.Rd @@ -41,6 +41,35 @@ See documentation below for method initialize() for details of creating an MCMCc ... )}}{For internal use. Adds default MCMC samplers to the specified nodes.} +\item{\code{addDerivedQuantity( + 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.} + \item{\code{addMonitors(..., ind = 1, print = TRUE)}}{Adds variables to the list of monitors. Arguments: @@ -67,6 +96,14 @@ Details: See the initialize() function } +\item{\code{addOneDerivedQuantity( + thisDerivedName, + derivedFunction, + interval, + thisControlList, + print +)}}{For internal use only} + \item{\code{addOneSampler( thisSamplerName, samplerFunction, @@ -135,6 +172,20 @@ Invisibly returns a list of the current sampler configurations, which are sample The second usage of 'multivariateNodesAsScalars' occurs when 'default' is TRUE, and therefore samplers are assigned according to the default logic of configureMCMC, which is further controlled by the arguments 'useConjugacy', 'onlyRW', 'onlySlice' and 'multivariateNodesAsScalars'. In this second usage, if 'multivariateNodesAsScalars' is TRUE, then multivariate nodes will be decomposed into their scalar components, and separate samplers assigned to each scalar element. Note, however, that multivariate nodes appearing in conjugate relationships will still be assigned the corresponding conjugate sampler (provided 'useConjugacy' is TRUE), regardless of the value of this argument. If 'multivariateNodesAsScalars' is FALSE, then a single multivarate sampler will be assigned to update each multivariate node. The default value of this argument can be controlled using the nimble option 'MCMCmultivariateNodesAsScalars'.} +\item{\code{getDerivedQuantities(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.} + +\item{\code{getDerivedQuantityDefinition(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.} + \item{\code{getMonitors()}}{Returns a character vector of the current monitors Details: @@ -179,6 +230,10 @@ ind: A numeric vector or character vector. A numeric vector may be used to spec onlyRW = FALSE, onlySlice = FALSE, multivariateNodesAsScalars = getNimbleOption("MCMCmultivariateNodesAsScalars"), + samplePredictiveNodes = getNimbleOption("MCMCassignSamplersToPosteriorPredictiveNodes"), + mean = FALSE, + variance = FALSE, + logProb = FALSE, enableWAIC = getNimbleOption("MCMCenableWAIC"), controlWAIC = list(), print = TRUE, @@ -218,6 +273,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)`. @@ -226,6 +289,18 @@ print: A logical argument specifying whether to print the montiors and samplers. ...: Additional named control list elements for default samplers, or additional arguments to be passed to the autoBlock function when autoBlock = TRUE.} +\item{\code{printDerivedQuantities(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).} + \item{\code{printMonitors()}}{Prints all current monitors and monitors2 Details: @@ -250,7 +325,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). @@ -260,6 +335,18 @@ executionOrder: A logical argument, specifying whether to print the sampler func byType: A logical argument, specifying whether the nodes being sampled should be printed, sorted and organized according to the type of sampler (the sampling algorithm) which is acting on the nodes (default FALSE).} +\item{\code{removeDerivedQuantities(..., 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).} + +\item{\code{removeDerivedQuantity(...)}}{Alias for removeDerivedQuantities method} + \item{\code{removeSampler(...)}}{Alias for removeSamplers method} \item{\code{removeSamplers(..., ind, print = FALSE)}}{Removes one or more samplers from an MCMCconf object. diff --git a/packages/nimble/man/MultivariateNormal.Rd b/packages/nimble/man/MultivariateNormal.Rd index d8df5d7a1..bc6b209cf 100644 --- a/packages/nimble/man/MultivariateNormal.Rd +++ b/packages/nimble/man/MultivariateNormal.Rd @@ -2,10 +2,16 @@ % Please edit documentation in R/distributions_implementations.R \name{MultivariateNormal} \alias{MultivariateNormal} +\alias{dmnorm_inv_ld} +\alias{rmnorm_inv_ld} \alias{dmnorm_chol} \alias{rmnorm_chol} \title{The Multivariate Normal Distribution} \usage{ +dmnorm_inv_ld(x, mean, mat, inv_ld, prec_param = TRUE, log = FALSE) + +rmnorm_inv_ld(n = 1, mean, mat, inv_ld, prec_param = TRUE) + dmnorm_chol(x, mean, cholesky, prec_param = TRUE, log = FALSE) rmnorm_chol(n = 1, mean, cholesky, prec_param = TRUE) @@ -15,13 +21,13 @@ rmnorm_chol(n = 1, mean, cholesky, prec_param = TRUE) \item{mean}{vector of values giving the mean of the distribution.} -\item{cholesky}{upper-triangular Cholesky factor of either the precision matrix (when \code{prec_param} is TRUE) or covariance matrix (otherwise).} - \item{prec_param}{logical; if TRUE the Cholesky factor is that of the precision matrix; otherwise, of the covariance matrix.} \item{log}{logical; if TRUE, probability density is returned on the log scale.} \item{n}{number of observations (only \code{n=1} is handled currently).} + +\item{cholesky}{upper-triangular Cholesky factor of either the precision matrix (when \code{prec_param} is TRUE) or covariance matrix (otherwise).} } \value{ \code{dmnorm_chol} gives the density and \code{rmnorm_chol} generates random deviates. diff --git a/packages/nimble/man/buildMCMC.Rd b/packages/nimble/man/buildMCMC.Rd index 0f202bd58..08b5f29bc 100644 --- a/packages/nimble/man/buildMCMC.Rd +++ b/packages/nimble/man/buildMCMC.Rd @@ -37,7 +37,9 @@ The uncompiled \code{run} function will have arguments: \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). \code{progressBar}: Boolean specifying whether to display a progress bar during MCMC execution (default = TRUE). The progress bar can be permanently disabled by setting the system option \code{nimbleOptions(MCMCprogressBar = FALSE)}. diff --git a/packages/nimble/man/buildMacro.Rd b/packages/nimble/man/buildMacro.Rd new file mode 100644 index 000000000..ad70752d3 --- /dev/null +++ b/packages/nimble/man/buildMacro.Rd @@ -0,0 +1,202 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/BUGS_macros.R +\name{buildMacro} +\alias{buildMacro} +\title{EXPERIMENTAL: Turn a function into a model macro} +\usage{ +buildMacro(fun, use3pieces = TRUE, unpackArgs = TRUE) +} +\arguments{ +\item{fun}{A function written to construct new lines of model code (see below).} + +\item{use3pieces}{logical indicating whether the arguments from the input +line be split into pieces for the LHS (left-hand side), RHS +(right-hand side, possibly further split depending on +\code{unpackArgs}), and \code{stoch} (\code{TRUE} if the line uses a +\code{~} and \code{FALSE} otherwise). (default = \code{TRUE})} + +\item{unpackArgs}{logical indicating whether arguments be passed as a list +(\code{FALSE}) or as separate arguments (\code{TRUE}). (default = \code{TRUE})} +} +\value{ +A list of class \code{model_macro} with one element called \code{process}, +which contains the macro function suitable for use by \code{nimbleModel}. +} +\description{ +A model macro expands one line of code in a nimbleModel into one or +more new lines. This supports compact programming by defining +re-usable modules. \code{buildMacro} takes as input a +function that constructs new lines of model code from the original +line of code. It returns a function suitable for internal use by +\code{nimbleModel} that arranges arguments for input function. Macros +are an experimental feature and are available only after setting +\code{nimbleOptions(enableMacros = TRUE)}. +} +\details{ +The arguments \code{use3pieces} and \code{unpackArgs} +indicate how \code{fun} expects to have arguments arranged from an +input line of code (processed by \code{nimbleModel}). + +Consider the defaults \code{use3pieces = TRUE} and \code{unpackArgs = +TRUE}, for a macro called \code{macro1}. In this case, the line of +model code \code{x ~ macro1(arg1 = z[1:10], arg2 = "hello")} will be +passed to \code{fun} as \code{fun(stoch = TRUE, LHS = x, arg1 = +z[1:10], arg2 = "hello")}. + +If \code{use3pieces = TRUE} but \code{unpackArgs = FALSE}, then the +RHS will be passed as is, without unpacking its arguments into +separate arguments to \code{fun}. In this case, \code{x ~ macro1(arg1 += z[1:10], arg2 = "hello")} will be passed to \code{fun} as +\code{fun(stoch = TRUE, LHS = x, RHS = macro1(arg1 = z[1:10], arg2 = +"hello"))}. + +If \code{use3pieces = FALSE} and \code{unpackArgs = FALSE}, the entire +line of code is passed as a single object. In this case, \code{x ~ +macro1(arg1 = z[1:10], arg2 = "hello")} will be passed to \code{fun} +as \code{fun(x ~ macro1(arg1 = z[1:10], arg2 = "hello"))}. It is also +possible in this case to pass a macro without using a \code{~} or +\code{<-}. For example, the line \code{macro1(arg1 = z[1:10], arg2 = +"hello")} will be passed to \code{fun} as \code{fun(macro1(arg1 = +z[1:10], arg2 = "hello"))}. + +If \code{use3pieces = FALSE} and \code{unpackArgs = TRUE}, it +won't make sense to anticipate a declaration using \code{~} or \code{<-}. +Instead, arguments from an arbitrary call will be passed as separate arguments. +For example, the line \code{macro1(arg1 = z[1:10], arg2 = "hello")} will be +passed to \code{fun} as \code{fun(arg1 = z[1:10], arg2 = "hello")}. + +In addition, the final two arguments of \code{fun} must be called \code{modelInfo} +and \code{.env} respectively. + +During macro processing, \code{nimbleModel} passes a named list to the \code{modelInfo} +argument of \code{fun} containing, among other things, elements called +\code{constants} and \code{dimensions}. Macro developers can modify these +two elements (for example, to add a new constant needed for a macro) and +these changes will be reflected in the final model object. Note that currently +it is not possible for a macro to modify the data. Furthermore, if your macro add a new element to the +constants that \code{nimbleModel} then moves to the data, this new data will not be retained +in the final model object and thus will not be usable. + +\code{nimbleModel} passes the R environment from which \code{nimbleModel} was +called to the \code{.env} argument. + +The \code{fun} function must return a named list with two elements: +\code{code}, the replacement code, and \code{modelInfo}, the \code{modelInfo} +list described above. \code{modelInfo} must be in the output even if the macro +does not modify it. + +It is extremely useful to be familiar with processing R code as an +object to write \code{fun} correctly. Functions such as +\code{\link{substitute}} and \code{\link{as.name}} +(e.g. \code{as.name('~')}), \code{\link{quote}}, \code{\link{parse}} +and \code{\link{deparse}} are particularly handy. + +Multiple lines of new code should be contained in \code{ {} }. Extra +curly braces are not a problem. See example 2. + +Macro expansion is done recursively: One macro can return code that +invokes another macro. +} +\examples{ +nimbleOptions(enableMacros = TRUE) +nimbleOptions(enableMacroComments = FALSE) +nimbleOptions(verbose = FALSE) + +## Example 1: Say one is tired of writing "for" loops. +## This macro will generate a "for" loop with dnorm declarations +all_dnorm <- buildMacro( + function(stoch, LHS, RHSvar, start, end, sd = 1, modelInfo, .env) { + newCode <- substitute( + for(i in START:END) { + LHS[i] ~ dnorm(RHSvar[i], SD) + }, + list(START = start, + END = end, + LHS = LHS, + RHSvar = RHSvar, + SD = sd)) + list(code = newCode) + }, + use3pieces = TRUE, + unpackArgs = TRUE +) + +model1 <- nimbleModel( + nimbleCode( + { + ## Create a "for" loop of dnorm declarations by invoking the macro + x ~ all_dnorm(mu, start = 1, end = 10) + } + )) + +## show code from expansion of macro +model1$getCode() +## The result should be: +## { +## for (i in 1:10) { +## x[i] ~ dnorm(mu[i], 1) +## } +## } + +## Example 2: Say one is tired of writing priors. +## This macro will generate a set of priors in one statement +flat_normal_priors <- buildMacro( + function(..., modelInfo, .env) { + allVars <- list(...) + priorDeclarations <- lapply(allVars, + function(x) + substitute(VAR ~ dnorm(0, sd = 1000), + list(VAR = x))) + newCode <- quote({}) + newCode[2:(length(allVars)+1)] <- priorDeclarations + list(code = newCode) + }, + use3pieces = FALSE, + unpackArgs = TRUE +) + +model2 <- nimbleModel( + nimbleCode( + { + flat_normal_priors(mu, beta, gamma) + } + )) + +## show code from expansion of macro +model2$getCode() +## The result should be: +## { +## mu ~ dnorm(0, sd = 1000) +## beta ~ dnorm(0, sd = 1000) +## gamma ~ dnorm(0, sd = 1000) +## } + +## Example 3: Macro that modifies constants +new_constant <- buildMacro( + function(stoch, LHS, RHS, modelInfo, .env) { + # number of elements + n <- as.numeric(length(modelInfo$constants[[deparse(LHS)]])) + code <- substitute({ + for (i in 1:N){ + L[i] ~ dnorm(mu[i], 1) + } + }, list(L = LHS, N = n)) + + # Add a new constant mu + modelInfo$constants$mu <- rnorm(n, 0, 1) + + list(code = code, modelInfo = modelInfo) + }, + use3pieces = TRUE, + unpackArgs = TRUE +) + +const <- list(y = rnorm(10)) +code <- nimbleCode({ + y ~ new_constant() +}) + +mod <- nimbleModel(code = code, constants=const) +mod$getCode() +mod$getConstants() # new constant is here +} diff --git a/packages/nimble/man/configureMCMC.Rd b/packages/nimble/man/configureMCMC.Rd index a39b786a6..8d0e68c18 100644 --- a/packages/nimble/man/configureMCMC.Rd +++ b/packages/nimble/man/configureMCMC.Rd @@ -16,6 +16,10 @@ configureMCMC( onlyRW = FALSE, onlySlice = FALSE, multivariateNodesAsScalars = getNimbleOption("MCMCmultivariateNodesAsScalars"), + samplePredictiveNodes = getNimbleOption("MCMCassignSamplersToPosteriorPredictiveNodes"), + mean = FALSE, + variance = FALSE, + logProb = FALSE, enableWAIC = getNimbleOption("MCMCenableWAIC"), controlWAIC = list(), print = getNimbleOption("verbose"), @@ -57,6 +61,14 @@ The default value is an empty character vector, i.e. no values will be recorded. \item{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.} +\item{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.} + +\item{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}.} + +\item{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}.} + +\item{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).} + \item{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}}.} \item{controlWAIC}{A named list of inputs that control the behavior of the WAIC calculation. See \code{help(waic)}.} diff --git a/packages/nimble/man/derived.Rd b/packages/nimble/man/derived.Rd new file mode 100644 index 000000000..74e104f86 --- /dev/null +++ b/packages/nimble/man/derived.Rd @@ -0,0 +1,42 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/MCMC_derived.R +\name{derived_BASE} +\alias{derived_BASE} +\alias{derived_mean} +\alias{derived_variance} +\alias{derived_logProb} +\alias{derived_predictive} +\alias{derived} +\title{MCMC Derived Quantities} +\usage{ +derived_BASE() + +derived_mean(model, mcmc, interval, control) + +derived_variance(model, mcmc, interval, control) + +derived_logProb(model, mcmc, interval, control) + +derived_predictive(model, mcmc, interval, control) +} +\description{ +Details of the NIMBLE MCMC engine handles derived quantities, which are deterministic functions that can be calaculated and recorded after each MCMC sampling iteration. +} +\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 +} diff --git a/packages/nimble/man/expmA.Rd b/packages/nimble/man/expmA.Rd new file mode 100644 index 000000000..6afb6e469 --- /dev/null +++ b/packages/nimble/man/expmA.Rd @@ -0,0 +1,44 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/miscAlgorithms.r +\name{expmA} +\alias{expmA} +\title{Matrix Exponential} +\usage{ +expmA(A, tol = 1e-08) +} +\arguments{ +\item{A}{Infinitesimal Generator Matrix.} + +\item{tol}{level of accuracy required (default = 1e-8).} +} +\value{ +\code{expmA} gives a matrix that is ans = exp(A). +} +\description{ +Compute the the matrix exponential expm(A) + by scaling and squaring. +} +\details{ +A, the generator matrix must have negative value diagonals. + +This function follows the scaling and squaring algorithm from Sherlock (2021), except that we compute the full +matrix exponential. It differs from the standard Taylor scaling and squaring algorithm reviewed by Ruiz et al (2016) and found in common texts. +If using the matrix exponential to create a transition probability matrix in a HMM context just once, +this function is good. If dimension is large, we recommend avoiding the matrix exponential and using `expmAv` instead. +Note that for computation efficiency matrix uniformization is always done by A* = A + rho I, where rho = max(abs(diag(A))). +Derivatives are currently not supported. +} +\examples{ +A <- rbind(c(-1, 0.25, 0.75), c(0, -2, 2), c(0.25, 0.25, -0.5)) +Lambda <- diag(c(0.25, 0.1, 0)) +expmA((A-Lambda)*2.5) +} +\references{ +Sherlock, C. (2021). Direct statistical inference for finite Markov jump processes via the matrix exponential. Computational Statistics, 36(4), 2863-2887. + +Ruiz, P., Sastre, J., Ibáñez, J., & Defez, E. (2016). High performance computing of the matrix exponential. +Journal of computational and applied mathematics, 291, 370-379. +} +\author{ +Paul van Dam-Bates +} diff --git a/packages/nimble/man/expmAv.Rd b/packages/nimble/man/expmAv.Rd new file mode 100644 index 000000000..91b0178de --- /dev/null +++ b/packages/nimble/man/expmAv.Rd @@ -0,0 +1,52 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/miscAlgorithms.r +\name{expmAv} +\alias{expmAv} +\title{Matrix Exponential times a vector} +\usage{ +expmAv(A, v, tol = 1e-08, rescaleFreq = 10, Nmax = 10000, sparse = FALSE) +} +\arguments{ +\item{A}{Infinitesimal Generator Matrix.} + +\item{v}{vector to multiply by the matrix exponential exp(A) %*% v.} + +\item{tol}{level of accuracy required (default = 1e-8).} + +\item{rescaleFreq}{How frequently should the terms be scaled to avoid underflow/overflow (default = 10).} + +\item{Nmax}{maximum number of iterations to compute (default = 10000).} + +\item{sparse}{(logical) specify if the matrix may be sparse and to do sparse computation (default = FALSE).} +} +\value{ +\code{expmAv} gives a vector that is ans = exp(A) %*% v. +} +\description{ +Compute the combined term expm(A) %*% v + to avoid a full matrix exponentiation. +} +\details{ +For large matrix exponentials it is much more efficient to compute exp(A) %*% v, than to actually compute the entire matrix exponential. A, the generator matrix must have +negative value diagonals. + +This function follows the function `expAv` from the R package RTMB (Kristensen, 2025), and theory outlined in Sherlock (2021). It is developed for working with continuous times +Markov chains. If using the matrix exponential to create a transition probability matrix in a HMM context just once, +this function may be slower than the one time call to compute the full matrix exponentiation. If a full matrix exponentiation is required, refer to `nimbleRcall` and +functions such as `pracma::expm` or `expm::expm` to compute. Choosing sparse = TRUE will then check which matrix A values are non-zero and do sparse linear +algebra on the non-zero terms only. Note that for computation efficiency matrix uniformization is always done by A* = A + rho I, where rho = max(abs(diag(A))), see Algorithm 2#' in Sherlock (2021). +} +\examples{ +A <- rbind(c(-1, 0.25, 0.75), c(0, -2, 2), c(0.25, 0.25, -0.5)) +v <- c(0.35, 0.25, 0.1) +expmAv(A, v) +} +\references{ +Sherlock, C. (2021). Direct statistical inference for finite Markov jump processes via the matrix exponential. Computational Statistics, 36(4), 2863-2887. + +Kristensen K (2025). _RTMB: 'R' Bindings for 'TMB'_. R package version 1.7, commit 6bd7a16403ccb4d3fc13ff7526827540bf27b352, +. +} +\author{ +Paul van Dam-Bates +} diff --git a/packages/nimble/man/getNimbleOption.Rd b/packages/nimble/man/getNimbleOption.Rd index 03f5d2a1d..5fe248183 100644 --- a/packages/nimble/man/getNimbleOption.Rd +++ b/packages/nimble/man/getNimbleOption.Rd @@ -17,7 +17,7 @@ Allow the user to get the value of a global _option_ that affects the way in which NIMBLE operates } \examples{ -getNimbleOption('verifyConjugatePosteriors') +getNimbleOption('MCMCverifyConjugatePosteriors') } \author{ Christopher Paciorek diff --git a/packages/nimble/man/getSamplesDPmeasure.Rd b/packages/nimble/man/getSamplesDPmeasure.Rd index f2bd5f71c..791daa037 100644 --- a/packages/nimble/man/getSamplesDPmeasure.Rd +++ b/packages/nimble/man/getSamplesDPmeasure.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/BNP_samplers.R \name{getSamplesDPmeasure} \alias{getSamplesDPmeasure} -\title{Get posterior samples for a Dirichlet process measure} +\title{Get posterior samples for a Dirichlet process distribution (measure)} \usage{ getSamplesDPmeasure( MCMC, @@ -14,25 +14,25 @@ getSamplesDPmeasure( \arguments{ \item{MCMC}{an MCMC class object, either compiled or uncompiled.} -\item{epsilon}{used for determining the truncation level of the representation of the random measure.} +\item{epsilon}{used for determining the truncation level of the representation of the random distribution (measure).} \item{setSeed}{Logical or numeric argument. If a single numeric value is provided, R's random number seed will be set to this value. In the case of a logical value, if \code{TRUE}, then R's random number seed will be set to \code{1}. Note that specifying the argument \code{setSeed = 0} does not prevent setting the RNG seed, but rather sets the random number generation seed to \code{0}. Default value is \code{FALSE}.} \item{progressBar}{Logical specifying whether to display a progress bar during execution (default = TRUE). The progress bar can be permanently disabled by setting the system option \code{nimbleOptions(MCMCprogressBar = FALSE)}} } \description{ -This function obtains posterior samples from a Dirichlet process distributed random measure of a model specified using the \code{dCRP} distribution. +This function obtains posterior samples from a Dirichlet process distributed random distribution (measure) of a model specified using the \code{dCRP} distribution. } \details{ -This function provides samples from a random measure having a Dirichlet process prior. Realizations are almost surely discrete and represented by a (finite) stick-breaking representation (Sethuraman, 1994), whose atoms (or point masses) are independent and identically distributed. This sampler can only be used with models containing a \code{dCRP} distribution. +This function provides samples from a random distribution (measure) having a Dirichlet process prior. Realizations are almost surely discrete and represented by a (finite) stick-breaking representation (Sethuraman, 1994), whose atoms (or point masses) are independent and identically distributed. This sampler can only be used with models containing a \code{dCRP} distribution. -The \code{MCMC} argument is an object of class MCMC provided by \code{buildMCMC}, or its compiled version. The MCMC should already have been run, as \code{getSamplesDPmeasure} uses the posterior samples to generate samples of the random measure. Note that the monitors associated with that MCMC must include the cluster membership variable (which has the \code{dCRP} distribution), the cluster parameter variables, all variables directly determining the \code{dCRP} concentration parameter, and any stochastic parent variables of the cluster parameter variables. See \code{help(configureMCMC)} or \code{help(addMonitors)} for information on specifying monitors for an MCMC. +The \code{MCMC} argument is an object of class MCMC provided by \code{buildMCMC}, or its compiled version. The MCMC should already have been run, as \code{getSamplesDPmeasure} uses the posterior samples to generate samples of the random distribution (measure). Note that the monitors associated with that MCMC must include the cluster membership variable (which has the \code{dCRP} distribution), the cluster parameter variables, all variables directly determining the \code{dCRP} concentration parameter, and any stochastic parent variables of the cluster parameter variables. See \code{help(configureMCMC)} or \code{help(addMonitors)} for information on specifying monitors for an MCMC. The \code{epsilon} argument is optional and used to determine the truncation level of the random measure. \code{epsilon} is the tail probability of the random measure, which together with posterior samples of the concentration parameter, determines the truncation level. The default value is 1e-4. -The output is a list of matrices. Each matrix represents a sample from the random measure. In order to reduce the output's dimensionality, the weights of identical atoms are added up. The stick-breaking weights are named \code{weights} and the atoms are named based on the cluster variables in the model. +The output is a list of matrices. Each matrix represents a sample from the random distribution (measure). In order to reduce the output's dimensionality, the weights of identical atoms are added up. The stick-breaking weights are named \code{weights} and the atoms are named based on the cluster variables in the model. -For more details about sampling the random measure and determining its truncation level, see Section 3 in Gelfand, A.E. and Kottas, A. 2002. +For more details about sampling the random distribution (measure) and determining its truncation level, see Section 3 in Gelfand, A.E. and Kottas, A. 2002. } \examples{ \dontrun{ diff --git a/packages/nimble/man/nimDerivs.Rd b/packages/nimble/man/nimDerivs.Rd index 2e1d1ebaf..e40f3354b 100644 --- a/packages/nimble/man/nimDerivs.Rd +++ b/packages/nimble/man/nimDerivs.Rd @@ -23,10 +23,13 @@ all arguments to \code{nimFxn}.} corresponding to whether the function value, Jacobian, and Hessian should be returned respectively. Defaults to \code{c(0, 1, 2)}.} -\item{model}{(optional) for derivatives of a nimbleFunction that involves model. -calculations, the uncompiled model that is used. This is needed in order +\item{model}{(optional) the uncompiled model that is used, if taking derivatives +of a nimbleFunction that involves model calculations. This is needed in order to be able to correctly restore values into the model when \code{order} does not -include 0 (or in all cases when double-taping).} +include 0 (or in all cases when double-taping). IMPORTANT: if \code{model} +is included, one should also include the arguments \code{updateNodes} and +\code{constantNodes} using the output obtained from running +\code{makeModelDerivsInfo}.} \item{...}{additional arguments intended for internal use only.} } @@ -43,6 +46,12 @@ Derivatives for uncompiled nimbleFunctions are calculated using the \code{numDeriv} package. If this package is not installed, an error will be issued. Derivatives for matrix valued arguments will be returned in column-major order. + +As discussed above with the \code{model} argument, if taking derivatives +of a nimbleFunction that involves model calculations (rather than directly +taking derivatives of `calculate`), care needs to be taken to provide +\code{model}, \code{updateNodes}, and \code{calcNodes} arguments. See +Section 16.7.2 of the User Manual for more details. } \examples{ diff --git a/packages/nimble/man/nimbleMCMC.Rd b/packages/nimble/man/nimbleMCMC.Rd index 5a7c77a68..5e948db43 100644 --- a/packages/nimble/man/nimbleMCMC.Rd +++ b/packages/nimble/man/nimbleMCMC.Rd @@ -22,7 +22,11 @@ nimbleMCMC( samples = TRUE, samplesAsCodaMCMC = FALSE, summary = FALSE, - WAIC = FALSE + mean = FALSE, + variance = FALSE, + logProb = FALSE, + WAIC = FALSE, + userEnv = parent.frame() ) } \arguments{ @@ -58,8 +62,13 @@ nimbleMCMC( \item{samplesAsCodaMCMC}{Logical argument. If \code{TRUE}, then a \code{coda} \code{mcmc} object is returned instead of an R matrix of samples, or when \code{nchains > 1} a \code{coda} \code{mcmc.list} object is returned containing \code{nchains} \code{mcmc} objects. This argument is only used when \code{samples} is \code{TRUE}. Default value is \code{FALSE}. See details.} -\item{summary}{Logical argument. When \code{TRUE}, summary statistics for the posterior samples of each parameter are also returned, for each MCMC chain. This may be returned in addition to the posterior samples themselves. Default value is \code{FALSE}. See details. -z} +\item{summary}{Logical argument. When \code{TRUE}, summary statistics for the posterior samples of each parameter are also returned, for each MCMC chain. This may be returned in addition to the posterior samples themselves. Default value is \code{FALSE}. See details.} + +\item{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}.} + +\item{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}.} + +\item{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).} \item{WAIC}{Logical argument. When \code{TRUE}, the WAIC (Watanabe, 2010) of the model is calculated and returned. If multiple chains are run, then a single WAIC value is calculated using the posterior samples from all chains. Default value is \code{FALSE}. Note that the version of WAIC used is the default WAIC conditional on random effects/latent states and without any grouping of data nodes. See \code{help(waic)} for more details. 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}.} } diff --git a/packages/nimble/man/nimbleOptions.Rd b/packages/nimble/man/nimbleOptions.Rd index f7709ff23..cd1078891 100644 --- a/packages/nimble/man/nimbleOptions.Rd +++ b/packages/nimble/man/nimbleOptions.Rd @@ -27,7 +27,7 @@ with no arguments to see a list of available opions. } \examples{ # Set one option: -nimbleOptions(verifyConjugatePosteriors = FALSE) +nimbleOptions(MCMCverifyConjugatePosteriors = FALSE) # Compactly print all options: str(nimbleOptions(), max.level = 1) diff --git a/packages/nimble/man/runMCMC.Rd b/packages/nimble/man/runMCMC.Rd index f04cc8493..fce9ff8f7 100644 --- a/packages/nimble/man/runMCMC.Rd +++ b/packages/nimble/man/runMCMC.Rd @@ -17,6 +17,8 @@ runMCMC( samples = TRUE, samplesAsCodaMCMC = FALSE, summary = FALSE, + derivedQuantities = (mcmc$getNumDerived() > 0) && + getNimbleOption("MCMCreturnDerivedQuantities"), WAIC = FALSE, perChainWAIC = FALSE ) @@ -46,6 +48,8 @@ runMCMC( \item{summary}{Logical argument. When \code{TRUE}, summary statistics for the posterior samples of each parameter are also returned, for each MCMC chain. This may be returned in addition to the posterior samples themselves. Default value is \code{FALSE}. See details.} +\item{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.} + \item{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)}.} \item{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}.} diff --git a/packages/nimble/man/samplers.Rd b/packages/nimble/man/samplers.Rd index f2ae019cf..bcf4e4941 100644 --- a/packages/nimble/man/samplers.Rd +++ b/packages/nimble/man/samplers.Rd @@ -125,7 +125,7 @@ sampler_slice_CRP_base_param(model, mvSaved, target, control) \item{control}{named list that controls the precise behavior of the sampler, with elements specific to \code{samplertype}. The default values for control list are specified in the setup code of each sampling algorithm. Descriptions of each sampling algorithm, and the possible customizations for each sampler (using the \code{control} argument) appear below.} } \description{ -Details of the MCMC sampling algorithms provided with the NIMBLE MCMC engine; HMC samplers are in the \code{nimbleHMC} package and particle filter samplers are in the \code{nimbleSMC} package. +Details of the MCMC sampling algorithms provided with the NIMBLE MCMC engine; HMC samplers are in the \code{nimbleHMC} package and particle filter samplers are in the \code{nimbleSMC} package. Additional details, including some recommendations for samplers that may perform better than the samplers that NIMBLE assigns by default are provided in Section 7.11 of the User Manual. } \section{Hamiltonian Monte Carlo samplers}{ diff --git a/packages/nimble/tests/testthat/test-miscAlgorithms.R b/packages/nimble/tests/testthat/test-miscAlgorithms.R new file mode 100644 index 000000000..3cbf83d12 --- /dev/null +++ b/packages/nimble/tests/testthat/test-miscAlgorithms.R @@ -0,0 +1,143 @@ +source(system.file(file.path('tests', 'testthat', 'test_utils.R'), package = 'nimble')) + +# Tests of Misc Algorithms which currently includes expAv and expm +test_that("expAv works.", { + + ## Note that pracma was not accurate enough even for a small example + ## expm::expm(A) %*% v + A <- matrix(c(-0.120, 100, 0.9, -150), 2, 2) + v <- c(120, 0) + eAv <- c(192.793313160, 128.120518022) + + ## Wrapper as derivs = TRUE with no setup. + expAv_wrap <- nimbleFunction( + run = function(A = double(2), v = double(1), + tol = double(0, default = 1e-8), + rescaleFreq = double(0, default = 10), + Nmax = integer(0, default = 10000), + sparse = logical(0, default = FALSE)){ + ans <- expAv(A, v, tol, rescaleFreq, Nmax, sparse) + returnType(double(1)) + return(ans) + } + ) + expAvc <- compileNimble(expAv_wrap) + + eAvNimR <- expAv(A, v) + eAvNim <- expAvc(A, v) + expect_equal(eAv, eAvNim, tol = 1e-8) + expect_equal(eAv, eAvNimR, tol = 1e-8) + + ## Increase tolerance. + eAvNim2 <- expAvc(A, v, tol = 1e-16) + expect_equal(eAv, eAvNim2, tol = 1e-10) + + ## Compare sparse version with regular version: + set.seed(1) + A <- as.matrix(Matrix::rsparsematrix(10, 10, .5)) + A <- A - diag(diag(A)) + v <- rnorm(10) + d <- -rowSums(A) - 8 + diag(A) <- d + eAvdense <- expAvc(A, v, sparse = FALSE) + eAvsparse <- expAvc(A, v, sparse = TRUE) + expect_equal(eAvdense, eAvsparse, tol = 1e-15) + + ## Increase tolerance: + eAvsparse2 <- expAvc(as.matrix(A), v, sparse = TRUE, tol = 1e-16) + expect_equal(eAvsparse2, eAvsparse, tol = 1e-8) ## Comparing 1e-8 and 1e-16 + eAvsparse3 <- expAvc(as.matrix(A), v, sparse = TRUE, tol = 1e-5) + expect_equal(eAvsparse3, eAvsparse2, tol = 1e-5) ## Comparing 1e-8 and 1e-16 + + ## Check rescaling: + eAvsparse_scale1 <- expAvc(as.matrix(A), v, sparse = TRUE, tol = 1e-8, rescaleFreq = 1) + eAvsparse_scale10 <- expAvc(as.matrix(A), v, sparse = TRUE, tol = 1e-8, rescaleFreq = 10) + eAvsparse_scale50 <- expAvc(as.matrix(A), v, sparse = TRUE, tol = 1e-8, rescaleFreq = 50) + expect_equal(eAvsparse_scale1, eAvsparse, tol = 1e-15) ## Approximation is best with frequent rescaling. + expect_equal(eAvsparse_scale10, eAvsparse, tol = 1e-12) + expect_equal(eAvsparse_scale50, eAvsparse, tol = 1e-10) + + ## Test larger matrix that has positive diagonals, does not require negative. + set.seed(2) + A <- as.matrix(Matrix::rsparsematrix(50, 50, 0.25)) + v <- rnorm(50) + eAvdense <- expAvc(A, v, sparse = FALSE) + eAvsparse <- expAvc(A, v, sparse = TRUE) + expect_equal(eAvdense, eAvsparse, tol = 1e-15) + + # expAv_rtmb <- RTMB::expAv(A, v) + ## Indices 1,5, 10 Check sprintf("%.16f",expAv_rtmb[c(1,5,10)]) + expAv_rtmb <- c(-10.4849004612874221, 5.7463833745583406, -4.3508466614873917) + expect_equal(eAvsparse[c(1,5,10)], expAv_rtmb, tol = 1e-8) + + ## Larger numerical test with bigger time period, leads to overflow for both expm package and ours, with same result. + # expAv_rtmb <- RTMB::expAv(A*10, v) + # Indices 1,5, 10 Check sprintf("%.16f",expAv_rtmb[c(1,5,10)]) + vals <- c(8591151029674.1933593750000000, -2283297500254.1206054687500000, -16239759412070.7285156250000000) + eAv <- expAvc(A*10, v, tol = 1e-8) + expect_equal(eAv[c(1,5,10)], vals, tol = 1e-8) +}) + +# Test full matrix exponential using scaling and squaring. +test_that("Matrix Exponential Works", { + + ## expm::expm(A) %*% v + A <- matrix(c(-0.120, 100, 0.9, -150), 2, 2) + v <- c(120, 0) + eAv <- c(192.793313160, 128.120518022) + + ## Wrapper as derivs = TRUE with no setup. + expm_wrap <- nimbleFunction( + run = function(A = double(2), tol = double(0, default = 1e-8)){ + returnType(double(2)) + ans <- expm(A, tol) + return(ans) + } + ) + expmc <- compileNimble(expm_wrap) + + eAR <- expm(A) + eAC <- expmc(A) + eAvR <- (eAR %*% v)[,1] + eAvC <- (eAC %*% v)[,1] + expect_equal(eAv, eAvR, tol = 1e-6) + expect_equal(eAv, eAvC, tol = 1e-6) + + ## Increase tolerance. + eAvNim2 <- (expmc(A, tol = 1e-16) %*% v)[,1] + expect_equal(eAv, eAvNim2, tol = 1e-10) + + ## Bigger Matrix: + set.seed(101) + n <- 10 + qd <- rgamma(n, 1, 1) + A <- matrix(0, n, n) + A[1,2] <- qd[1] + for( i in 2:(n-1) ){ + A[i, i+1] <- qd[i]/2 + A[i, i-1] <- qd[i]/2 + } + A[n, n-1] <- qd[n] + diag(A) <- -qd + v <- rgamma(n, 5, 1) + eAv <- expAv(A, v, tol = 1e-16) + eAv2 <- (expmc(A, tol = 1e-16) %*% v)[,1] + expect_equal(eAv, eAv2, tol = 1e-14) + + ## Now a much larger matrix and compare with expm::expm... + set.seed(2) + A <- as.matrix(Matrix::rsparsematrix(50, 50, 0.25)) + + # test <- expm::expm(A) + # sprintf("%.16f",test[cbind(c(1,10,20,32),c(1,15,25,32))]) + vals <- c(1.4785161123142350, 2.9154566237879371, -0.9980758066172498, 1.1608012543866268) + eAv <- expm(A, tol = 1e-14) + expect_equal(eAv[cbind(c(1,10,20,32),c(1,15,25,32))], vals, tol = 1e-14) + + ## Larger numerical test with bigger time period, leads to overflow for both expm package and ours, with same result. + # test <- expm::expm(A*10) + # sprintf("%.16f",test[cbind(c(1,10,20,32),c(1,15,25,32))]) + vals <- c(2247213364689.7807617187500000, 1321534261508.6760253906250000, 81388538896.9135742187500000, -666846013929.7093505859375000) + eAv <- expm(A*10, tol = 1e-14) + expect_equal(eAv[cbind(c(1,10,20,32),c(1,15,25,32))], vals, tol = 1e-12) +})