Package 'nimble'

Description

This function is used in a method of a nimbleFunction that has derivatives enabled. It returns its value but breaks tracking of derivatives.

Usage

ADbreak(x)

Arguments

Details

This funcion only works with scalars.

Data type for the return value of nimDerivs

Description

nimbleList definition for the type of nimbleList returned by nimDerivs.

Usage

ADNimbleList

Format

An object of class list of length 1.

Fields

value

The value of the function evaluated at the given input arguments.

jacobian

The Jacobian of the function evaluated at the given input arguments.

hessian

The Hessian of the function evaluated at the given input arguments.

See Also

nimDerivs

Description

create an ADproxyModelClass object. For internal use.

Arguments

Details

This is a proxy model for model_AD. The class needs just enough pieces to be used like a model for purposes of nodeFunction compilation. The model will contain an ADproxyModel and then the nodeFunction setup code will extract it. The model interface will population the proxy model's CobjectInterface

Author(s)

NIMBLE development team

Description

NIMBLE language functions that can be used in either compiled or uncompiled nimbleFunctions to detect if there are any NA or NaN values in a vector.

Usage

any_na(x)

any_nan(x)

Arguments

Author(s)

NIMBLE Development Team

Description

This will convert alternate representations of CAR process structure into (adj, weights, num) form required by dcar_normal.

Usage

as.carAdjacency(...)

Arguments

Details

Two alternate representations are handled:

A single matrix argument will be interpreted as a matrix of symmetric unnormalized weights.

Two lists will be interpreted as (the first) a list of numeric vectors specifying the adjacency (neighboring) indices of each CAR process component, and (the second) a list of numeric vectors giving the unnormalized weights for each of these neighboring relationships.

Author(s)

Daniel Turek

See Also

CAR-Normal

Description

Convert weights vector to C and M parameters of dcar_proper distribution

Usage

as.carCM(adj, weights, num)

Arguments

Details

Given a symmetric matrix of unnormalized weights, this function will calculate corresponding values for the C and M arguments suitable for use in the dcar_proper distribution. This function can be used to transition between usage of dcar_normal and dcar_proper, since dcar_normal uses the adj, weights, and num arguments, while dcar_proper requires adj, num, and also the C and M as returned by this function.

Here, C is a sparse vector representation of the row-normalized adjacency matrix, and M is a vector containing the conditional variance for each region. The resulting values of C and M are guaranteed to satisfy the symmetry constraint imposed on $C$ and $M$, that $M^{-1} C$ is symmetric, where $M$ is a diagonal matrix and $C$ is the row-normalized adjacency matrix.

Value

A named list with elements C and M. These may be used as the C and M arguments to the dcar_proper distribution.

Author(s)

Daniel Turek

See Also

CAR-Normal, CAR-Proper

Description

Turns a numeric vector into a matrix that has 1 row or 1 column. Part of NIMBLE language.

Usage

asRow(x)

asCol(x)

Arguments

Details

In the NIMBLE language, some automatic determination of how to turn vectors into single-row or single-column matrices is done. For example, in A %*% x, where A is a matrix and x a vector, x will be turned into a single-column matrix unless it is known at compile time that A is a single column, in which case x will be turned into a single-row matrix. However, if it is desired that x be turned into a single row but A cannot be determined at compile time to be a single column, then one can use A %*% asRow(x) to force this conversion.

Author(s)

Perry de Valpine

Description

The automated parameter blocking algorithm is no longer actively maintained. In some cases, it may not operate correctly with more recent system features and/or distributions.

Usage

autoBlock(
  Rmodel,
  autoIt = 20000,
  run = list("all", "default"),
  setSeed = TRUE,
  verbose = FALSE,
  makePlots = FALSE,
  round = TRUE
)

Arguments

Details

Runs NIMBLE's automated blocking procedure for a given model object, to dynamically determine a blocking scheme of the continuous-valued model nodes. This blocking scheme is designed to produce efficient MCMC sampling (defined as number of effective samples generated per second of algorithm runtime). See Turek, et al (2015) for details of this algorithm. This also (optionally) compares this blocked MCMC against several static MCMC algorithms, including all univariate sampling, blocking of all continuous-valued nodes, NIMBLE's default MCMC configuration, and custom-specified blockings of parameters.

This method allows for fine-tuned usage of the automated blocking procedure. However, the main entry point to the automatic blocking procedure is intended to be through either buildMCMC(..., autoBlock = TRUE), or configureMCMC(..., autoBlock = TRUE).

Value

Returns a named list containing elements:

• summary: A data frame containing a numerical summary of the performance of all MCMC algorithms (including that from automated blocking)

• autoGroups: A list specifying the parameter blockings converged on by the automated blocking procedure

• conf: A NIMBLE MCMC configuration object corresponding to the results of the automated blocking procedure

Author(s)

Daniel Turek

References

Turek, D., de Valpine, P., Paciorek, C., and Anderson-Bergman, C. (2015). Automated Parameter Blocking for Efficient Markov-Chain Monte Carlo Sampling. <arXiv:1503.05621>.

See Also

configureMCMC buildMCMC

Description

BUGSdeclClass contains the information extracted from one BUGS declaration

Description

Create quadrature grid for use in AGHQuad methods in Nimble.

Arguments

Details

This function is used by used by buildOneAGHQuad1D and buildOneAGHQuad create the quadrature grid using adaptive Gauss-Hermite quadrature. Handles single or multiple dimension grids and computes both grid locations and weights. Additionally, acts as a cache system to do transformations, and return marginalized log density.

Any of the input node vectors, when provided, will be processed using nodes <- model$expandNodeNames(nodes), where nodes may be paramNodes, randomEffectsNodes, and so on. This step allows any of the inputs to include node-name-like syntax that might contain multiple nodes. For example, paramNodes = 'beta[1:10]' can be provided if there are actually 10 scalar parameters, 'beta[1]' through 'beta[10]'. The actual node names in the model will be determined by the exapndNodeNames step.

Available methods include

• buildAGHQ. Builds a adaptive Gauss-Hermite quadrature grid in d dimensions. Calls buildAGHQOne to build the one dimensional grid and then expands in each dimension. Some numerical issues occur in Eigen decomposition making the grid weights only accurate up to 35 quadrature nodes.

• Options to get internally cached values are getGridSize, getModeIndex for when there are an odd number of quadrature nodes, getLogDensity for the cached values, getAllNodes for the quadrature grids, getNodes for getting a single indexed nodes, getAllNodesTransformed for nodes transformed to the parameter scale, getNodesTransformed for a single transformed node, getAllWeights to get all quadrature weights, getWeights single indexed weight.

• transformGrid(cholNegHess, inner_mode, method) transforms the grid using either cholesky trasnformations, as default, or spectral that makes use of the Eigen decomposition. For a single dimension transformGrid1D is used.

• As the log density is evaluated externally, it is saved via saveLogDens, which then is summed via quadSum.

• buildGrid builds the grid the initial time and is only run once in code. After, the user must choose to setGridSize to update the grid size.

• check. If TRUE (default), a warning is issued if paramNodes, randomEffectsNodes and/or calcNodes are provided but seek to have missing elements or unnecessary elements based on some default inspection of the model. If unnecessary warnings are emitted, simply set check=FALSE.

• innerOptimControl. A list of control parameters for the inner optimization of Laplace approximation using optim. See 'Details' of optim for further information.

• innerOptimMethod. Optimization method to be used in optim for the inner optimization. See 'Details' of optim. Currently optim in NIMBLE supports: "Nelder-Mead", "BFGS", "CG", and "L-BFGS-B". By default, method "CG" is used when marginalizing over a single (scalar) random effect, and "BFGS" is used for multiple random effects being jointly marginalized over.

• innerOptimStart. Choice of starting values for the inner optimization. This could be "last", "last.best", or a vector of user provided values. "last" means the most recent random effects values left in the model will be used. When finding the MLE, the most recent values will be the result of the most recent inner optimization for Laplace. "last.best" means the random effects values corresponding to the largest Laplace likelihood (from any call to the calcLaplace or calcLogLik method, including during an MLE search) will be used (even if it was not the most recent Laplace likelihood). By default, the initial random effects values will be used for inner optimization.

• outOptimControl. A list of control parameters for maximizing the Laplace log-likelihood using optim. See 'Details' of optim for further information.

References

Golub, G. H. and Welsch, J. H. (1969). Calculation of Gauss Quadrature Rules. Mathematics of Computation 23 (106): 221-230.

Liu, Q. and Pierce, D. A. (1994). A Note on Gauss-Hermite Quadrature. Biometrika, 81(3) 624-629.

Jackel, P. (2005). A note on multivariate Gauss-Hermite quadrature. London: ABN-Amro. Re.

Description

This function has been moved to the 'nimbleSMC' package.

Usage

buildAuxiliaryFilter(...)

Arguments

Description

This function has been moved to the 'nimbleSMC' package.

Usage

buildBootstrapFilter(...)

Arguments

Description

This function has been moved to the 'nimbleSMC' package.

Usage

buildEnsembleKF(...)

Arguments

Description

This function has been moved to the 'nimbleSMC' package.

Usage

buildIteratedFilter2(...)

Arguments

Description

Build a Laplace or AGHQ approximation algorithm for a given NIMBLE model.

Usage

buildLaplace(
  model,
  paramNodes,
  randomEffectsNodes,
  calcNodes,
  calcNodesOther,
  control = list()
)

buildAGHQ(
  model,
  nQuad = 1,
  paramNodes,
  randomEffectsNodes,
  calcNodes,
  calcNodesOther,
  control = list()
)

Arguments

buildLaplace

buildLaplace creates an object that can run Laplace approximation and for a given model or part of a model. buildAGHQ creates an object that can run adaptive Gauss-Hermite quadrature (AGHQ, sometimes called "adaptive Gaussian quadrature") for a given model or part of a model. Laplace approximation is AGHQ with one quadrature point, hence 'buildLaplace' simply calls 'buildAGHQ' with 'nQuad=1'. These methods approximate the integration over continuous random effects in a hierarchical model to calculate the (marginal) likelihood.

buildAGHQ and buildLaplace will by default (unless changed manually via 'control$split') determine from the model which random effects can be integrated over (marginalized) independently. For example, in a GLMM with a grouping factor and an independent random effect intercept for each group, the random effects can be marginalized as a set of univariate approximations rather than one multivariate approximation. On the other hand, correlated or nested random effects would require multivariate marginalization.

Maximum likelihood estimation is available for Laplace approximation ('nQuad=1') with univariate or multivariate integrations. With 'nQuad > 1', maximum likelihood estimation is available only if all integrations are univariate (e.g., a set of univariate random effects). If there are multivariate integrations, these can be calculated at chosen input parameters but not maximized over parameters. For example, one can find the MLE based on Laplace approximation and then increase 'nQuad' (using the 'updateSettings' method below) to check on accuracy of the marginal log likelihood at the MLE.

Beware that quadrature will use 'nQuad^k' quadrature points, where 'k' is the dimension of each integration. Therefore quadrature for 'k' greater that 2 or 3 can be slow. As just noted, 'buildAGHQ' will determine independent dimensions of quadrature, so it is fine to have a set of univariate random effects, as these will each have k=1. Multivariate quadrature (k>1) is only necessary for nested, correlated, or otherwise dependent random effects.

The recommended way to find the maximum likelihood estimate and associated outputs is by calling runLaplace or runAGHQ. The input should be the compiled Laplace or AGHQ algorithm object. This would be produced by running compileNimble with input that is the result of buildLaplace or buildAGHQ.

For more granular control, see below for methods findMLE and summary. See function summaryLaplace for an easier way to call the summary method and obtain results that include node names. These steps are all done within runLaplace and runAGHQ.

The NIMBLE User Manual at r-nimble.org also contains an example of Laplace approximation.

How input nodes are processed

buildLaplace and buildAGHQ make good tries at deciding what to do with the input model and any (optional) of the node arguments. However, random effects (over which approximate integration will be done) can be written in models in multiple equivalent ways, and customized use cases may call for integrating over chosen parts of a model. Hence, one can take full charge of how different parts of the model will be used.

Any of the input node vectors, when provided, will be processed using nodes <- model$expandNodeNames(nodes), where nodes may be paramNodes, randomEffectsNodes, and so on. This step allows any of the inputs to include node-name-like syntax that might contain multiple nodes. For example, paramNodes = 'beta[1:10]' can be provided if there are actually 10 scalar parameters, 'beta[1]' through 'beta[10]'. The actual node names in the model will be determined by the exapndNodeNames step.

In many (but not all) cases, one only needs to provide a NIMBLE model object and then the function will construct reasonable defaults necessary for Laplace approximation to marginalize over all continuous latent states (aka random effects) in a model. The default values for the four groups of nodes are obtained by calling setupMargNodes, whose arguments match those here (except for a few arguments which are taken from control list elements here).

setupMargNodes tries to give sensible defaults from any combination of paramNodes, randomEffectsNodes, calcNodes, and calcNodesOther that are provided. For example, if you provide only randomEffectsNodes (perhaps you want to marginalize over only some of the random effects in your model), setupMargNodes will try to determine appropriate choices for the others.

setupMargNodes also determines which integration dimensions are conditionally independent, i.e., which can be done separately from each other. For example, when possible, 10 univariate random effects will be split into 10 univariate integration problems rather than one 10-dimensional integration problem.

The defaults make general assumptions such as that randomEffectsNodes have paramNodes as parents. However, The steps for determining defaults are not simple, and it is possible that they will be refined in the future. It is also possible that they simply don't give what you want for a particular model. One example where they will not give desired results can occur when random effects have no prior parameters, such as 'N(0,1)' nodes that will be multiplied by a scale factor (e.g. sigma) and added to other explanatory terms in a model. Such nodes look like top-level parameters in terms of model structure, so you must provide a randomEffectsNodes argument to indicate which they are.

It can be helpful to call setupMargNodes directly to see exactly how nodes will be arranged for Laplace approximation. For example, you may want to verify the choice of randomEffectsNodes or get the order of parameters it has established to use for making sense of the MLE and results from the summary method. One can also call setupMargNodes, customize the returned list, and then provide that to buildLaplace as paramNodes. In that case, setupMargNodes will not be called (again) by buildLaplace.

If setupMargNodes is emitting an unnecessary warning, simply use control=list(check=FALSE).

Managing parameter transformations that may be used internally

If any paramNodes (parameters) or randomEffectsNodes (random effects / latent states) have constraints on the range of valid values (because of the distribution they follow), they will be used on a transformed scale determined by parameterTransform. This means the Laplace approximation itself will be done on the transformed scale for random effects and finding the MLE will be done on the transformed scale for parameters. For parameters, prior distributions are not included in calculations, but they are used to determine valid parameter ranges and hence to set up any transformations. For example, if sigma is a standard deviation, you can declare it with a prior such as sigma ~ dhalfflat() to indicate that it must be greater than 0.

For default determination of when transformations are needed, all parameters must have a prior distribution simply to indicate the range of valid values. For a param p that has no constraint, a simple choice is p ~ dflat().

Understanding inner and outer optimizations

Note that there are two numerical optimizations when finding maximum likelihood estimates with a Laplace or (1D) AGHQ algorithm: (1) maximizing the joint log-likelihood of random effects and data given a parameter value to construct the approximation to the marginal log-likelihood at the given parameter value; (2) maximizing the approximation to the marginal log-likelihood over the parameters. In what follows, the prefix 'inner' refers to optimization (1) and 'outer' refers to optimization (2). Currently both optimizations default to using method "BFGS". However, one can use other optimizers or simply run optimization (2) manually from R; see the example below. In some problems, choice of inner and/or outer optimizer can make a big difference for obtaining accurate results, especially for standard errors. Hence it is worth experimenting if one is in doubt.

control list arguments

The control list allows additional settings to be made using named elements of the list. Most (or all) of these can be updated later using the 'updateSettings' method. Supported elements include:

• split. If TRUE (default), randomEffectsNodes will be split into conditionally independent sets if possible. This facilitates more efficient Laplace or AGHQ approximation because each conditionally independent set can be marginalized independently. If FALSE, randomEffectsNodes will be handled as one multivariate block, with one multivariate approximation. If split is a numeric vector, randomEffectsNodes will be split by calling split(randomEffectsNodes, control$split). The last option allows arbitrary control over how randomEffectsNodes are blocked.

• check. If TRUE (default), a warning is issued if paramNodes, randomEffectsNodes and/or calcNodes are provided but seem to have missing or unnecessary elements based on some default inspections of the model. If unnecessary warnings are emitted, simply set check=FALSE.

• innerOptimControl. A list (either an R list or a 'optimControlNimbleList') of control parameters for the inner optimization of Laplace approximation using nimOptim. See 'Details' of nimOptim for further information. Default is 'nimOptimDefaultControl()'.

• innerOptimMethod. Optimization method to be used in nimOptim for the inner optimization. See 'Details' of nimOptim. Currently nimOptim in NIMBLE supports: "Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "nlminb", and user-provided optimizers. By default, method "BFGS" is used for both univariate and multivariate cases. For "nlminb" or user-provided optimizers, only a subset of elements of the innerOptimControlList are supported. (Note that control over the outer optimization method is available as an argument to 'findMLE'). Choice of optimizers can be important and so can be worth exploring.

• innerOptimStart. Method for determining starting values for the inner optimization. Options are:

  • "zero" (default): use all zeros;

  • "last": use the result of the last inner optimization;

  • "last.best": use the result of the best inner optimization so far for each conditionally independent part of the approximation;

  • "constant": always use the same values, determined by innerOptimStartValues;

  • "random": randomly draw new starting values from the model (i.e., from the prior);

  • "model": use values for random effects stored in the model, which are determined from the first call.

  Note that "model" and "zero" are shorthand for "constant" with particular choices of innerOptimStartValues. Note that "last" and "last.best" require a choice for the very first values, which will come from innerOptimStartValues. The default is innerOptimStart="zero" and may change in the future.

• innerOptimStartValues. Values for some of innerOptimStart approaches. If a scalar is provided, that value is used for all elements of random effects for each conditionally independent set. If a vector is provided, it must be the length of *all* random effects. If these are named (by node names), the names will be used to split them correctly among each conditionally independent set of random effects. If they are not named, it is not always obvious what the order should be because it may depend on the conditionally independent sets of random effects. It should match the order of names returned as part of 'summaryLaplace'.

• innerOptimWarning. If FALSE (default), do not emit warnings from the inner optimization. Optimization methods may sometimes emit a warning such as for bad parameter values encountered during the optimization search. Often, a method can recover and still find the optimum. In the approximations here, sometimes the inner optimization search can fail entirely, yet the outer optimization see this as one failed parameter value and can recover. Hence, it is often desirable to silence warnings from the inner optimizer, and this is done by default. Set innerOptimWarning=TRUE to see all warnings.

• useInnerCache. If TRUE (default), use caching system for efficiency of inner optimizations. The caching system records one set of previous parameters and uses the corresponding results if those parameters are used again (e.g., in a gradient call). This should generally not be modified.

• outerOptimControl. A list of control parameters for maximizing the Laplace log-likelihood using nimOptim. See 'Details' of nimOptim for further information.

• computeMethod. There are three approaches available for internal details of how the approximations, and specifically derivatives involved in their calculation, are handled. These are labeled simply 1, 2, and 3, and the default is 2. The relative performance of the methods will depend on the specific model. Users wanting to explore efficiency can try switching from method 2 (default) to methods 1 or 3 and comparing performance. The first Laplace approximation with each method will be (much) slower than subsequent Laplace approximations. Further details are not provided at this time.

• gridType (relevant only nQuad>1). For multivariate AGHQ, a grid must be constructed based on the Hessian at the inner mode. Options include "cholesky" (default) and "spectral" (i.e., eigenvectors and eigenvalues) for the corresponding matrix decompositions on which the grid can be based.

# end itemize

Available methods

The object returned by buildLaplace is a nimbleFunction object with numerous methods (functions). Here these are described in three tiers of user relevance.

Most useful methods

The most relevant methods to a user are:

• calcLogLik(p, trans=FALSE). Calculate the approximation to the marginal log-likelihood function at parameter value p, which (if trans is FALSE) should match the order of paramNodes. For any non-scalar nodes in paramNodes, the order within the node is column-major. The order of names can be obtained from method getNodeNamesVec(TRUE). Return value is the scalar (approximate, marginal) log likelihood.

  If trans is TRUE, then p is the vector of parameters on the transformed scale, if any, described above. In this case, the parameters on the original scale (as the model was written) will be determined by calling the method pInverseTransform(p). Note that the length of the parameter vector on the transformed scale might not be the same as on the original scale (because some constraints of non-scalar parameters result in fewer free transformed parameters than original parameters).

• calcLaplace(p, trans). This is the same as calcLogLik but requires that the approximation be Laplace (i.e nQuad is 1), and results in an error otherwise.

• findMLE(pStart, method, hessian). Find the maximum likelihood estimates of parameters using the approximated marginal likelihood. This can be used if nQuad is 1 (Laplace case) or if nQuad>1 and all marginalizations involve only univariate random effects. Arguments include pStart: initial parameter values (defaults to parameter values currently in the model); method: (outer) optimization method to use in nimOptim (defaults to "BFGS", although some problems may benefit from other choices); and hessian: whether to calculate and return the Hessian matrix (defaults to TRUE, which is required for subsequent use of 'summary' method). Second derivatives in the Hessian are determined by finite differences of the gradients obtained by automatic differentiation (AD). Return value is a nimbleList of type optimResultNimbleList, similar to what is returned by R's optim. See help(nimOptim). Note that parameters ('par') are returned for the natural parameters, i.e. how they are defined in the model. But the 'hessian', if requested, is computed for the parameters as transformed for optimization if necessary. Hence one must be careful interpreting 'hessian' if any parameters have constraints, and the safest next step is to use the 'summary' method or 'summaryLaplace' function.

• summary(MLEoutput, originalScale, randomEffectsStdError, jointCovariance). Summarize the maximum likelihood estimation results, given object MLEoutput that was returned by findMLE. The summary can include a covariance matrix for the parameters, the random effects, or both), and these can be returned on the original parameter scale or on the (potentially) transformed scale(s) used in estimation. It is often preferred instead to call function (not method) 'summaryLaplace' because this will attach parameter and random effects names (i.e., node names) to the results.

  In more detail, summary accepts the following optional arguments:

  • originalScale. Logical. If TRUE, the function returns results on the original scale(s) of parameters and random effects; otherwise, it returns results on the transformed scale(s). If there are no constraints, the two scales are identical. Defaults to TRUE.

  • randomEffectsStdError. Logical. If TRUE, standard errors of random effects will be calculated. Defaults to FALSE.

  • jointCovariance. Logical. If TRUE, the joint variance-covariance matrix of the parameters and the random effects will be returned. If FALSE, the variance-covariance matrix of the parameters will be returned. Defaults to FALSE.

  The object returned by summary is an AGHQuad_summary nimbleList with elements:

  • params. A nimbleList that contains estimates and standard errors of parameters (on the original or transformed scale, as chosen by originalScale).

  • randomEffects. A nimbleList that contains estimates of random effects and, if requested (randomEffectsStdError=TRUE) their standard errors, on original or transformed scale. Standard errors are calculated following the generalized delta method of Kass and Steffey (1989).

  • vcov. If requested (i.e. jointCovariance=TRUE), the joint variance-covariance matrix of the parameters and random effects, on original or transformed scale. If jointCovariance=FALSE, the covariance matrix of the parameters, on original or transformed scale.

  • scale. "original" or "transformed", the scale on which results were requested.

Methods for more advanced uses

Additional methods to access or control more details of the Laplace approximation include:

• updateSettings. This provides a single function through which many of the settings described above (mostly for the control list) can be later changed. Options that can be changed include: innerOptimMethod, innerOptimStart, innerOptimStartValues, useInnerCache, nQuad, gridType, innerOptimControl, outerOptimControl, and computeMethod. For innerOptimStart, method "zero" cannot be specified but can be achieved by choosing method "constant" with innerOptimStartValues=0. Only provided options will be modified. The exceptions are innerOptimControl, outerOptimControl, which are replaced only replace_innerOptimControl=TRUE or replace_outerOptimControl=TRUE, respectively.

• getNodeNamesVec(returnParams). Return a vector (>1) of names of parameters/random effects nodes, according to returnParams = TRUE/FALSE. Use this if there is more than one node.

• getNodeNameSingle(returnParams). Return the name of a single parameter/random effect node, according to returnParams = TRUE/FALSE. Use this if there is only one node.

• checkInnerConvergence(message). Checks whether all internal optimizers converged. Returns a zero if everything converged and one otherwise. If message = TRUE, it will print more details about convergence for each conditionally independent set.

• gr_logLik(p, trans). Gradient of the (approximated) marginal log-likelihood at parameter value p. Argument trans is similar to that in calcLaplace. If there are multiple parameters, the vector p is given in the order of parameter names returned by getNodeNamesVec(returnParams=TRUE).

• gr_Laplace(p, trans). This is the same as gr_logLik.

• otherLogLik(p). Calculate the calcNodesOther nodes, which returns the log-likelihood of the parts of the model that are not included in the Laplace or AGHQ approximation.

• gr_otherLogLik(p). Gradient (vector of derivatives with respect to each parameter) of otherLogLik(p). Results should match gr_otherLogLik_internal(p) but may be more efficient after the first call.

Internal or development methods

Some methods are included for calculating the (approximate) marginal log posterior density by including the prior distribution of the parameters. This is useful for finding the maximum a posteriori probability (MAP) estimate. Currently these are provided for point calculations without estimation methods.

• calcPrior_p(p). Log density of prior distribution.

• calcPrior_pTransformed(pTransform). Log density of prior distribution on transformed scale, includes the Jacobian.

• calcPostLogDens(p). Marginal log posterior density in terms of the parameter p.

• calcPostLogDens_pTransformed (pTransform). Marginal log posterior density in terms of the transformed parameter, which includes the Jacobian transformation.

• gr_postLogDens_pTransformed(pTransform). Graident of marginal log posterior density on the transformed scale. Other available options that are used in the derivative for more flexible include logDetJacobian(pTransform) and gr_logDeJacobian(pTransform), as well as gr_prior(p).

Finally, methods that are primarily for internal use by other methods include:

• gr_logLik_pTransformed. Gradient of the Laplace approximation (calcLogLik_pTransformed(pTransform)) at transformed (unconstrained) parameter value pTransform.

• pInverseTransform(pTransform). Back-transform the transformed parameter value pTransform to original scale.

• derivs_pInverseTransform(pTransform, order). Derivatives of the back-transformation (i.e. inverse of parameter transformation) with respect to transformed parameters at pTransform. Derivative order is given by order (any of 0, 1, and/or 2).

• reInverseTransform(reTrans). Back-transform the transformed random effects value reTrans to original scale.

• derivs_reInverseTransform(reTrans, order). Derivatives of the back-transformation (i.e. inverse of random effects transformation) with respect to transformed random effects at reTrans. Derivative order is given by order (any of 0, 1, and/or 2).

• optimRandomEffects(pTransform). Calculate the optimized random effects given transformed parameter value pTransform. The optimized random effects are the mode of the conditional distribution of random effects given data at parameters pTransform, i.e. the calculation of calcNodes.

• inverse_negHess(p, reTransform). Calculate the inverse of the negative Hessian matrix of the joint (parameters and random effects) log-likelihood with respect to transformed random effects, evaluated at parameter value p and transformed random effects reTransform.

• hess_logLik_wrt_p_wrt_re(p, reTransform). Calculate the Hessian matrix of the joint log-likelihood with respect to parameters and transformed random effects, evaluated at parameter value p and transformed random effects reTransform.

• one_time_fixes(). Users never need to run this. Is is called when necessary internally to fix dimensionality issues if there is only one parameter in the model.

• calcLogLik_pTransformed(pTransform). Laplace approximation at transformed (unconstrained) parameter value pTransform. To make maximizing the Laplace likelihood unconstrained, an automated transformation via parameterTransform is performed on any parameters with constraints indicated by their priors (even though the prior probabilities are not used).

• gr_otherLogLik_internal(p). Gradient (vector of derivatives with respect to each parameter) of otherLogLik(p). This is obtained using automatic differentiation (AD) with single-taping. First call will always be slower than later calls.

• cache_outer_logLik(logLikVal). Save the marginal log likelihood value to the inner Laplace mariginlization functions to track the outer maximum internally.

• reset_outer_inner_logLik(). Reset the internal saved maximum marginal log likelihood.

• get_inner_cholesky(atOuterMode = integer(0, default = 0)). Returns the cholesky of the negative Hessian with respect to the random effects. If atOuterMode = 1 then returns the value at the overall best marginal likelihood value, otherwise atOuterMode = 0 returns the last.

• get_inner_mode(atOuterMode = integer(0, default = 0)). Returns the mode of the random effects for either the last call to the innner quadrature functions (atOuterMode = 0 ), or the last best value for the marginal log likelihood, atOuterMode = 1.

Author(s)

Wei Zhang, Perry de Valpine, Paul van Dam-Bates

References

Kass, R. and Steffey, D. (1989). Approximate Bayesian inference in conditionally independent hierarchical models (parametric empirical Bayes models). Journal of the American Statistical Association, 84(407), 717-726.

Liu, Q. and Pierce, D. A. (1994). A Note on Gauss-Hermite Quadrature. Biometrika, 81(3) 624-629.

Jackel, P. (2005). A note on multivariate Gauss-Hermite quadrature. London: ABN-Amro. Re.

Skaug, H. and Fournier, D. (2006). Automatic approximation of the marginal likelihood in non-Gaussian hierarchical models. Computational Statistics & Data Analysis, 56, 699-709.

Examples

pumpCode <- nimbleCode({ 
  for (i in 1:N){
    theta[i] ~ dgamma(alpha, beta)
    lambda[i] <- theta[i] * t[i]
    x[i] ~ dpois(lambda[i])
  }
  alpha ~ dexp(1.0)
  beta ~ dgamma(0.1, 1.0)
})
pumpConsts <- list(N = 10, t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5))
pumpData <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22))
pumpInits <- list(alpha = 0.1, beta = 0.1, theta = rep(0.1, pumpConsts$N))
pump <- nimbleModel(code = pumpCode, name = "pump", constants = pumpConsts, 
                    data = pumpData, inits = pumpInits, buildDerivs = TRUE)
                    
# Build Laplace approximation
pumpLaplace <- buildLaplace(pump)

## Not run: 
# Compile the model
Cpump <- compileNimble(pump)
CpumpLaplace <- compileNimble(pumpLaplace, project = pump)
# Calculate MLEs of parameters
MLEres <- CpumpLaplace$findMLE()
# Calculate estimates and standard errors for parameters and random effects on original scale
allres <- CpumpLaplace$summary(MLEres, randomEffectsStdError = TRUE)

# Change the settings and also illustrate runLaplace
CpumpLaplace$updateSettings(innerOptimMethod = "nlminb", outerOptimMethod = "nlminb")
newres <- runLaplace(CpumpLaplace)

# Illustrate use of the component log likelihood and gradient functions to
# run an optimizer manually from R.
# Use nlminb to find MLEs
MLEres.manual <- nlminb(c(0.1, 0.1),
                        function(x) -CpumpLaplace$calcLogLik(x),
                        function(x) -CpumpLaplace$gr_Laplace(x))

## End(Not run)

Description

This function has been moved to the 'nimbleSMC' package.

Usage

buildLiuWestFilter(...)

Arguments

Description

Takes a NIMBLE model (with some missing data, aka random effects or latent state nodes) and builds a Monte Carlo Expectation Maximization (MCEM) algorithm for maximum likelihood estimation. The user can specify which latent nodes are to be integrated out in the E-Step, or default choices will be made based on model structure. All other stochastic non-data nodes will be maximized over. The E-step is done with a sample from a nimble MCMC algorithm. The M-step is done by a call to optim.

Usage

buildMCEM(
  model,
  paramNodes,
  latentNodes,
  calcNodes,
  calcNodesOther,
  control = list(),
  ...
)

Arguments

Details

buildMCEM is a nimbleFunction that creates an MCEM algorithm for a model and choices (perhaps default) of nodes in different roles in the model. The MCEM can then be compiled for fast execution with a compiled model.

Note that buildMCEM was re-written for nimble version 1.2.0 and is not backward-compatible with previous versions. The new version is considered to be in beta testing.

Denote data by Y, latent states (or missing data) by X, and parameters by T. MCEM works by the following steps, starting from some T:

1. Draw a sample of size M from P(X | Y, T) using MCMC.

2. Update T to be the maximizer of E[log P(X, Y | T)] where the expectation is approximated as a Monte Carlo average over the sample from step(1)

3. Repeat until converged.

The default version of MCEM is the ascent-based MCEM of Caffo et al. (2015). This attempts to update M when necessary to ensure that step 2 really moves uphill given that it is maximizing a Monte Carlo approximation and could accidently move downhill on the real surface of interest due to Monte Carlo error. The main tuning parameters include alpha, beta, gamma, Mfactor, C, and tol (tolerance).

If the model supports derivatives via nimble's automatic differentiation (AD) (and buildDerivs=TRUE in nimbleModel), the maximization step can use gradients from AD. You must manually set useDerivs=FALSE in the control list if derivatives aren't supported or if you don't want to use them.

In the ascent-based method, after maximization in step 2, the Monte Carlo standard error of the uphill movement is estimated. If the standardized uphill step is bigger than 0 with Type I error rate alpha, the iteration is accepted and the algorithm continues. Otherwise, it is not certain that step 2 really moved uphill due to Monte Carlo error, so the MCMC sample size M is incremented by a fixed factor (e.g. 0.33 or 0.5, called Mfactor in the control list), the additional samples are added by continuing step 1, and step 2 is tried again. If the Monte Carlo noise still overwhelms the magnitude of uphill movement, the sample size is increased again, and so on. alpha should be between 0 and 0.5. A larger value than usually used for inference is recommended so that there is an easy threshold to determine uphill movement, which avoids increasing M prematurely. M will never be increased above maxM.

Convergence is determined in a similar way. After a definite move uphill, we determine if the uphill increment is less than tol, with Type I error rate gamma. (But if M hits a maximum value, the convergence criterion changes. See below.)

beta is used to help set M to a minimal level based on previous iterations. This is a desired Type II error rate, assuming an uphill move and standard error based on the previous iteration. Set adjustM=FALSE in the control list if you don't want this behavior.

There are some additional controls on convergence for practical purposes. Set C in the control list to be the number of times the convergence criterion mut be satisfied in order to actually stop. E.g setting C=2 means there will always be a restart after the first convergence.

One problem that can occur with ascent-based MCEM is that the final iteration can be very slow if M must become very large to satisfy the convergence criterion. Indeed, if the algorithm starts near the MLE, this can occur. Set maxM in the control list to set the MCMC sample size that should never be exceeded.

If M==maxM, a softer convergence criterion is used. This second convergence criterion is to stop if we can't be sure we moved uphill using Type I error rate delta. This is a soft criterion because for small delta, Type II errors will be common (e.g. if we really did move uphill but can't be sure from the Monte Carlo sample), allowing the algorithm to terminate. One can continue the algorithm from where it stopped, so it is helpful to not have it get stuck when having a very hard time with the first (stricter) convergence criterion.

All of alpha, beta, delta, and gamma are utilized based on asymptotic arguments but in practice must be chosen heuristically. In other words, their theoretical basis does not exactly yield practical advice on good choices for efficiency and accuracy, so some trial and error will be needed.

It can also be helpful to set a minimum and maximum of allowed iterations (of steps 1 and 2 above). Setting minIter>1 in the control list can sometimes help avoid a false convergence on the first iteration by forcing at least one more iteration. Setting maxIter provides a failsafe on a stuck run.

If you don't want the ascent-based method at all and simply want to run a set of iterations, set ascent=FALSE in the control list. This will use the second (softer) convergence criterion.

Parameters to be maximized will by default be handled in an unconstrained parameter space, transformed if necessary by a parameterTransform object. In that case, the default optim method will be "BFGS" and can can be changed by setting optimMehod in the control list. Set useTransform=FALSE in the control list if you don't want the parameters transformed. In that case the default optimMethod will be "L-BFGS-B" if there are any actual constraints, and you can provide a list of boxConstraints in the control list. (Constraints may be determined by priors written in the model for parameters, even though their priors play no other role in MLE. E.g. sigma ~ halfflat() indicates sigma > 0).

Most of the control list elements can be overridden when calling the findMLE method. The findMLE argument continue=TRUE results in attempting to continue the algorithm where the previous call finished, including whatever settings were in use.

See setupMargNodes (which is called with the given arguments for paramNodes, calcNodes, and calcNodesOther; and with allowDiscreteLatent=TRUE, randomEffectsNodes=latentNodes, and check=check) for more about how the different groups of nodes are determined. In general, you can provide none, one, or more of the different kinds of nodes and setupMargNodes will try to determine the others in a sensible way. However, note that this cannot work for all ways of writing a model. One key example is that if random (latent) nodes are written as top-level nodes (e.g. following N(0,1)), they appear structurally to be parameters and you must tell buildMCEM that they are latentNodes. The various "Nodes" arguments will all be passed through model$expandNodeNames, allowing for example simply "x" to be provided when there are many nodes within "x".

Estimating the Monte Carlo standard error of the uphill step is not trivial because the sample was obtained by MCMC and so likely is autocorrelated. This is done by calling whatever function in R's global environment is called "MCEM_mcse", which is required to take two arguments: samples (which will be a vector of the differences in log(P(Y, X | T)) between the new and old values of T, across the sample of X) and m, the sample size. It must return an estimate of the standard error of the mean of the sample. NIMBLE provides a default version (exported from the package namespace), which calls mcmcse::mcse with method "obm". Simply provide a different function with this name in your R session to override NIMBLE's default.

Control list details

The control list accepts the following named elements:

• initM initial MCMC sample size, M. Default=1000.

• Mfactor Factor by which to increase MCMC sample size when step 2 results in noise overwhelming the uphill movement. The new M will be 1+Mfactor)*M (rounded up). Mfactor is 1/k of Caffo et al. (2015). Default=1/3.

• maxM Maximum allowed value of M (see above). Default=initM*20.

• burnin Number of burn-in iterations for the MCMC in step 1. Note that the initial states of one MCMC will be the last states from the previous MCMC, so they will often be good initial values after multiple iterations. Default=500.

• thin Thinning interval for the MCMC in step 1. Default=1.

• alpha Type I error rate for determining when step 2 has moved uphill. See above. Default=0.25.

• beta Used for determining a minimal value of $M$ based on previous iteration, if adjustM is TRUE. beta is a desired Type II error rate for determining uphill moves. Default=0.25.

• delta Type I error rate for the soft convergence approach (second approach above). Default=0.25.

• gamma Type I error rate for determining when step 2 has moved less than tol uphill, in which case ascent-based convergence is achieved (first approach above). Default=0.05.

• buffer A small amount added to lower box constraints and substracted from upper box constraints for all parameters, relevant only if useTransform=FALSE and some parameters do have boxConstraints set or have bounds that can be determined from the model. Default=1e-6.

• tol Ascent-based convergence tolerance. Default=0.001.

• ascent Logical to determine whether to use the ascent-based method of Caffo et al. Default=TRUE.

• C Number of convergences required to actually stop the algorithm. Default = 1.

• maxIter Maximum number of MCEM iterations to run.

• minIter Minimum number of MCEM iterations to run.

• adjustM Logical indicating whether to see if M needs to be increased based on statistical power argument in each iteration (using beta). Default=TRUE.

• verbose Logical indicating whether verbose output is desired. Default=TRUE.

• MCMCprogressBar Logical indicating whether MCMC progress bars should be shown for every iteration of step 1. This argument is passed to configureMCMC, or to config if provided. Default=TRUE.

• derivsDelta If AD derivatives are not used, then the method vcov must use finite difference derivatives to implement the method of Louis (1982). The finite differences will be delta or delta/2 for various steps. This is the same for all dimensions. Default=0.0001.

• mcmcControl This is passed to configureMCMC, or config if provided, as the control argument. i.e. control=mcmcControl.

• boxContrainst List of box constraints for the nodes that will be maximized over, only relevant if useTransform=FALSE and forceNoConstraints=FALSE (and ignored otherwise). Each constraint is a list in which the first element is a character vector of node names to which the constraint applies and the second element is a vector giving the lower and upper limits. Limits of -Inf or Inf are allowed. Any nodes that are not given constrains will have their constraints automatically determined by NIMBLE. See above. Default=list().

• forceNoConstraints Logical indicating whether to force ignoring constraints even if they might be necessary. Default=FALSE.

• useTransform Logical indicating whether to use a parameter transformation (see parameterTransform) to create an unbounded parameter space for the paramNodes. This allows unconstrained maximization algorithms to be used. Default=TRUE.

• check Logical passed as the check argument to setupMargNodes. Default=TRUE.

• useDerivs Logical indicating whether to use AD. If TRUE, the model must have been build with 'buildDerivs=TRUE'. It is not automatically determined from the model whether derivatives are supported. Default=TRUE.

• config Function to create the MCMC configuration used for step 1. The MCMC configuration is created by calling

  config(model, nodes = latentNodes, monitors = latentNodes,
  thin = thinDefault, control = mcmcControl, print = FALSE) 

  The default for config (if it is missing) is configureMCMC, which is nimble's general default MCMC configuration function.

Methods in the returned algorithm

The object returned by buildMCEM is a nimbleFunction object with the following methods

• findMLE is the main method of interest, launching the MCEM algorithm. It takes the following arguments:

  • pStart. Vector of initial parameter values. If omitted, the values currently in the model object are used.

  • returnTrans. Logical indicating whether to return parameters in the transformed space, if a parameter transformation is in use. Default=FALSE.

  • continue. Logical indicating whether to continue the MCEM from where the last call stopped. In addition, if TRUE, any other control setting provided in the last call will be used again. If FALSE, all control settings are reset to the values provided when buildMCEM was called. Any control settings provided in the same call as continue=TRUE will over-ride these behaviors and be used in the continued run.

  • All run-time control settings available in the control list for buildMCEM (except for buffer, boxConstraints, forceNoConstraints, useTransform, and useDerivs) are accepted as individual arguments to over-ride the values provided in the control list.

  findMLE returns on object of class optimResultNimbleList with the results of the final optimization of step 2. The par element of this list is the vector of maximum likelihood (MLE) parameters.

• vcov computes the approximate variance-covariance matrix of the MLE using the method of Louis (1982). It takes the following arguments:

  • params. Vector of parameters at which to compute the Hessian matrix used to obtain the vcov result. Typically this will be MLE$par, if MLE is the output of findMLE.

  • trans. Logical indicating whether params is on the transformed prameter scale, if a parameter transformation is in use. Typically this should be the same as the returnTrans argument to findMLE. Default=FALSE.

  • returnTrans. Logical indicting whether the vcov result should be for the transformed parameter space. Default matches trans.

  • M. Number of MCMC samples to obtain if resetSamples=TRUE. Default is the final value of M from the last call to findMLE. It can be helpful to increase M to obtain a more accurate vcov result (i.e. with less Monte Carlo noise).

  • resetSamples. Logical indicating whether to generate a new MCMC sample from P(X | Y, T), where T is params. If FALSE, the last sample from findMLE will be used. If MLE convergence was reasonable, this sample can be used. However, if the last MCEM step made a big move in parameter space (e.g. if convergence was not achieved), the last MCMC sample may not be accurate for obtaining vcov. Default=FALSE.

  • atMLE. Logical indicating whether you believe the params represents the MLE. If TRUE, one part of the computation will be skipped because it is expected to be 0 at the MLE. If there are parts of the model that are not connected to the latent nodes, i.e. of calcNodesOther is not empty, then atMLE will be ignored and set to FALSE. Default=FALSE. It is not really worth using TRUE unless you are confident and the time saving is meaningful, which is not very likely. In other words, this argument is provided for technical completeness.

  vcov returns a matrix that is the inverse of the negative Hessian of the log likelihood surface, i.e. the usual asymptotic approximation of the parameter variance-covariance matrix.

• doMCMC. This method runs the MCMC to sample from P(X | Y, T). One does not need to call this, as it is called via the MCEM algorithm in findMLE. This method is provided for users who want to use the MCMC for latent states directly. Samples should be retrieved by as.matrix(MCEM$mvSamples), where MCEM is the (compiled or uncompiled) MCEM algorithm object. This method takes the following arguments:

  • M. MCMC sample size.

  • thin. MCMC thinning interval.

  • reset. Logical indicating whether to reset the MCMC (passed to the MCMC run method as reset).

• transform and inverseTransform. Convert a parameter vector to an unconstrained parameter space and vice-versa, if useTransform=TRUE in the call to buildDerivs.

• resetControls. Reset all control arguments to the values provided in the call to buildMCEM. The user does not normally need to call this.

Author(s)

Perry de Valpine, Clifford Anderson-Bergman and Nicholas Michaud

References

Caffo, Brian S., Wolfgang Jank, and Galin L. Jones (2005). Ascent-based Monte Carlo expectation-maximization. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 67(2), 235-251.

Louis, Thomas A (1982). Finding the Observed Information Matrix When Using the EM Algorithm. Journal of the Royal Statistical Society. Series B (Statistical Methodology), 44(2), 226-233.

Examples

## Not run: 
pumpCode <- nimbleCode({
 for (i in 1:N){
     theta[i] ~ dgamma(alpha,beta);
     lambda[i] <- theta[i]*t[i];
     x[i] ~ dpois(lambda[i])
 }
 alpha ~ dexp(1.0);
 beta ~ dgamma(0.1,1.0);
})

pumpConsts <- list(N = 10,
              t = c(94.3, 15.7, 62.9, 126, 5.24,
                31.4, 1.05, 1.05, 2.1, 10.5))

pumpData <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22))

pumpInits <- list(alpha = 1, beta = 1,
             theta = rep(0.1, pumpConsts$N))
pumpModel <- nimbleModel(code = pumpCode, name = 'pump', constants = pumpConsts,
                         data = pumpData, inits = pumpInits,
                         buildDerivs=TRUE)

pumpMCEM <- buildMCEM(model = pumpModel)

CpumpModel <- compileNimble(pumpModel)

CpumpMCEM <- compileNimble(pumpMCEM, project=pumpModel)

MLE <- CpumpMCEM$findMLE()
vcov <- CpumpMCEM$vcov(MLE$par)


## End(Not run)

Description

First required argument, which may be of class MCMCconf (an MCMC configuration object), or inherit from class modelBaseClass (a NIMBLE model object). Returns an uncompiled executable MCMC object. See details.

Usage

buildMCMC(conf, print, ...)

Arguments

Details

Calling buildMCMC(conf) will produce an uncompiled MCMC object. The object contains several methods, including the main run function for running the MCMC, a getTimes function for determining the computation time spent in each sampler (see 'getTimes' section below), and functions related to WAIC (getWAIC, getWAICdetails, calculateWAIC (see help(waic)).

The uncompiled run function will have arguments:

niter: The number of iterations to run the MCMC.

nburnin: Number of initial, pre-thinning, MCMC iterations to discard (default = 0).

thin: The thinning interval for the monitors that were specified in the MCMC configuration. If this argument is provided at MCMC runtime, it will take precedence over the thin interval that was specified in the MCMC configuration. If omitted, the thin interval from the MCMC configuration will be used.

thin2: The thinning interval for the second set of monitors (monitors2) that were specified in the MCMC configuration. If this argument is provided at MCMC runtime, it will take precedence over the thin2 interval that was specified in the MCMC configuration. If omitted, the thin2 interval from the MCMC configuration will be used.

reset: Boolean specifying whether to reset the internal MCMC sampling algorithms to their initial state (in terms of self-adapting tuning parameters), and begin recording posterior sample chains anew. Specifying reset = FALSE allows the MCMC algorithm to continue running from where it left off, appending additional posterior samples to the already existing sample chains. Generally, reset = FALSE should only be used when the MCMC has already been run (default = TRUE).

resetMV: Boolean specifying whether to begin recording posterior sample chains anew. This argument is only considered when using reset = FALSE. Specifying reset = FALSE, resetMV = TRUE allows the MCMC algorithm to continue running from where it left off, but without appending the new posterior samples to the already existing samples, i.e. all previously obtained samples will be erased. This option can help reduce memory usage during burn-in (default = FALSE).

resetWAIC: Boolean specifying whether to reset the WAIC summary statistics to their initial states and thereby begin the WAIC calculation anew (default = TRUE). Specifying resetWAIC = FALSE allows the WAIC calculation to continue running from where it left off.

initializeModel: Boolean specifying whether to run the initializeModel routine on the underlying model object, prior to beginning MCMC sampling (default = TRUE).

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). time: Boolean specifying whether to record runtimes of the individual internal MCMC samplers. When time = TRUE, a vector of runtimes (measured in seconds) can be extracted from the MCMC using the method mcmc$getTimes() (default = FALSE).

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 nimbleOptions(MCMCprogressBar = FALSE).

Samples corresponding to the monitors and monitors2 from the MCMCconf are stored into the interval variables mvSamples and mvSamples2, respectively. These may be accessed and converted into R matrix or list objects via: as.matrix(mcmc$mvSamples) as.list(mcmc$mvSamples) as.matrix(mcmc$mvSamples2) as.list(mcmc$mvSamples2)

The uncompiled MCMC function may be compiled to a compiled MCMC object, taking care to compile in the same project as the R model object, using: Cmcmc <- compileNimble(Rmcmc, project = Rmodel)

The compiled object will function identically to the uncompiled object except acting on the compiled model object.

Timing the MCMC samplers

If you want to obtain the computation time spent in each sampler, you can set time=TRUE as a run-time argument to run() and then use the method getTimes() to obtain the times.

Calculating WAIC

Please see help(waic) for more information.

Author(s)

Daniel Turek

References

Watanabe, S. (2010). Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory. Journal of Machine Learning Research 11: 3571-3594.

Gelman, A., Hwang, J. and Vehtari, A. (2014). Understanding predictive information criteria for Bayesian models. Statistics and Computing 24(6): 997-1016.

Ariyo, O., Quintero, A., Munoz, J., Verbeke, G. and Lesaffre, E. (2019). Bayesian model selection in linear mixed models for longitudinal data. Journal of Applied Statistics 47: 890-913.

See Also

configureMCMC runMCMC nimbleMCMC

Examples

## Not run: 
code <- nimbleCode({
    mu ~ dnorm(0, 1)
    x ~ dnorm(mu, 1)
    y ~ dnorm(x, 1)
})
Rmodel <- nimbleModel(code, data = list(y = 0))
conf <- configureMCMC(Rmodel, monitors = c('mu', 'x'), enableWAIC = TRUE)
Rmcmc <- buildMCMC(conf)
Cmodel <- compileNimble(Rmodel)
Cmcmc <- compileNimble(Rmcmc, project=Rmodel)

## Running the MCMC with `run`
Cmcmc$run(10000)
samples <- as.matrix(Cmcmc$mvSamples)
samplesAsList <- as.list(Cmcmc$mvSamples)
head(samples)

## Getting WAIC
waicInfo <- Cmcmc$getWAIC()
waicInfo$WAIC
waicInfo$pWAIC

## Timing the samplers (must set `time = TRUE` when running the MCMC)
Cmcmc$run(10000, time = TRUE)
Cmcmc$getTimes()

## End(Not run)

Description

In addition to the core online algorithm, NIMBLE implements an offline WAIC algorithm that can be computed on the results of an MCMC. In contrast to NIMBLE's built-in online WAIC, offline WAIC can compute only conditional WAIC and does not allow for grouping data nodes.

Usage

calculateWAIC(mcmc, model, nburnin = 0, thin = 1)

Arguments

Details

The ability to calculate WAIC post hoc after all MCMC sampling has been done has certain advantages (e.g., allowing a user to calculate WAIC from MCMC chains run separately) in addition to providing compatibility with versions of NIMBLE before 0.12.0. This functionality includes the ability to call the calculateWAIC function on an MCMC object or matrix of samples after running an MCMC and without setting up the MCMC initially to use WAIC.

Important: The necessary variables to compute WAIC (all stochastic parent nodes of the data nodes) must have been monitored when setting up the MCMC.

Also note that while the model argument can be either a compiled or uncompiled model, the model must have been compiled prior to calling calculateWAIC.

See help(waic) for details on using NIMBLE's recommended online algorithm for WAIC.

Offline WAIC (WAIC computed after MCMC sampling)

As an alternative to online WAIC, NIMBLE also provides a function, calculateWAIC, that can be called on an MCMC object or a matrix of samples, after running an MCMC. This function does not require that one set enableWAIC = TRUE nor WAIC = TRUE when calling runMCMC. The function checks that the necessary variables were monitored in the MCMC and returns an error if they were not. This function behaves identically to the calculateWAIC method of an MCMC object. Note that to use this function when using nimbleMCMC one would need to build the model outside of nimbleMCMC.

The calculateWAIC function requires either an MCMC object or a matrix (or dataframe) of posterior samples plus a model object. In addition, one can provide optional burnin and thin arguments.

In addition, for compatibility with older versions of NIMBLE (prior to v0.12.0), one can also use the calculateWAIC method of the MCMC object to calculate WAIC after all sampling has been completed.

The calculateWAIC() method accepts a single argument, nburnin, equivalent to the nburnin argument of the calculateWAIC function described above.

The calculateWAIC method can only be used if the enableWAIC argument to configureMCMC or to buildMCMC is set to TRUE, or if the NIMBLE option enableWAIC is set to TRUE. If a user attempts to call calculateWAIC without having set enableWAIC = TRUE (either in the call to configureMCMC, or buildMCMC, or as a NIMBLE option), an error will occur.

The calculateWAIC function and method calculate the WAIC based on Equations 5, 12, and 13 in Gelman et al. (2014) (i.e., using pWAIC2).

Note that there is not a unique value of WAIC for a model. The calculateWAIC function and method only provide the conditional WAIC, namely the version of WAIC where all parameters directly involved in the likelihood are treated as $theta$ for the purposes of Equation 5 from Gelman et al. (2014). As a result, the user must set the MCMC monitors (via the monitors argument) to include all stochastic nodes that are parents of any data nodes; by default the MCMC monitors are only the top-level nodes of the model. For more detail on the use of different predictive distributions, see Section 2.5 from Gelman et al. (2014) or Ariyo et al. (2019). Also note that WAIC relies on a partition of the observations, i.e., 'pointwise' prediction. In calculateWAIC the sum over log pointwise predictive density values treats each data node as contributing a single value to the sum. When a data node is multivariate, that data node contributes a single value to the sum based on the joint density of the elements in the node. Note that if one wants the WAIC calculation via calculateWAIC to be based on the joint predictive density for each group of observations (e.g., grouping the observations from each person or unit in a longitudinal data context), one would need to use a multivariate distribution for the observations in each group (potentially by writing a user-defined distribution).

For more control over and flexibility in how WAIC is calculated, see help(waic).

Author(s)

Joshua Hug and Christopher Paciorek

References

Watanabe, S. (2010). Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory. Journal of Machine Learning Research 11: 3571-3594.

Gelman, A., Hwang, J. and Vehtari, A. (2014). Understanding predictive information criteria for Bayesian models. Statistics and Computing 24(6): 997-1016.

Ariyo, O., Quintero, A., Munoz, J., Verbeke, G. and Lesaffre, E. (2019). Bayesian model selection in linear mixed models for longitudinal data. Journal of Applied Statistics 47: 890-913.

Vehtari, A., Gelman, A. and Gabry, J. (2017). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing 27: 1413-1432.

Hug, J.E. and Paciorek, C.J. (2021). A numerically stable online implementation and exploration of WAIC through variations of the predictive density, using NIMBLE. arXiv e-print <arXiv:2106.13359>.

See Also

waic configureMCMC buildMCMC runMCMC nimbleMCMC

Examples

code <- nimbleCode({
  for(j in 1:J) {
    for(i in 1:n) 
      y[j, i] ~ dnorm(mu[j], sd = sigma)
    mu[j] ~ dnorm(mu0, sd = tau)
  }
  tau ~ dunif(0, 10)
  sigma ~ dunif(0, 10)
})
J <- 5
n <- 10
y <- matrix(rnorm(J*n), J, n)
Rmodel <- nimbleModel(code, constants = list(J = J, n = n), data = list(y = y),
                      inits = list(tau = 1, sigma = 1))

## Make sure the needed variables are monitored.
## Only conditional WAIC without data grouping is available via this approach.
conf <- configureMCMC(Rmodel, monitors = c('mu', 'sigma'))
## Not run: 
Cmodel <- compileNimble(Rmodel)
Rmcmc <- buildMCMC(conf)
Cmcmc <- compileNimble(Rmcmc, project = Rmodel)
output <- runMCMC(Cmcmc, niter = 1000)
calculateWAIC(Cmcmc)           # Can run on the MCMC object
calculateWAIC(output, Rmodel)  # Can run on the samples directly

## Apply additional burnin (additional to any burnin already done in the MCMC.
calculateWAIC(Cmcmc, burnin = 500)

## End(Not run)

Description

Calculate number of islands (distinct connected groups) based on a CAR adjacency matrix.

Usage

CAR_calcNumIslands(adj, num)

Arguments

Author(s)

Daniel Turek

See Also

CAR-Normal

Description

Density function and random generation for the improper (intrinsic) Gaussian conditional autoregressive (CAR) distribution.

Usage

dcar_normal(
  x,
  adj,
  weights = adj/adj,
  num,
  tau,
  c = CAR_calcNumIslands(adj, num),
  zero_mean = 0,
  log = FALSE
)

rcar_normal(
  n = 1,
  adj,
  weights = adj/adj,
  num,
  tau,
  c = CAR_calcNumIslands(adj, num),
  zero_mean = 0
)

Arguments

Details

When specifying a CAR distribution in BUGS model code, the zero_mean parameter should be specified as either 0 or 1 (rather than TRUE or FALSE).

Note that because the distribution is improper, rcar_normal does not generate a sample from the distribution. However, as discussed in Rue and Held (2005), it is possible to generate a sample from the distribution under constraints imposed based on the eigenvalues of the precision matrix that are zero.

Value

dcar_normal gives the density, while rcar_normal returns the current process values, since this distribution is improper.

Author(s)

Daniel Turek

References

Banerjee, S., Carlin, B.P., and Gelfand, A.E. (2015). Hierarchical Modeling and Analysis for Spatial Data, 2nd ed. Chapman and Hall/CRC.

Rue, H. and L. Held (2005). Gaussian Markov Random Fields, Chapman and Hall/CRC.

See Also

CAR-Proper, Distributions for other standard distributions

Examples

x <- c(1, 3, 3, 4)
num <- c(1, 2, 2, 1)
adj <- c(2, 1,3, 2,4, 3)
weights <- c(1, 1, 1, 1, 1, 1)
lp <- dcar_normal(x, adj, weights, num, tau = 1)

Description

Density function and random generation for the proper Gaussian conditional autoregressive (CAR) distribution.

Usage

dcar_proper(
  x,
  mu,
  C = CAR_calcC(adj, num),
  adj,
  num,
  M = CAR_calcM(num),
  tau,
  gamma,
  evs = CAR_calcEVs3(C, adj, num),
  log = FALSE
)

rcar_proper(
  n = 1,
  mu,
  C = CAR_calcC(adj, num),
  adj,
  num,
  M = CAR_calcM(num),
  tau,
  gamma,
  evs = CAR_calcEVs3(C, adj, num)
)

Arguments

Details

If both C and M are omitted, then all weights are taken as one, and corresponding values of C and M are generated.

The C and M parameters must jointly satisfy a symmetry constraint: that M^(-1) %*% C is symmetric, where M is a diagonal matrix and C is the full weight matrix that is sparsely represented by the parameter vector C.

For a proper CAR model, the value of gamma must lie within the inverse minimum and maximum eigenvalues of M^(-0.5) %*% C %*% M^(0.5), where M is a diagonal matrix and C is the full weight matrix. These bounds can be calculated using the deterministic functions carMinBound(C, adj, num, M) and carMaxBound(C, adj, num, M), or simultaneously using carBounds(C, adj, num, M). In the case where C and M are omitted (all weights equal to one), the bounds on gamma are necessarily (-1, 1).

Value

dcar_proper gives the density, and rcar_proper generates random deviates.

Author(s)

Daniel Turek

References

Banerjee, S., Carlin, B.P., and Gelfand, A.E. (2015). Hierarchical Modeling and Analysis for Spatial Data, 2nd ed. Chapman and Hall/CRC.

See Also

CAR-Normal, Distributions for other standard distributions

Examples

x <- c(1, 3, 3, 4)
mu <- rep(3, 4)
adj <- c(2, 1,3, 2,4, 3)
num <- c(1, 2, 2, 1)
 
## omitting C and M uses all weights = 1
dcar_proper(x, mu, adj = adj, num = num, tau = 1, gamma = 0.95)
 
## equivalent to above: specifying all weights = 1,
## then using as.carCM to generate C and M arguments
weights <- rep(1, 6)
CM <- as.carCM(adj, weights, num)
C <- CM$C
M <- CM$M
dcar_proper(x, mu, C, adj, num, M, tau = 1, gamma = 0.95)
 
## now using non-unit weights
weights <- c(2, 2, 3, 3, 4, 4)
CM2 <- as.carCM(adj, weights, num)
C2 <- CM2$C
M2 <- CM2$M
dcar_proper(x, mu, C2, adj, num, M2, tau = 1, gamma = 0.95)

Description

Calculate the lower and upper bounds for the gamma parameter of the dcar_proper distribution

Usage

carBounds(C, adj, num, M)

Arguments

Details

Bounds for gamma are the inverse of the minimum and maximum eigenvalues of: $M^(-0.5) C M^(0.5)$. The lower and upper bounds are returned in a numeric vector.

Value

A numeric vector containing the bounds (minimum and maximum allowable values) for the gamma parameter of the dcar_proper distribution.

Author(s)

Daniel Turek

See Also

CAR-Proper, carMinBound, carMaxBound

Description

Calculate the upper bound for the gamma parameter of the dcar_proper distribution

Usage

carMaxBound(C, adj, num, M)

Arguments

Details

Bounds for gamma are the inverse of the minimum and maximum eigenvalues of $M^(-0.5) C M^(0.5)$.

Value

The upper bound (maximum allowable value) for the gamma parameter of the dcar_proper distribution.

Author(s)

Daniel Turek

See Also

CAR-Proper, carMinBound, carBounds

Description

Calculate the lower bound for the gamma parameter of the dcar_proper distribution

Usage

carMinBound(C, adj, num, M)

Arguments

Details

Bounds for gamma are the inverse of the minimum and maximum eigenvalues of: $M^(-0.5) C M^(0.5)$.

Value

The lower bound (minimum allowable value) for the gamma parameter of the dcar_proper distribution.

Author(s)

Daniel Turek

See Also

CAR-Proper, carMaxBound, carBounds

Description

Density and random generation for the categorical distribution

Usage

dcat(x, prob, log = FALSE)

rcat(n = 1, prob)

Arguments

Details

See the BUGS manual for mathematical details.

Value

dcat gives the density and rcat generates random deviates.

Author(s)

Christopher Paciorek

See Also

Distributions for other standard distributions

Examples

probs <- c(1/4, 1/10, 1 - 1/4 - 1/10)
x <- rcat(n = 30, probs)
dcat(x, probs)

Description

Check for interrupt (e.g. Ctrl-C) during nimbleFunction execution. Part of the NIMBLE language.

Usage

checkInterrupt()

Details

During execution of nimbleFunctions that take a long time, it is nice to occassionally check if the user has entered an interrupt and bail out of execution if so. This function does that. During uncompiled nimbleFunction execution, it does nothing. During compiled execution, it calls R_checkUserInterrupt() of the R headers.

Author(s)

Perry de Valpine

Description

Density and random generation for the Chinese Restaurant Process distribution.

Usage

dCRP(x, conc = 1, size, log = 0)

rCRP(n, conc = 1, size)

Arguments

Details

The Chinese restaurant process distribution is a distribution on the space of partitions of the positive integers. The distribution with concentration parameter $\alpha$ equal to conc has probability function

$f(x_i \mid x_1, \ldots, x_{i-1})=\frac{1}{i-1+\alpha}\sum_{j=1}^{i-1}\delta_{x_j}+ \frac{\alpha}{i-1+\alpha}\delta_{x^{new}},$

where $x^{new}$ is a new integer not in $x_1, \ldots, x_{i-1}$.

If conc is not specified, it assumes the default value of 1. The conc parameter has to be larger than zero. Otherwise, NaN are returned.

Value

dCRP gives the density, and rCRP gives random generation.

Author(s)

Claudia Wehrhahn

References

Blackwell, D., and MacQueen, J. B. (1973). Ferguson distributions via Pólya urn schemes. The Annals of Statistics, 1: 353-355.

Aldous, D. J. (1985). Exchangeability and related topics. In École d'Été de Probabilités de Saint-Flour XIII - 1983 (pp. 1-198). Springer, Berlin, Heidelberg.

Pitman, J. (1996). Some developments of the Blackwell-MacQueen urn scheme. IMS Lecture Notes-Monograph Series, 30: 245-267.

Examples

x <- rCRP(n=1, conc = 1, size=10)
dCRP(x, conc = 1, size=10)

Description

Clear all compiled objects from a project and unload the shared library produced by the C++ compiler. Has no effect on Windows.

Usage

clearCompiled(obj)

Arguments

Details

This will clear all compiled objects associated with your NIMBLE project. For example, if cModel is a compiled model, clearCompiled(cModel) will clear both the model and all associated nimbleFunctions such as compiled MCMCs that use that model.

Use of this function can be dangerous. There is some risk that if you have copies of the R objects that interfaced to compiled C++ objects that have been removed, and you attempt to use those R objects after clearing their compiled counterparts, you will crash R. We have tried to minimize that risk, but we can't guarantee safe behavior.

Description

Classes used internally in NIMBLE and not expected to be called directly by users.

Description

Classes used internally in NIMBLE and not expected to be called directly by users.

Description

Classes used internally in NIMBLE and not expected to be called directly by users.

Description

compile a collection of models and nimbleFunctions: generate C++, compile the C++, load the result, and return an interface object

Usage

compileNimble(
  ...,
  project,
  dirName = NULL,
  projectName = "",
  control = list(),
  resetFunctions = FALSE,
  showCompilerOutput = getNimbleOption("showCompilerOutput")
)

Arguments

Details

This is the main function for calling the NIMBLE compiler. A set of compiler calls and output will be seen. Compiling in NIMBLE does 4 things: 1. It generates C++ code files for all the model and nimbleFunction components. 2. It calls the system's C++ compiler. 3. It loads the compiled object(s) into R using dyn.load. And 4. it generates R objects for using the compiled model and nimbleFunctions.

When the units for compilation provided in ... include multiple models and/or nimbleFunctions, models are compiled first, in the order in which they are provided. Groups of nimbleFunctions that were specialized from the same nimbleFunction generator (the result of a call to nimbleFunction, which then takes setup arguments and returns a specialized nimbleFunction) are then compiled as a group, in the order of first appearance.

The behavior of adding new compilation units to an existing project is limited. For example, one can compile a model in one call to compileNimble and then compile a nimbleFunction that uses the model (i.e. was given the model as a setup argument) in a second call to compileNimble, with the model provided as the project argument. Either the uncompiled or compiled model can be provided. However, compiling a second nimbleFunction and adding it to the same project will only work in limited circumstances. Basically, the limitations occur because it attempts to re-use already compiled pieces, but if these do not have all the necessary information for the new compilation, it gives up. An attempt has been made to give up in a controlled manner and provide somewhat informative messages.

When compilation is not allowed or doesn't work, try using resetFunctions = TRUE, which will force recompilation of all nimbleFunctions in the new call. Previously compiled nimbleFunctions will be unaffected, and their R interface objects should continue to work. The only cost is additional compilation time for the current compilation call. If that doesn't work, try re-creating the model and/or the nimbleFunctions from their generators. An alternative possible fix is to compile multiple units in one call, rather than sequentially in multiple calls.

The control list can contain the following named elements, each with TRUE or FALSE: debug, which sets a debug mode for the compiler for development purposes; debugCpp, which inserts an output message before every line of C++ code for debugging purposes; compileR, which determines whether the R-only steps of compilation should be executed; writeCpp, which determines whether the C++ files should be generated; compileCpp, which determines whether the C++ should be compiled; loadSO, which determines whether the DLL or shared object should be loaded and interfaced; and returnAsList, which determines whether calls to the compiled nimbleFunction should return only the returned value of the call (returnAsList = FALSE) or whether a list including the input arguments, possibly modified, should be returned in a list with the returned value of the call at the end (returnAsList = TRUE). The control list is mostly for developer use, although returnAsArgs may be useful to a user. An example of developer use is that one can have the compiler write the C++ files but not compile them, then modify them by hand, then have the C++ compiler do the subsequent steps without over-writing the files.

See the NIMBLE User Manual Manual for examples

Value

If there is only one compilation unit (one model or nimbleFunction), an R interface object is returned. This object can be used like the uncompiled model or nimbleFunction, but execution will call the corresponding compiled objects or functions. If there are multiple compilation units, they will be returned as a list of interface objects, in the order provided. If names were included in the arguments, or in a list if any elements of ... are lists, those names will be used for the corresponding element of the returned list. Otherwise an attempt will be made to generate names from the argument code. For example compileNimble(A = fun1, B = fun2, project = myModel) will return a list with named elements A and B, while compileNimble(fun1, fun2, project = myModel) will return a list with named elements fun1 and fun2.

Author(s)

Perry de Valpine

Description

Creates a default MCMC configuration for a given model.

Usage

configureMCMC(
  model,
  nodes,
  control = list(),
  monitors,
  thin = 1,
  monitors2 = character(),
  thin2 = 1,
  useConjugacy = getNimbleOption("MCMCuseConjugacy"),
  onlyRW = FALSE,
  onlySlice = FALSE,
  multivariateNodesAsScalars = getNimbleOption("MCMCmultivariateNodesAsScalars"),
  enableWAIC = getNimbleOption("MCMCenableWAIC"),
  controlWAIC = list(),
  print = getNimbleOption("verbose"),
  autoBlock = FALSE,
  oldConf,
  ...
)

Arguments

Details

See MCMCconf for details on how to manipulate the MCMCconf object

Author(s)

Daniel Turek

See Also

buildMCMC runMCMC nimbleMCMC

Description

Modifies an MCMC configuration object to perform a reversible jump MCMC sampling for variable selection, using a univariate normal proposal distribution. Users can control the mean and scale of the proposal. This function supports two different types of model specification: with and without indicator variables.

Usage

configureRJ(
  conf,
  targetNodes,
  indicatorNodes = NULL,
  priorProb = NULL,
  control = list(mean = NULL, scale = NULL, fixedValue = NULL)
)

Arguments

Details

This function modifies the samplers in MCMC configuration object for each of the nodes provided in the targetNodes argument. To these elements two samplers are assigned: a reversible jump sampler to transition the variable in/out of the model, and a modified version of the original sampler, which performs updates only when the target node is already in the model.

configureRJ can handle two different ways of writing a NIMBLE model, either with or without indicator variables. When using indicator variables, the indicatorNodes argument must be provided. Without indicator variables, the priorProb argument must be provided. In the latter case, the user can provide a non-zero value for fixedValue if desired.

Note that this functionality is intended for variable selection in regression-style models but may be useful for other situations as well. At the moment, setting a variance component to zero and thereby removing a set of random effects that are explicitly part of a model will not work because MCMC sampling in that case would need to propose values for multiple parameters (the random effects), whereas the current functionality only proposes adding/removing a single model node.

Value

NULL configureRJ modifies the input MCMC configuration object in place.

Author(s)

Sally Paganin, Perry de Valpine, Daniel Turek

References

Peter J. Green. (1995). Reversible jump Markov chain Monte Carlo computation and Bayesian model determination. Biometrika, 82(4), 711-732.

See Also

samplers configureMCMC

Examples

## Not run: 

## Linear regression with intercept and two covariates, using indicator variables

code <- nimbleCode({
  beta0 ~ dnorm(0, sd = 100)
  beta1 ~ dnorm(0, sd = 100)
  beta2 ~ dnorm(0, sd = 100)
  sigma ~ dunif(0, 100) 
  z1 ~ dbern(psi)   ## indicator variable associated with beta1
  z2 ~ dbern(psi)   ## indicator variable associated with beta2
  psi ~ dunif(0, 1) ## hyperprior on inclusion probability
  for(i in 1:N) {
    Ypred[i] <- beta0 + beta1 * z1 * x1[i] + beta2 * z2 * x2[i]
    Y[i] ~ dnorm(Ypred[i], sd = sigma)
  }
})

## simulate some data
set.seed(1)
N <- 100
x1 <- runif(N, -1, 1)
x2 <- runif(N, -1, 1) ## this covariate is not included
Y <- rnorm(N, 1 + 2.5 * x1, sd = 1)

## build the model
rIndicatorModel <- nimbleModel(code, constants = list(N = N),
                               data = list(Y = Y, x1 = x1, x2 = x2), 
                               inits = list(beta0 = 0, beta1 = 0, beta2 = 0, sigma = sd(Y),
                               z1 = 1, z2 = 1, psi = 0.5))

indicatorModelConf <- configureMCMC(rIndicatorModel)

## Add reversible jump  
configureRJ(conf = indicatorModelConf,        ## model configuration
            targetNodes = c("beta1", "beta2"), ## coefficients for selection
            indicatorNodes = c("z1", "z2"),    ## indicators paired with coefficients
            control = list(mean = 0, scale = 2))

indicatorModelConf$addMonitors("beta1", "beta2", "z1", "z2")

rIndicatorMCMC <- buildMCMC(indicatorModelConf)
cIndicatorModel <- compileNimble(rIndicatorModel)
cIndicatorMCMC <- compileNimble(rIndicatorMCMC, project = rIndicatorModel)

set.seed(1)
samples <- runMCMC(cIndicatorMCMC, 10000, nburnin = 6000)

## posterior probability to be included in the mode
mean(samples[ , "z1"])
mean(samples[ , "z2"])

## posterior means when in the model
mean(samples[ , "beta1"][samples[ , "z1"] != 0])
mean(samples[ , "beta2"][samples[ , "z2"] != 0])


## Linear regression with intercept and two covariates, without indicator variables

code <- nimbleCode({
  beta0 ~ dnorm(0, sd = 100)
  beta1 ~ dnorm(0, sd = 100)
  beta2 ~ dnorm(0, sd = 100)
  sigma ~ dunif(0, 100)
  for(i in 1:N) {
    Ypred[i] <- beta0 + beta1 * x1[i] + beta2 * x2[i]
    Y[i] ~ dnorm(Ypred[i], sd = sigma)
  }
})

rNoIndicatorModel <- nimbleModel(code, constants = list(N = N),
                                 data = list(Y = Y, x1 = x1, x2 = x2), 
                                 inits=  list(beta0 = 0, beta1 = 0, beta2 = 0, sigma = sd(Y)))

noIndicatorModelConf <- configureMCMC(rNoIndicatorModel)

## Add reversible jump  
configureRJ(conf = noIndicatorModelConf,      ## model configuration
            targetNodes = c("beta1", "beta2"), ## coefficients for selection   
            priorProb = 0.5,                   ## prior probability of inclusion
            control = list(mean = 0, scale = 2))

## add monitors
noIndicatorModelConf$addMonitors("beta1", "beta2")
rNoIndicatorMCMC <- buildMCMC(noIndicatorModelConf) 

cNoIndicatorModel <- compileNimble(rNoIndicatorModel)
cNoIndicatorMCMC <- compileNimble(rNoIndicatorMCMC, project = rNoIndicatorModel)

set.seed(1)
samples <- runMCMC(cNoIndicatorMCMC, 10000, nburnin = 6000)

## posterior probability to be included in the mode
mean(samples[ , "beta1"] != 0)
mean(samples[ , "beta2"] != 0)

## posterior means when in the model
mean(samples[ , "beta1"][samples[ , "beta1"] != 0])
mean(samples[ , "beta2"][samples[ , "beta2"] != 0])

## End(Not run)

Description

Calculations to handle censoring

Usage

dconstraint(x, cond, log = FALSE)

rconstraint(n = 1, cond)

Arguments

Details

Used for working with constraints in BUGS code. See the NIMBLE manual for additional details.

Value

dconstraint gives the density and rconstraint generates random deviates, but these are unusual as the density is 1 if x matches cond and 0 otherwise and the deviates are simply the value of cond

Author(s)

Christopher Paciorek

See Also

Distributions for other standard distributions

Examples

constr <- 3 > 2 && 4 > 0
x <- rconstraint(1, constr)
dconstraint(x, constr)
dconstraint(0, 3 > 4)
dconstraint(1, 3 > 4)
rconstraint(1, 3 > 4)

Description

This function returns a logical TRUE/FALSE value, indicating whether the proposed transition should be accepted (TRUE) or rejected (FALSE).

Usage

decide(logMetropolisRatio)

Arguments

Details

The Metropolis-Hastings accept/reject decisions is made as follows. If logMetropolisRatio is greater than 0, accept (return TRUE). Otherwise draw a uniform random number between 0 and 1 and accept if it is less that exp(logMetropolisRatio. The proposed transition will be rejected (return FALSE). If logMetropolisRatio is NA, NaN, or -Inf, a reject (FALSE) decision will be returned.

Author(s)

Daniel Turek

Description

This nimbleFunction generator must be specialized to three required arguments: a model, a modelValues, and a character vector of node names.

Usage

decideAndJump(model, mvSaved, target, UNUSED)

Arguments

Details

Calling decideAndJump(model, mvSaved, target) will generate a specialized nimbleFunction with four required numeric arguments:

modelLP1: The model log-probability associated with the newly proposed value(s)

modelLP0: The model log-probability associated with the original value(s)

propLP1: The log-probability associated with the proposal forward-transition

propLP0: The log-probability associated with the proposal reverse-tranisiton

Executing this function has the following effects: – Calculate the (log) Metropolis-Hastings ratio, as logMHR = modelLP1 - modelLP0 - propLP1 + propLP0 – Make the proposal acceptance decision based upon the (log) Metropolis-Hastings ratio – If the proposal is accepted, the values and associated logProbs of all calcNodes are copied from the model object into the mvSaved object – If the proposal is rejected, the values and associated logProbs of all calcNodes are copied from the mvSaved object into the model object – Return a logical value, indicating whether the proposal was accepted

Author(s)

Daniel Turek

Description

Explicitly declare a variable in run-time code of a nimbleFunction, for cases when its dimensions cannot be inferred before it is used. Works in R and NIMBLE.

Usage

declare(name, def)

Arguments

Details

In a run-time function of a nimbleFunction (either the run function or a function provided in methods when calling nimbleFunction), the dimensionality and numeric type of a variable is inferred when possible from the statement first assigning into it. E.g. A <- B + C infers that A has numeric types, dimensions and sizes taken from B + C. However, if the first appearance of A is e.g. A[i] <- 5, A must have been explicitly declared. In this case, declare(A, double(1)) would make A a 1-dimensional (i.e. vector) double.

When sizes are not set, they can be set by a call to setSize or by assignment to the whole object. Sizes are not automatically extended if assignment is made to elements beyond the current sizes. In compiled nimbleFunctions doing so can cause a segfault and crash the R session.

This part of the NIMBLE language is needed for compilation, but it also runs in R. When run in R, is works by the side effect of creating or modifying name in the calling environment.

Author(s)

NIMBLE development team

Examples

declare(A, logical())             ## scalar logical, the only kind allowed
declare(B, integer(2, c(10, 10))) ## 10 x 10 integer matrix
declare(C, double(3))             ## 3-dimensional double array with no sizes set.

Description

Deregister distributional information originally supplied by the user for use in BUGS model code

Usage

deregisterDistributions(distributionsNames, userEnv = parent.frame())

Arguments

Author(s)

Christopher Paciorek

Description

Density and random generation for the Dirichlet distribution

Usage

ddirch(x, alpha, log = FALSE)

rdirch(n = 1, alpha)

Arguments

Details

See Gelman et al., Appendix A or the BUGS manual for mathematical details.

Value

ddirch gives the density and rdirch generates random deviates.

Author(s)

Christopher Paciorek

References

Gelman, A., Carlin, J.B., Stern, H.S., and Rubin, D.B. (2004) Bayesian Data Analysis, 2nd ed. Chapman and Hall/CRC.

See Also

Distributions for other standard distributions

Examples

alpha <- c(1, 10, 30)
x <- rdirch(1, alpha)
ddirch(x, alpha)

Description

Give information about each BUGS distribution

Usage

getDistributionInfo(dist)

isUserDefined(dist)

pqDefined(dist)

getType(
  dist,
  params = NULL,
  valueOnly = is.null(params) && !includeParams,
  includeParams = !is.null(params)
)

getParamNames(dist, includeValue = TRUE)

Arguments

Details

NIMBLE provides various functions to give information about a BUGS distribution. In some cases, functions of the same name and similar functionality operate on the node(s) of a model as well (see help(modelBaseClass)).

getDistributionInfo returns an internal data structure (a reference class object) providing various information about the distribution. The output is not very user-friendly, but does contain all of the information that NIMBLE has about the distribution.

isDiscrete tests if a BUGS distribution is a discrete distribution.

isUserDefined tests if a BUGS distribution is a user-defined distribution.

pqAvail tests if a BUGS distribution provides distribution ('p') and quantile ('q') functions.

getDimension provides the dimension of the value and/or parameters of a BUGS distribution. The return value is a numeric vector with an element for each parameter/value requested.

getType provides the type (numeric, logical, integer) of the value and/or parameters of a BUGS distribution. The return value is a character vector with an element for each parameter/value requested. At present, all quantities are stored as numeric (double) values, so this function is of little practical use but could be exploited in the future.

getParamNames provides the value and/or parameter names of a BUGS distribution.

Author(s)

Christopher Paciorek

Examples

distInfo <- getDistributionInfo('dnorm')
distInfo
distInfo$range

isDiscrete('dbin')

isUserDefined('dbin')

pqDefined('dgamma')
pqDefined('dmnorm')

getDimension('dnorm')
getDimension('dnorm', includeParams = TRUE)
getDimension('dnorm', c('var', 'sd'))
getDimension('dcat', includeParams = TRUE)
getDimension('dwish', includeParams = TRUE)

getType('dnorm')
getType('dnorm', includeParams = TRUE)
getType('dnorm', c('var', 'sd'))
getType('dcat', includeParams = TRUE)
getType('dwish', includeParams = TRUE)

getParamNames('dnorm', includeValue = FALSE)
getParamNames('dmnorm')

Description

Density, distribution function, quantile function and random generation for the double exponential distribution, allowing non-zero location, mu, and non-unit scale, sigma, or non-unit rate, tau

Usage

ddexp(x, location = 0, scale = 1, rate = 1/scale, log = FALSE)

rdexp(n, location = 0, scale = 1, rate = 1/scale)

pdexp(
  q,
  location = 0,
  scale = 1,
  rate = 1/scale,
  lower.tail = TRUE,
  log.p = FALSE
)

qdexp(
  p,
  location = 0,
  scale = 1,
  rate = 1/scale,
  lower.tail = TRUE,
  log.p = FALSE
)

Arguments

Details

See Gelman et al., Appendix A or the BUGS manual for mathematical details.

Value

ddexp gives the density, pdexp gives the distribution function, qdexp gives the quantile function, and rdexp generates random deviates.

Author(s)

Christopher Paciorek

References

Gelman, A., Carlin, J.B., Stern, H.S., and Rubin, D.B. (2004) Bayesian Data Analysis, 2nd ed. Chapman and Hall/CRC.

See Also

Distributions for other standard distributions

Examples

x <- rdexp(50, location = 2, scale = 1)
ddexp(x, 2, 1)

Description

nimbleList definition for the type of nimbleList returned by nimEigen.

Usage

eigenNimbleList

Format

An object of class list of length 1.

Author(s)

NIMBLE development team

See Also

nimEigen

Description

Density, distribution function, quantile function and random generation for the exponential distribution with rate (i.e., mean of 1/rate) or scale parameterizations.

Usage

dexp_nimble(x, rate = 1/scale, scale = 1, log = FALSE)

rexp_nimble(n = 1, rate = 1/scale, scale = 1)

pexp_nimble(q, rate = 1/scale, scale = 1, lower.tail = TRUE, log.p = FALSE)

qexp_nimble(p, rate = 1/scale, scale = 1, lower.tail = TRUE, log.p = FALSE)

Arguments

Details

NIMBLE's exponential distribution functions use Rmath's functions under the hood, but are parameterized to take both rate and scale and to use 'rate' as the core parameterization in C, unlike Rmath, which uses 'scale'. See Gelman et al., Appendix A or the BUGS manual for mathematical details.

Value

dexp_nimble gives the density, pexp_nimble gives the distribution function, qexp_nimble gives the quantile function, and rexp_nimble generates random deviates.

Author(s)

Christopher Paciorek

References

Gelman, A., Carlin, J.B., Stern, H.S., and Rubin, D.B. (2004) Bayesian Data Analysis, 2nd ed. Chapman and Hall/CRC.

See Also

Distributions for other standard distributions

Examples

x <- rexp_nimble(50, scale = 3)
dexp_nimble(x, scale = 3)

Description

Extract named elements from MCMC sampler control list

Usage

extractControlElement(controlList, elementName, defaultValue, error)

Arguments

Value

The element of controlList whose name matches elementName. If no controlList name matches elementName, then defaultValue is returned.

Author(s)

Daniel Turek

Description

Improper flat distribution for use as a prior distribution in BUGS models

Usage

dflat(x, log = FALSE)

rflat(n = 1)

dhalfflat(x, log = FALSE)

rhalfflat(n = 1)

Arguments

Value

dflat gives the pseudo-density value of 1, while rflat and rhalfflat return NaN, since one cannot simulate from an improper distribution. Similarly, dhalfflat gives a pseudo-density value of 1 when x is non-negative.

Author(s)

Christopher Paciorek

See Also

Distributions for other standard distributions

Examples

dflat(1)

Description

Get the value of the lower or upper bound for a single stochastic node in a model.

Usage

getBound(model, node, bound, nodeFunctionIndex)

Arguments

Details

Standard usage is as a method of a model, in the form model$getBound(node, bound), but the usage as a simple function with the model as the first argument as above is also allowed.

For nodes that do not involve truncation of the distribution this will return the lower or upper bound of the distribution, which may be a constant or for a limited number of distributions a parameter or functional of a parameter (at the moment in NIMBLE, the only case where a bound is a parameter is for the uniform distribution. For nodes that are truncated, this will return the desired bound, which may be a functional of other quantities in the model or may be a constant.

Description

NIMBLE comes with some of the classic BUGS examples. getBUGSexampleDir looks up the location of an example from its name.

Usage

getBUGSexampleDir(example)

Arguments

Value

Character string of the fully pathed directory of the BUGS example.

Author(s)

Christopher Paciorek

See Also

readBUGSmodel for usage in creating a model from a classic BUGS example

Description

A conditionally independent set of nodes is such that the joint probability (density) of nodes in the set will not change even if any non-given node outside the set is changed. Default given nodes are data nodes and parameter nodes (aka "top-level" nodes, i.e. nodes with no parent nodes), but this can be controlled.

Usage

getConditionallyIndependentSets(
  model,
  nodes,
  givenNodes,
  omit = integer(),
  explore = c("both", "down", "up"),
  unknownAsGiven = TRUE,
  returnType = "names",
  returnScalarComponents = FALSE,
  endAsGiven = FALSE
)

Arguments

Details

This function returns sets of conditionally independent nodes. Multiple input nodes might be in the same set or different sets.

The nodes input and the returned sets include only stochastic nodes because conditional independence is a property of random variables. Deterministic nodes are considered in determining the sets. givenNodes may contain stochastic or deterministic nodes.

Value

List of nodes that are in conditionally independent sets. With each set, nodes are returned in topologically sorted order. The sets themselves are returned in topologically sorted order of their first nodes.

Other nodes (not in nodes) may be included in the output if unknownAsGiven=FALSE.

Author(s)

Perry de Valpine

Description

Returns a list containing the nimbleFunction definition components (setup function, run function, and other member methods) for the supplied nimbleFunction generator or specialized instance.

Usage

getDefinition(nf)

Arguments

Author(s)

Daniel Turek

Description

Get a list of all parameter names (or certain categories of parameters) generated by model macros in the model code.

Usage

getMacroParameters(
  model,
  includeLHS = TRUE,
  includeRHS = TRUE,
  includeDeterm = TRUE,
  includeStoch = TRUE,
  includeIndexPars = FALSE
)

Arguments

Details

Some model macros will generate new parameters to be included in the output code. NIMBLE automatically detects these new parameters and records them in the model object. This function allows easy access to this stored list, or subsets of it (see arguments).

Value

A named list of generated parameters, with the element names corresponding to the original source macro.

Examples

nimbleOptions(enableModelMacros = TRUE)
nimbleOptions(enableMacroComments = FALSE)
nimbleOptions(verbose = FALSE)

testMacro <- list(process = function(code, modelInfo, .env){
  code <- quote({
    for (i_ in 1:n){
      mu[i_] <- alpha + beta
      y[i_] ~ dnorm(0, sigma)
    }
    alpha ~ dnorm(0, 1)
  })
  list(code = code, modelInfo=modelInfo)
})
 class(testMacro) <- "model_macro"

code <- nimbleCode({
  y[1:n] ~ testMacro()
})

const <- list(y = rnorm(10), n = 10)

mod <- nimbleModel(code, constants=const)

mod$getMacroParameters()
# should be list(testMacro = list(c("mu", "alpha", "beta", "sigma")))

mod$getMacroParameters(includeRHS = FALSE)
# should be list(testMacro = list(c("mu", "alpha")))

Description

Allow the user to get the value of a global _option_ that affects the way in which NIMBLE operates

Usage

getNimbleOption(x)

Arguments

Value

The value of the option.

Author(s)

Christopher Paciorek

Examples

getNimbleOption('verifyConjugatePosteriors')

Description

Get the value of a parameter for any single stochastic node in a model.

Usage

getParam(model, node, param, nodeFunctionIndex, warn = TRUE)

Arguments

Details

Standard usage is as a method of a model, in the form model$getParam(node, param), but the usage as a simple function with the model as the first argument as above is also allowed.

For example, suppose node 'x[1:5]' follows a multivariate normal distribution (dmnorm) in a model declared by BUGS code. model$getParam('x[1:5]', 'mean') would return the current value of the mean parameter (which may be determined from other nodes). The parameter requested does not have to be part of the parameterization used to declare the node. Rather, it can be any parameter known to the distribution. For example, one can request the scale or rate parameter of a gamma distribution, regardless of which one was used to declare the node.

Description

This function obtains posterior samples from a Dirichlet process distributed random measure of a model specified using the dCRP distribution.

Usage

getSamplesDPmeasure(
  MCMC,
  epsilon = 1e-04,
  setSeed = FALSE,
  progressBar = getNimbleOption("MCMCprogressBar")
)

Arguments

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 dCRP distribution.

The MCMC argument is an object of class MCMC provided by buildMCMC, or its compiled version. The MCMC should already have been run, as 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 dCRP distribution), the cluster parameter variables, all variables directly determining the dCRP concentration parameter, and any stochastic parent variables of the cluster parameter variables. See help(configureMCMC) or help(addMonitors) for information on specifying monitors for an MCMC.

The epsilon argument is optional and used to determine the truncation level of the random measure. 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 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.

Author(s)

Claudia Wehrhahn and Christopher Paciorek

References

Sethuraman, J. (1994). A constructive definition of Dirichlet priors. Statistica Sinica, 639-650.

Gelfand, A.E. and Kottas, A. (2002). A computational approach for full nonparametric Bayesian inference under Dirichlet process mixture models. Journal of Computational and Graphical Statistics, 11(2), 289-305.

See Also

buildMCMC, configureMCMC,

Examples

## Not run: 
  conf <- configureMCMC(model)
  mcmc <- buildMCMC(conf)
  cmodel <- compileNimble(model)
  cmcmc <- compileNimble(mcmc, project = model)
  runMCMC(cmcmc, niter = 1000)
  outputG <- getSamplesDPmeasure(cmcmc)

## End(Not run)

Description

Returns the number of rows of NIMBLE modelValues object. Works in R and NIMBLE.

Usage

getsize(container)

Arguments

Details

See the User Manual or help(modelValuesBaseClass) for information about modelValues objects

Author(s)

Clifford Anderson-Bergman

Examples

mvConf <- modelValuesConf(vars = 'a', types = 'double', sizes = list(a = 1) )
mv <- modelValues(mvConf)
 resize(mv, 10)
getsize(mv)

Description

Returns a d-by-d identity matrix (square matrix of 0's, with 1's on the main diagnol).

Usage

identityMatrix(d)

Arguments

Details

This function can be used in NIMBLE run code. It is deprecated because now one can use diag(d) instead.

Value

A d-by-d identity matrix

Author(s)

Daniel Turek

Examples

Id <- identityMatrix(d = 3)

Description

Performs initialization of nimble model node values and log probabilities

Usage

initializeModel(model, silent = FALSE)

Arguments

Details

This nimbleFunction may be used at the beginning of nimble algorithms to perform model initialization. The intended usage is to specialize an instance of this nimbleFunction in the setup function of an algorithm, then execute that specialied function at the beginning of the algorithm run function. The specialized function takes no arguments.

Executing this function ensures that all right-hand-side only nodes have been assigned real values, that all stochastic nodes have a real value, or otherwise have their simulate() method called, that all deterministic nodes have their simulate() method called, and that all log-probabilities have been calculated with the current model values. An error results if model initialization encounters a problem, for example a missing right-hand-side only node value.

Author(s)

Daniel Turek

Examples

myNewAlgorithm <- nimbleFunction(
   setup = function(model, ...) {
      my_initializeModel <- initializeModel(model)
      ....
   },
   run = function(...) {
      my_initializeModel()
      ....
   }
)

Description

Calculations to handle censoring

Usage

dinterval(x, t, c, log = FALSE)

rinterval(n = 1, t, c)

Arguments

Details

Used for working with censoring in BUGS code. Taking c to define the endpoints of two or more intervals (with implicit endpoints of plus/minus infinity), x (or the return value of rinterval) gives the non-negative integer valued index of the interval in which t falls. See the NIMBLE manual for additional details.

Value

dinterval gives the density and rinterval generates random deviates, but these are unusual as the density is 1 if x indicates the interval in which t falls and 0 otherwise and the deviates are simply the interval(s) in which t falls.

Author(s)

Christopher Paciorek

See Also

Distributions for other standard distributions

Examples

endpoints <- c(-3, 0, 3)
vals <- c(-4, -1, 1, 5)
x <- rinterval(4, vals, endpoints)
dinterval(x, vals, endpoints)
dinterval(c(1, 5, 2, 3), vals, endpoints)

Description

Density, distribution function, quantile function and random generation for the inverse gamma distribution with rate or scale (mean = scale / (shape - 1)) parameterizations.

Usage

dinvgamma(x, shape, scale = 1, rate = 1/scale, log = FALSE)

rinvgamma(n = 1, shape, scale = 1, rate = 1/scale)

pinvgamma(
  q,
  shape,
  scale = 1,
  rate = 1/scale,
  lower.tail = TRUE,
  log.p = FALSE
)

qinvgamma(
  p,
  shape,
  scale = 1,
  rate = 1/scale,
  lower.tail = TRUE,
  log.p = FALSE
)

Arguments

Details

The inverse gamma distribution with parameters shape $=\alpha$ and scale $=\sigma$ has density

$f(x)= \frac{s^a}{\Gamma(\alpha)} {x}^{-(\alpha+1)} e^{-\sigma/x}$

for $x \ge 0$, $\alpha > 0$ and $\sigma > 0$. (Here $\Gamma(\alpha)$ is the function implemented by R's gamma() and defined in its help.

The mean and variance are $E(X) = \frac{\sigma}{\alpha}-1$ and $Var(X) = \frac{\sigma^2}{(\alpha-1)^2 (\alpha-2)}$, with the mean defined only for $\alpha > 1$ and the variance only for $\alpha > 2$.

See Gelman et al., Appendix A or the BUGS manual for mathematical details.

Value

dinvgamma gives the density, pinvgamma gives the distribution function, qinvgamma gives the quantile function, and rinvgamma generates random deviates.

Author(s)

Christopher Paciorek

References

Gelman, A., Carlin, J.B., Stern, H.S., and Rubin, D.B. (2004) Bayesian Data Analysis, 2nd ed. Chapman and Hall/CRC.

See Also

Distributions for other standard distributions

Examples

x <- rinvgamma(50, shape = 1, scale = 3)
dinvgamma(x, shape = 1, scale = 3)

Description

Density and random generation for the Inverse Wishart distribution, using the Cholesky factor of either the scale matrix or the rate matrix.

Usage

dinvwish_chol(x, cholesky, df, scale_param = TRUE, log = FALSE)

rinvwish_chol(n = 1, cholesky, df, scale_param = TRUE)

Arguments

Details

See Gelman et al., Appendix A for mathematical details. The rate matrix as used here is defined as the inverse of the scale matrix, $S^{-1}$, given in Gelman et al.

Value

dinvwish_chol gives the density and rinvwish_chol generates random deviates.

Author(s)

Christopher Paciorek

References

Gelman, A., Carlin, J.B., Stern, H.S., and Rubin, D.B. (2004) Bayesian Data Analysis, 2nd ed. Chapman and Hall/CRC.

See Also

Distributions for other standard distributions

Examples

df <- 40
ch <- chol(matrix(c(1, .7, .7, 1), 2))
x <- rwish_chol(1, ch, df = df)
dwish_chol(x, ch, df = df)

Description

Checks an object to determine if it is a nimbleFunction (i.e., a function created by nimbleFunction using setup code).

Usage

is.nf(f, inputIsName = FALSE, where = -1)

Arguments

See Also

nimbleFunction for how to create a nimbleFunction

Description

Checks an object to determine if it is a nimbleList (i.e., a list created by nlDef$new()).

Usage

is.nl(l)

Arguments

See Also

nimbleList for how to create a nimbleList

Description

Density and random generation for the LKJ distribution for the Cholesky factor of a correlation matrix.

Usage

dlkj_corr_cholesky(x, eta, p, log = FALSE)

rlkj_corr_cholesky(n = 1, eta, p)

Arguments

Details

See Stan Development Team for mathematical details.

Value

dlkj_corr_cholesky gives the density and rlkj_corr_cholesky generates random deviates.

Author(s)

Christopher Paciorek

References

Stan Development Team. Stan Reference Functions, version 2.27.

See Also

Distributions for other standard distributions

Examples

eta <- 3
x <- rlkj_corr_cholesky(1, eta, 5)
dlkj_corr_cholesky(x, eta, 5)

Description

Creates a simple getBound_info object, which has a list with a boundID and a type. Unlike makeParamInfo this is more bare-bones, but keeping it for parallelism with getParam.

Usage

makeBoundInfo(model, nodes, bound)

Arguments

Details

This is used internally by getBound. It is not intended for direct use by a user or even a nimbleFunction programmer.

Description

Inspect structure of a nimble model to determine nodes needed as "update" and/or "constant" entries in usage of nimDerivs. This will typically be used in the setup code of a nimbleFunction.

Usage

makeModelDerivsInfo(model, wrtNodes, calcNodes, dataAsConstantNodes = TRUE)

Arguments

Details

In the compilable parts of a nimbleFunction (i.e. run or other method code, not setup code), a call like nimDerivs(foo(x), ...) records derivatives of foo(x). If foo contains any calls to model$calculate(calcNodes), it may be necessary to provide auxiliary information about the model in further arguments to nimDerivs, specifically the model, updateNodes and constantNodes arguments. 'makeModelDerivsInfo' is a utility to set up that information for typical use cases. It returns a list with elements updateNodes and constantNodes to be passed as arguments of the same name to nimDerivs (along with passing the model as the model argument).

The reason auxiliary information is needed is that recording of derivatives uses a different model than for regular calculations. Together, updateNodes and constantNodes should contain all nodes whose values are needed for the model calculations being recorded and that are not part of wrtNodes. These may include parents of nodes that are in calcNodes but are not themselves in calcNodes, as well as the values of stochastic nodes in calcNodes, which are needed to calculate the corresponding log probabilities. updateNodes will have their values updated from the regular model every time that recorded derivative calculations are used. constantNodes will not be updated every time, which means their values will be permanently fixed either the first time the call to 'nimDerivs' is invoked or on any subsequent call that has reset=TRUE. Use of constantNodes can be slightly more efficient, but one must be careful to be aware that values will not be updated unless reset=TRUE. See the automatic differentiation section of the User Manual for more information.

In the above explanation, care must be taken to understand what should be included in wrtNodes. In a typical use case, some arguments to foo are put into the model using values(model, nodes) <- some_foo_arguments. Next there is typically a call to model$calculate(calcNodes). Here the nodes are considered "with-respect-to" nodes because derivative tracking will follow the arguments of foo, including when they are put into a model and hence used in model$calculate. Therefore these nodes should be the wrtNodes for makeModelDerivsInfo.

Value

A list with elements updateNodes and constantNodes. These shouls be provided as the same-named arguments to nimDerivs (same as derivs).

When using double-taping of derivatives (i.e. foo contains another call to nimDerivs), both calls to nimDerivs should include the model, updateNodes, and constantNodes arguments.

Description

Creates a simple getParam_info object, which has a list with a paramID and a type

Usage

makeParamInfo(model, nodes, param, vector = FALSE)

Arguments

Details

This is used internally by getParam. It is not intended for direct use by a user or even a nimbleFunction programmer.

Description

Objects of this class configure an MCMC algorithm, specific to a particular model. Objects are normally created by calling configureMCMC. Given an MCMCconf object, the actual MCMC function can be built by calling buildMCMC(conf). See documentation below for method initialize() for details of creating an MCMCconf object.

Methods

addDefaultSampler( nodes = character(), control = list(), useConjugacy = getNimbleOption("MCMCuseConjugacy"), onlyRW = FALSE, onlySlice = FALSE, multivariateNodesAsScalars = getNimbleOption("MCMCmultivariateNodesAsScalars"), print = TRUE, ... )

For internal use. Adds default MCMC samplers to the specified nodes.

addMonitors(..., ind = 1, print = TRUE)

Adds variables to the list of monitors.

Arguments:

...: One or more character vectors of indexed nodes, or variables, which are to be monitored. These are added onto the current monitors list.

print: A logical argument specifying whether to print all current monitors (default TRUE).

Details:

See the initialize() function

addMonitors2(..., print = TRUE)

Adds variables to the list of monitors2.

Arguments:

...: One or more character vectors of indexed nodes, or variables, which are to be monitored. These are added onto the current monitors2 list.

print: A logical argument specifying whether to print all current monitors (default TRUE).

Details:

See the initialize() function

addOneSampler( thisSamplerName, samplerFunction, targetOne, thisControlList, allowData, print )

For internal use only

addSampler( target = character(), type = "RW", control = list(), print = NULL, name, targetByNode = FALSE, multivariateNodesAsScalars = getNimbleOption("MCMCmultivariateNodesAsScalars"), expandTarget, scalarComponents, silent = FALSE, default = FALSE, useConjugacy = getNimbleOption("MCMCuseConjugacy"), onlyRW = FALSE, onlySlice = FALSE, allowData = FALSE, ... )

Adds a sampler to the list of samplers contained in the MCMCconf object.

Arguments:

target: The target node or nodes to be sampled. This may be specified as a character vector of model node and/or variable names. For univariate samplers, only a single target node should be provided (unless 'targetByNode' is TRUE). For multivariate samplers, one instance of the multivariate sampler will be assigned to all nodes specified. Nodes are specified in combination with the 'targetByNode' and 'multivariateNodesAsScalars' arguments.

type: When 'default' is FALSE, specifies the type of sampler to add, specified as either a character string or a nimbleFunction object. If the character argument type='newSamplerType', then either newSamplerType or sampler_newSamplerType must correspond to a nimbleFunction (i.e. a function returned by nimbleFunction, not a specialized 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 sampler. The default value is 'RW' which specifies scalar adaptive Metropolis-Hastings sampling with a normal proposal distribution. This default will result in an error if 'target' specifies more than one target node (unless 'targetByNode' is TRUE). This argument is not used when the 'default' argument is TRUE.

control: An optional list of control arguments to sampler functions. These will override those specified in the control list argument to configureMCMC. If a control list is provided, the elements will be provided to all sampler functions which utilize the named elements given. For example, the standard Metropolis-Hastings random walk sampler (sampler_RW) utilizes control list elements 'adaptive', 'adaptInterval', 'scale'. The default values for control list arguments for samplers (if not otherwise provided as an argument to configureMCMC or addSampler) are contained in the setup code of each sampling algorithm.

print: Logical argument, specifying whether to print the details of newly added sampler(s).

name: Optional character string name for the sampler, which is used by the printSamplers method. If 'name' is not provided, the 'type' argument is used to generate the sampler name.

targetByNode: Logical argument, with default FALSE. This arguments controls whether separate instances of the specified sampler 'type' should be assigned to each node contained in 'target'. When FALSE, a single instance of sampler 'type' is assigned to operate on 'target'. When TRUE, potentially multiple instances of sampler 'type' will be added to the MCMC configuration, operating on the distinct nodes which compose 'target'. For example, if 'target' is a vector of distinct node names, then a separate sampler will be assigned to each node in this vector. If 'target' is a model variable which itself is comprised of multiple distinct nodes, then a separate sampler is assigned to each node composing the 'target' variable. Additional control of the handling of multivariate nodes is provided using the 'multivariateNodesAsScalars' argument.

multivariateNodesAsScalars: Logical argument, with default value FALSE. This argument is used in two ways. Functionally, both uses result in separate instances of samplers being added to the scalar components which compose multivariate nodes. See details below.

silent: Logical argument, specifying whether to print warning messages when assigning samplers.

default: Logical argument, with default value FALSE. When FALSE, the 'type' argument dictates what sampling algorithm is assigned to the specified nodes. When TRUE, default samplers will be assigned to the specified nodes following the same logic as the configureMCMC method, and also using the 'useConjugacy', 'onlyRW', 'onlySlice' and 'multivariateNodesAsScalars' arguments.

useConjugacy: Logical argument, with default value TRUE. If specified as FALSE, then no conjugate samplers will be used, even when a node is determined to be in a conjugate relationship. This argument is only used when the 'default' argument is TRUE.

onlyRW: Logical argument, with default value FALSE. If specified as TRUE, then Metropolis-Hastings random walk samplers will be assigned for all non-terminal continuous-valued nodes nodes. Discrete-valued nodes are assigned a slice sampler, and terminal nodes are assigned a posterior_predictive sampler. This argument is only used when the 'default' argument is TRUE.

onlySlice: 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. This argument is only used when the 'default' argument is TRUE.

allowData: Logical argument, with default value FALSE. When FALSE, samplers will not be assigned to operate on data nodes, even if data nodes are included in 'target'. When TRUE, samplers will be assigned to 'target' without regard to whether nodes are designated as data.

...: Additional named arguments passed through ... will be used as additional control list elements.

Details:

Samplers are added to the end of the list of samplers for this MCMCconf object, and do not replace any existing samplers. Samplers are removed using the removeSamplers method.

Invisibly returns a list of the current sampler configurations, which are samplerConf reference class objects.

'multivariateNodesAsScalars' has two usages. The first usage occurs when 'targetByNode' is TRUE and therefore separate instances of sampler 'type' are assigned to each node which compose 'target'. In this first usage, this argument controls how multivariate nodes (those included in the 'target') are handled. If FALSE, any multivariate nodes in 'target' have a single instance of sampler 'type' assigned. If TRUE, any multivariate nodes appearing in 'target' are themselves decomposed into their scalar elements, and a separate instance of sampler 'type' is assigned to operate on each scalar element.

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'.

getMonitors()

Returns a character vector of the current monitors

Details:

See the initialize() function

getMonitors2()

Returns a character vector of the current monitors2

Details:

See the initialize() function

getSamplerDefinition(ind, print = FALSE)

Returns the nimbleFunction definition of an MCMC sampler.

Arguments:

ind: A numeric vector or character vector. A numeric vector may be used to specify the index of the sampler definition to return, or a character vector may be used to indicate a target node for which the sampler acting on this nodes will be printed. For example, getSamplerDefinition('x[2]') will return the definition of the sampler whose target is model node 'x[2]'. If more than one sampler 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 nimbleFunction sampler.

getSamplerExecutionOrder()

Returns a numeric vector, specifying the ordering of sampler function execution.

The indices of execution specified in this numeric vector correspond to the enumeration of samplers printed by printSamplers(), or returned by getSamplers().

getSamplers(ind)

Returns a list of samplerConf objects.

Arguments:

ind: A numeric vector or character vector. A numeric vector may be used to specify the indices of the samplerConf objects to return, 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 returned. For example, getSamplers('x') will return all samplerConf objects whose target is model node 'x', or whose targets are contained (entirely or in part) in the model variable 'x'. If omitted, then all samplerConf objects in this MCMC configuration object are returned.

initialize( model, nodes, control = list(), monitors, thin = 1, monitors2 = character(), thin2 = 1, useConjugacy = getNimbleOption("MCMCuseConjugacy"), onlyRW = FALSE, onlySlice = FALSE, multivariateNodesAsScalars = getNimbleOption("MCMCmultivariateNodesAsScalars"), enableWAIC = getNimbleOption("MCMCenableWAIC"), controlWAIC = list(), print = TRUE, ... )

Creates a MCMC configuration for a given model. The resulting object is suitable as an argument to buildMCMC.

Arguments:

model: A NIMBLE model object, created from nimbleModel(...)

nodes: An optional character vector, specifying the nodes for which samplers should be created. Nodes may be specified in their indexed form, 'y[1, 3]', or nodes specified without indexing will be expanded fully, e.g., 'x' will be expanded to 'x[1]', 'x[2]', etc. If missing, the default value is all non-data stochastic nodes. If NULL, then no samplers are added.

control: An optional list of control arguments to sampler functions. If a control list is provided, the elements will be provided to all sampler functions which utilize the named elements given. For example, the standard Metropolis-Hastings random walk sampler (sampler_RW) utilizes control list elements 'adaptive', 'adaptInterval', 'scale'. The default values for control list arguments for samplers (if not otherwise provided as an argument to configureMCMC() or addSampler()) are contained in the setup code of each sampling algorithm.

monitors: A character vector of node names or variable names, to record during MCMC sampling. This set of monitors will be recorded with thinning interval 'thin', and the samples will be stored into the 'mvSamples' object. The default value is all top-level stochastic nodes of the model – those having no stochastic parent nodes.

monitors2: A character vector of node names or variable names, to record during MCMC sampling. This set of monitors will be recorded with thinning interval 'thin2', and the samples will be stored into the 'mvSamples2' object. The default value is an empty character vector, i.e. no values will be recorded.

thin: The thinning interval for 'monitors'. Default value is one.

thin2: The thinning interval for 'monitors2'. Default value is one.

useConjugacy: A logical argument, with default value TRUE. If specified as FALSE, then no conjugate samplers will be used, even when a node is determined to be in a conjugate relationship.

onlyRW: A logical argument, with default value FALSE. If specified as TRUE, then Metropolis-Hastings random walk samplers will be assigned for all non-terminal continuous-valued nodes nodes. Discrete-valued nodes are assigned a slice sampler, and terminal nodes are assigned a posterior_predictive sampler.

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.

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.

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)'.

print: A logical argument specifying whether to print the montiors and samplers. Default is TRUE.

...: Additional named control list elements for default samplers, or additional arguments to be passed to the autoBlock function when autoBlock = TRUE.

printMonitors()

Prints all current monitors and monitors2

Details:

See the initialize() function

printSamplers( ..., ind, type, displayControlDefaults = FALSE, displayNonScalars = FALSE, displayConjugateDependencies = FALSE, executionOrder = FALSE, byType = FALSE )

Prints details of the MCMC samplers.

Arguments:

...: Character node or variable names, or numeric indices. Numeric indices may be used to specify the indices of the samplers to print, or character strings 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.

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.

displayConjugateDependencies: A logical argument, specifying whether to display the dependency lists of conjugate samplers (default FALSE).

displayNonScalars: A logical argument, specifying whether to display the values of non-scalar control list elements (default FALSE).

executionOrder: A logical argument, specifying whether to print the sampler functions in the (possibly modified) order of execution (default FALSE).

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).

removeSampler(...)

Alias for removeSamplers method

removeSamplers(..., ind, print = FALSE)

Removes one or more samplers from an MCMCconf object.

This function also has the side effect of resetting the sampler execution ordering so as to iterate over the remaining set of samplers, sequentially, executing each sampler once.

Arguments:

...: Character node names or numeric indices. Character node names specify the node names for samplers to remove, or numeric indices can provide the indices of samplers to remove.

ind: A numeric vector or character vector specifying the samplers to remove. A numeric vector may specify the indices of the samplers to be removed. Alternatively, a character vector may be used to specify a set of model nodes and/or variables, and all samplers whose 'target' is among these nodes will be removed. If omitted, then all samplers are removed.

print: A logical argument specifying whether to print the current list of samplers once the removal has been done (default FALSE).

replaceSampler(...)

Alias for replaceSamplers method

replaceSamplers(...)

Replaces one or more samplers from an MCMCconf object with newly specified sampler(s). Operation and arguments are identical to the 'addSampler' method, with the additional side effect of first removing any existing samplers which operate on the specified node(s).

This function also has the side effect of resetting the sampler execution ordering so as to iterate over the remaining set of samplers, sequentially, executing each sampler once.

See 'addSamplers' for a description of the arguments.

This function also has the side effect of resetting the sampler execution ordering so as to iterate over the newly specified set of samplers, sequentially, executing each sampler once.

resetMonitors()

Resets the current monitors and monitors2 lists to nothing.

Details:

See the initialize() function

setMonitors(..., ind = 1, print = TRUE)

Sets new variables to the list of monitors.

Arguments:

...: One or more character vectors of indexed nodes, or variables, which are to be monitored. These replace the current monitors list.

print: A logical argument specifying whether to print all current monitors (default TRUE).

Details:

See the initialize() function

setMonitors2(..., print = TRUE)

Sets new variables to the list of monitors2.

Arguments:

...: One or more character vectors of indexed nodes, or variables, which are to be monitored. These replace the current monitors2 list.

print: A logical argument specifying whether to print all current monitors (default TRUE).

Details:

See the initialize() function

setSampler(...)

Alias for setSamplers method

setSamplerExecutionOrder(order, print = FALSE)

Sets the ordering in which sampler functions will execute.

This allows some samplers to be "turned off", or others to execute multiple times in a single MCMC iteration. The ordering in which samplers execute can also be interleaved.

Arguments:

order: A numeric vector, specifying the ordering in which the sampler functions will execute. The indices of execution specified in this numeric vector correspond to the enumeration of samplers printed by printSamplers(), or returned by getSamplers(). If this argument is omitted, the sampler execution ordering is reset so as to sequentially execute each sampler once.

print: A logical argument specifying whether to print the current list of samplers in the modified order of execution (default FALSE).

setSamplers(..., ind, print = FALSE)

Sets the ordering of the list of MCMC samplers.

This function also has the side effect of resetting the sampler execution ordering so as to iterate over the specified set of samplers, sequentially, executing each sampler once.

Arguments:

...: Chracter strings or numeric indices. Character names may be used to specify the node names for samplers to retain. Numeric indices may be used to specify the indicies for the new list of MCMC samplers, in terms of the current ordered list of samplers.

ind: A numeric vector or character vector. A numeric vector may be used to specify the indicies for the new list of MCMC samplers, in terms of the current ordered list of samplers. For example, if the MCMCconf object currently has 3 samplers, then the ordering may be reversed by calling MCMCconf$setSamplers(3:1), or all samplers may be removed by calling MCMCconf$setSamplers(numeric(0)).

Alternatively, a character vector may be used to specify a set of model nodes and/or variables, and the sampler list will modified to only those samplers acting on these target nodes.

As another alternative, a list of samplerConf objects may be used as the argument, in which case this ordered list of samplerConf objects will define the samplers in this MCMC configuration object, completely over-writing the current list of samplers. No checking is done to ensure the validity of the contents of these samplerConf objects; only that all elements of the list argument are, in fact, samplerConf objects.

print: A logical argument specifying whether to print the new list of samplers (default FALSE).

setThin(thin, print = TRUE, ind = 1)

Sets the value of thin.

Arguments:

thin: The new value for the thinning interval 'thin'.

print: A logical argument specifying whether to print all current monitors (default TRUE).

Details:

See the initialize() function

setThin2(thin2, print = TRUE)

Sets the value of thin2.

Arguments:

thin2: The new value for the thinning interval 'thin2'.

print: A logical argument specifying whether to print all current monitors (default TRUE).

Details:

See the initialize() function

Author(s)

Daniel Turek

See Also

configureMCMC

Examples

code <- nimbleCode({
 mu ~ dnorm(0, 1)
 x ~ dnorm(mu, 1)
})
Rmodel <- nimbleModel(code)
conf <- configureMCMC(Rmodel)
conf$setSamplers(1)
conf$addSampler(target = 'x', type = 'slice', control = list(adaptInterval = 100))
conf$addMonitors('mu')
conf$addMonitors2('x')
conf$setThin(5)
conf$setThin2(10)
conf$printMonitors()
conf$printSamplers()

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. model_macro_builder 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 nimbleModel that arranges arguments for input function. Macros are an experimental feature and are available only after setting nimbleOptions(enableModelMacros = TRUE).

Usage

model_macro_builder(fun, use3pieces = TRUE, unpackArgs = TRUE)

Arguments

Details

The arguments use3pieces and unpackArgs indicate how fun expects to have arguments arranged from an input line of code (processed by nimbleModel).

Consider the defaults use3pieces = TRUE and unpackArgs = TRUE, for a macro called macro1. In this case, the line of model code x ~ macro1(arg1 = z[1:10], arg2 = "hello") will be passed to fun as fun(stoch = TRUE, LHS = x, arg1 = z[1:10], arg2 = "hello").

If use3pieces = TRUE but unpackArgs = FALSE, then the RHS will be passed as is, without unpacking its arguments into separate arguments to fun. In this case, x ~ macro1(arg1 = z[1:10], arg2 = "hello") will be passed to fun as fun(stoch = TRUE, LHS = x, RHS = macro1(arg1 = z[1:10], arg2 = "hello")).

If use3pieces = FALSE and unpackArgs = FALSE, the entire line of code is passed as a single object. In this case, x ~ macro1(arg1 = z[1:10], arg2 = "hello") will be passed to fun as fun(x ~ macro1(arg1 = z[1:10], arg2 = "hello")). It is also possible in this case to pass a macro without using a ~ or <-. For example, the line macro1(arg1 = z[1:10], arg2 = "hello") will be passed to fun as fun(macro1(arg1 = z[1:10], arg2 = "hello")).

If use3pieces = FALSE and unpackArgs = TRUE, it won't make sense to anticipate a declaration using ~ or <-. Instead, arguments from an arbitrary call will be passed as separate arguments. For example, the line macro1(arg1 = z[1:10], arg2 = "hello") will be passed to fun as fun(arg1 = z[1:10], arg2 = "hello").

In addition, the final two arguments of fun must be called modelInfo and .env respectively.

During macro processing, nimbleModel passes a named list to the modelInfo argument of fun containing, among other things, elements called constants and 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 nimbleModel then moves to the data, this new data will not be retained in the final model object and thus will not be usable.

nimbleModel passes the R environment from which nimbleModel was called to the .env argument.

The fun function must return a named list with two elements: code, the replacement code, and modelInfo, the modelInfo list described above. 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 fun correctly. Functions such as substitute and as.name (e.g. as.name('~')), quote, parse and deparse are particularly handy.

Multiple lines of new code should be contained in {} . Extra curly braces are not a problem. See example 2.

Macro expansion is done recursively: One macro can return code that invokes another macro.

Value

A list of class model_macro with one element called process, which contains the macro function suitable for use by nimbleModel.

Examples

nimbleOptions(enableModelMacros = 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 <- model_macro_builder(
    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 <- model_macro_builder(
    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 <- model_macro_builder(
   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