diff --git a/UserManual/src/chapter_AD.Rmd b/UserManual/src/chapter_AD.Rmd index 26620b1c7..9219b9269 100644 --- a/UserManual/src/chapter_AD.Rmd +++ b/UserManual/src/chapter_AD.Rmd @@ -2,7 +2,7 @@ ```{r, echo=FALSE} require(nimble) -``` +``` # Automatic Derivatives {#cha-AD} @@ -14,7 +14,7 @@ if(!require(nimble, warn.conflicts = FALSE, quietly = TRUE)) { require(methods, warn.conflicts = FALSE, quietly = TRUE) # seems to be needed, but why? require(igraph, warn.conflicts = FALSE, quietly = TRUE) # same question eval <- TRUE -``` +``` As of version 1.0.0, NIMBLE can automatically provide numerically accurate derivatives of potentially arbitrary order for most calculations in models and/or nimbleFunctions. This feature enables methods such as Hamiltonian Monte Carlo (HMC, see package `nimbleHMC`), Laplace (and AGHQ) approximation, and fast optimization with methods that use function gradients. @@ -248,7 +248,7 @@ Derivatives can be tracked through function calls, including from a model to a u For this situation, you can use `ADbreak` to stop derivative tracking. For example, if one has a nimbleFunction argument `x`, one would do `x_noDeriv <- ADbreak(x)` and then use `x_noDeriv` in subsequent calculations. This "breaks" derivative tracking in the sense that the derivative of `x_noDeriv` with respect to `x` is 0. However, it is important to realize that `x_noDeriv` is still a derivative-tracking type (unless it is in the `ignore` set), so it can be passed to other functions involved in the derivative-tracked operations. And, on the other hand, it can't be used as a for-loop index or as the start of end of a for-loop extent. As noted above, if the values of `x` or `x_noDeriv` are used in further calculations, their values will be "baked in" to the tape. -Including `"x"` in the `ignore` set won't work in this case because, if `x` is a function argument, its type will then be changed to be a non-derivative type. If a derivative-tracking type is then passed as the value for `x` there will be a mismatch in types. Since variables passed from models always allow derivative-tracking, one should not `ignore` the corresponding function arguments. +Including `"x"` in the `ignore` set won't work in this case because, if `x` is a function argument, its type will then be changed to be a non-derivative type. If a derivative-tracking type is then passed as the value for `x` there will be a mismatch in types. Since variables passed from models always allow derivative-tracking, one should not `ignore` the corresponding function arguments. Note that `ADbreak` only works for scalars. If you need to break AD tracking for a non-scalar, you have to write a for-loop to apply `ADbreak` for each element. For example: @@ -282,10 +282,10 @@ derivs_demo3 <- nimbleFunction( setup = function() {}, run = function(d = double(), x = double(1)) { ans <- exp(-d*x) - ans <- nf_sqrt(ans) - return(ans) + ans <- nf_sqrt(ans) + return(ans) returnType(double(1)) - }, + }, methods = list( derivsRun = function(d = double(), x = double(1)) { wrt <- 1:(1 + length(x)) # total length of d and x @@ -306,11 +306,11 @@ x <- c(2.1, 2.2) C_my_derivs_demo3$derivsRun(d, x)$jacobian ``` -Note that for a nimbleFunction without setup code, one can say `buildDerivs=TRUE`, `buildDerivs = 'run'`, or `buildDerivs = list(run = list())`. One can't take derivatives of `nf_sqrt` on its own, but it can be called by a method of a nimbleFunction that can have its derivatives taken (e.g. the `run` method of a `derivs_demo3` object). +Note that for a nimbleFunction without setup code, one can say `buildDerivs=TRUE`, `buildDerivs = 'run'`, or `buildDerivs = list(run = list())`. One can't take derivatives of `nf_sqrt` on its own, but it can be called by a method of a nimbleFunction that can have its derivatives taken (e.g. the `run` method of a `derivs_demo3` object). ### Understanding more about how AD works: *taping* of operations {#sec:understanding-more-AD} -At this point, it will be helpful to grasp more of how AD works and its implementation in nimble via the CppAD library [@bell-22]. AD methods work by following a set of calculations multiple times, sometimes in reverse order. +At this point, it will be helpful to grasp more of how AD works and its implementation in nimble via the CppAD library [@bell-22]. AD methods work by following a set of calculations multiple times, sometimes in reverse order. For example, the derivative of $y = \exp(x^2)$ can be calculated (by the chain rule of calculus) as the derivative of $y = \exp(z)$ (evaluated at $z = x^2$) times the derivative of $z = x^2$ (evaluated at $x$). This can be calculated by first going through the sequence of steps for the value, i.e. (i) $z = x^2$, (ii) $y = \exp(z)$, and then again for the derivatives, i.e. (i) $dz = 2x (dx)$, (ii) $dy = \exp(z) dz$. These steps determine the instantaneous change $dy$ that results from an instantaneous change $dx$. (It may be helpful to think of both values as a function of a third variable such as $t$. Then we are determining $dy/dt$ as a function of $dx/dt$.) In CppAD, the basic numeric object type (double) is replaced with a special type and corresponding versions of all the basic math operations so that those operations can track not just the values but also derivative information. @@ -395,7 +395,7 @@ If you are interested in measuring the performance of AD in nimble, please remem ## Advanced uses: double taping -Suppose you are interested only in the Jacobian, not in the value, and/or want only some elements of the Jacobian. You might still need the value, but obtaining the value alone from an AD tape (i.e. from `derivs`) will be slower than obtaining it by simply calling the function. On the other hand, AD methods need to calculate the value before calculating first order derivatives, so the value will be calculated anyway. However, in some cases, some of the steps of value calculations aren't really needed if one only wants, say, first-order derivatives. In addition, if not all elements of the Jacobian are wanted, then some unnecessary calculations will be done internally that one might want to avoid. +Suppose you are interested only in the Jacobian, not in the value, and/or want only some elements of the Jacobian. You might still need the value, but obtaining the value alone from an AD tape (i.e. from `derivs`) will be slower than obtaining it by simply calling the function. On the other hand, AD methods need to calculate the value before calculating first order derivatives, so the value will be calculated anyway. However, in some cases, some of the steps of value calculations aren't really needed if one only wants, say, first-order derivatives. In addition, if not all elements of the Jacobian are wanted, then some unnecessary calculations will be done internally that one might want to avoid. A way to cut out unnecessary calculations is to record a tape of a tape, which we call double taping. Let's see an example before explaining further. @@ -411,7 +411,7 @@ derivs_demo5 <- nimbleFunction( wrt = integer(1), reset = logical(0, default=FALSE)) { ans <- derivs(run(d, x), wrt = wrt, order = 1, reset = reset) - jac <- ans$jacobian[,1] # derivatives wrt 'd' only + jac <- ans$jacobian[,1] # derivatives wrt 'd' only return(jac) returnType(double(1)) }, @@ -458,7 +458,7 @@ Note that in any of the methods below, when order 0 (value) is requested, then a Recall that `model$calculate(nodes)` returns the sum of the log probabilities of all stochastic nodes in `nodes`. Deterministic calculations are also executed; they contribute 0 to the sum of probabilities but may be needed for inputs to subsequent calculations. Calculations are done in the order of `nodes`, which should be a valid order for the model, often obtained from `model$getDependencies`. -The simplest way to get derivatives for model calculations is to use `model$calculate` as the function taped by `derivs` (here shown by its alternative name `nimDerivs` for illustration). +The simplest way to get derivatives for model calculations is to use `model$calculate` as the function taped by `derivs` (here shown by its alternative name `nimDerivs` for illustration). ```{r eval=eval} derivs_nf <- nimbleFunction( @@ -488,11 +488,11 @@ Now we can make an instance of `derivs_nf`, compile the model and nimbleFunction ```{r, eval=eval} model_code <- nimbleCode({ - # priors + # priors intercept ~ dnorm(0, sd = 100) beta ~ dnorm(0, sd = 100) sigma ~ dhalfflat() - # random effects and data + # random effects and data for(i in 1:10) { # random effects ran_eff[i] ~ dnorm(0, sd = sigma) @@ -524,7 +524,7 @@ derivs_result <- cDerivs_all$run(order = 0:2) derivs_result ``` -As above, using `order = 0:2` results in the the value (0th order), Jacobian (1st order), and Hessian (2nd order). The function `model$calculate` is organized here to have inputs that are the current values of `intercept`, `beta`, and `sigma` in the model. It has output that is the summed log probability of the `calc_nodes`. The Jacobian columns are first derivatives with respect to `intercept`, `beta` and `sigma`, respectively. The first and second indices of the Hessian array follow the same ordering. For example `derivs_result$hessian[2,1,1]` is the second derivative with respect to `beta` (first index = 2) and `intercept` (second index = 1). +As above, using `order = 0:2` results in the the value (0th order), Jacobian (1st order), and Hessian (2nd order). The function `model$calculate` is organized here to have inputs that are the current values of `intercept`, `beta`, and `sigma` in the model. It has output that is the summed log probability of the `calc_nodes`. The Jacobian columns are first derivatives with respect to `intercept`, `beta` and `sigma`, respectively. The first and second indices of the Hessian array follow the same ordering. For example `derivs_result$hessian[2,1,1]` is the second derivative with respect to `beta` (first index = 2) and `intercept` (second index = 1). In the case of `model$calculate`, the first index of the Jacobian and the last index of the Hessian are always 1 because derivatives are of the first (and only) output value. @@ -542,7 +542,7 @@ We can see that all the results clearly match the compiled results to a reasonab ### Method 2: `nimDerivs` of a method that calls `model$calculate` -Sometimes one needs derivatives of calculations done in a nimbleFunction as well as in a model. And sometimes one needs to change values in a model before and/or after doing model calculations and have those changes recorded in the derivative tape. +Sometimes one needs derivatives of calculations done in a nimbleFunction as well as in a model. And sometimes one needs to change values in a model before and/or after doing model calculations and have those changes recorded in the derivative tape. For these reasons, it is possible to take derivatives of a method that includes a call to `model$calculate`. Continuing the above example, say we want to take derivatives with respect to the log of `sigma`, so the input vector will be treated as `intercept`, `beta`, and `log(sigma)`, in that order. We will convert `log(sigma)` to `sigma` before using it in the model, and we want the derivative of that transformation included in the tape (i.e. using the chain rule). @@ -556,7 +556,7 @@ derivs_nf2 <- nimbleFunction( derivsInfo <- makeModelDerivsInfo(model, wrt_nodes, calc_nodes) updateNodes <- derivsInfo$updateNodes constantNodes <- derivsInfo$constantNodes - n_wrt <- length(wrt_nodes) + n_wrt <- length(wrt_nodes) # If wrt_nodes might contain non-scalar nodes, the more general way # to determine the length of all scalar elements is: # length(model$expandNodeNames(wrt_nodes, returnScalarComponents = TRUE)) @@ -575,7 +575,7 @@ derivs_nf2 <- nimbleFunction( reset = logical(0, default=FALSE)) { wrt <- 1:n_wrt ans <- nimDerivs(run(x), wrt = wrt, order = order, reset = reset, - model = model, updateNodes = updateNodes, + model = model, updateNodes = updateNodes, constantNodes = constantNodes) return(ans) returnType(ADNimbleList()) @@ -617,7 +617,7 @@ In most cases, you can obtain `constantNodes` and `updateNodes` from `makeModelD - not in the `wrt` argument to `nimDerivs` (only relevant for Method 1), and - not assigned into the model using `values(model, nodes) <<- some_values` prior to `model$calculate` (only relevant for Method 2), and - not a deterministic node in `calc_nodes`. - + - `updateNodes` includes node names satisfying those conditions whose values might change between uses of the tape regardless of the `reset` argument. - `constantNodes` includes node names satisfying those conditions whose values will be baked into the tape (will not change) until the next call with `reset=TRUE`. @@ -684,7 +684,7 @@ The worked example next will illustrate use of the parameter transformation syst In this example, we will obtain maximum likelihood estimates of a Poisson GLM. It will be similar to the GLMM above, but without the random effects. Maximization of a likelihood function can be done with methods that don't require gradients, but it is much faster to use methods that do use gradients. -To illustrate the parameter transformation system, we will give `p` a uniform prior, constrained to values between 0 and 1. The maximum likelihood estimation below will not use the prior probability, but nevertheless the prior will be interpreted for its constraints on valid values in this case. This need for parameter transformation is somewhat artifical in the interest of simplicity. +To illustrate the parameter transformation system, we will give `p` a uniform prior, constrained to values between 0 and 1. The maximum likelihood estimation below will not use the prior probability, but nevertheless the prior will be interpreted for its constraints on valid values in this case. This need for parameter transformation is somewhat artificial in the interest of simplicity. ```{r eval=eval} model_code <- nimbleCode({ @@ -692,7 +692,7 @@ model_code <- nimbleCode({ p ~ dunif(0, 1) log_p <- log(p) beta ~ dnorm(0, sd = 100) - # random effects and data + # random effects and data for(i in 1:50) { # data y[i] ~ dpois(exp(log_p + beta*X[i])) @@ -770,7 +770,7 @@ logLikelihood_nf <- nimbleFunction( ) ``` -Sometimes you might want a call to `transformer$transform` or `transformer$inverseTransform` to be included in a derivative calculation. In this case, we don't. +Sometimes you might want a call to `transformer$transform` or `transformer$inverseTransform` to be included in a derivative calculation. In this case, we don't. Next we build and compile an instance of `logLikelihood_nf` specialized (by the `setup`` function) to our model and its parameters. @@ -779,12 +779,12 @@ ll_glm <- logLikelihood_nf(model, c("p", "beta")) Cll_glm <- compileNimble(ll_glm, project = model) ``` -Now let's see how the parameter transformation works for this model. For `beta`, it already has an unconstrained range of valid values, so it has no transformation (i.e. it has an indentity transformation). For `p`, it is constrained beween 0 and 1, so it will be logit-transformed. +Now let's see how the parameter transformation works for this model. For `beta`, it already has an unconstrained range of valid values, so it has no transformation (i.e. it has an identity transformation). For `p`, it is constrained between 0 and 1, so it will be logit-transformed. ```{r eval=eval} Cll_glm$transform(c(0.2, 1)) # The inverse transformation goes back to the original parameters, within rounding: -Cll_glm$inverseTransform(c(-1.386294, 1)) +Cll_glm$inverseTransform(c(-1.386294, 1)) ``` Now we are ready to use the methods provided in `Cll_glm` in a call to R's `optim` to find the MLE. diff --git a/UserManual/src/chapter_BNP.Rmd b/UserManual/src/chapter_BNP.Rmd index ecf5941fe..f346ccbd8 100644 --- a/UserManual/src/chapter_BNP.Rmd +++ b/UserManual/src/chapter_BNP.Rmd @@ -1,5 +1,5 @@ - @@ -7,12 +7,12 @@ ```{r, echo=FALSE} require(nimble) -``` +``` # Bayesian nonparametric models {#cha-bnp} ## Bayesian nonparametric mixture models {#sec:bnpmixtures} - + NIMBLE provides support for Bayesian nonparametric (BNP) mixture modeling. The current implementation provides support for hierarchical specifications involving Dirichlet process (DP) mixtures [@ferguson_73;@ferguson_74;@lo_84;@escobar_94;@escobar_west_95]. These allow one to avoid specifying a particular parametric distribution for a given node (parameter) in a model. Instead one can use a DP mixture as a much more general, nonparametric distribution. For example, a normal distribution for a random effect could be replaced by a DP mixture of normal distributions, with the number of components of the mixture being determined from the data and not fixed in advance. We'll first introduce the general, technical definition of a DP mixture model before describing the Chinese Restaurant Process representation, which may be more interpretable for many readers. More specifically, a DP mixture model for a random variable $y_i$ takes the form @@ -22,7 +22,7 @@ $$G \mid \alpha, G_0 \sim DP(\alpha, G_0),$$ where $h(\cdot \mid \theta)$ is a suitable kernel (i.e., probability density/mass function) with parameter $\theta$, and $\alpha$ and $G_0$ are the concentration and baseline distribution parameters of the DP, respectively. DP mixture models can be written with different levels of hierarchy, all being equivalent to the model above. While "y" would often be used as notation for a data value, it is used generically here, noting that often DP mixtures are used for random effects rather than directly for observations. -When the random distribution (also referred to as a random 'measure') $G$ is integrated out from the model, the DP mixture model can be written using latent or membership variables, $z_i$, following a Chinese Restaurant Process (CRP) distribution [@blackwell_mcqueen_73], discussed in Section \@ref(sec:crp). The model takes the form +When the random distribution (also referred to as a random 'measure') $G$ is integrated out from the model, the DP mixture model can be written using latent or membership variables, $z_i$, following a Chinese Restaurant Process (CRP) distribution [@blackwell_mcqueen_73], discussed in Section \@ref(sec:crp). The model takes the form $$y_i \mid \tilde{\boldsymbol{\theta}}, z_i \overset{ind}{\sim} h(\cdot \mid \tilde{\theta}_{z_i}),$$ $$\boldsymbol{z}\mid \alpha \sim \mbox{CRP}(\alpha),\hspace{0.5cm} \tilde{\theta}_j \overset{iid}{\sim}G_0,$$ @@ -33,7 +33,7 @@ If a stick-breaking representation [@sethuraman_94], discussed in section \@ref( $$y_i \mid {\boldsymbol{\theta}}^{\star}, \boldsymbol{v} \overset{ind}{\sim} \sum_{l=1}^{\infty}\left\{ v_l\prod_{m ```{r, mvConf} -mvConf = modelValuesConf(vars = c('a', 'b', 'c'), - type = c('double', 'int', 'double'), +mvConf = modelValuesConf(vars = c('a', 'b', 'c'), + type = c('double', 'int', 'double'), size = list(a = 2, b =c(2,2), c = 1) ) customMV = modelValues(mvConf, m = 2) @@ -99,7 +99,7 @@ can be used as arguments to setup code in nimbleFunctions. In the example above a modelValues object is passed to setup code, but a modelValues configuration can also be passed, with creation of -modelValues object(s) from the configuration done in setup code. +modelValues object(s) from the configuration done in setup code. ### Accessing contents of modelValues {#sec:access-cont-modelv} @@ -108,22 +108,22 @@ from R, and in fewer ways from NIMBLE. ```{r, mv-access} # sets the first row of a to (0, 1). R only. -customMV[['a']][[1]] <- c(0,1) +customMV[['a']][[1]] <- c(0,1) # sets the second row of a to (2, 3) -customMV['a', 2] <- c(2,3) +customMV['a', 2] <- c(2,3) # can access subsets of each row customMV['a', 2][2] <- 4 # accesses all values of 'a'. Output is a list. R only. -customMV[['a']] +customMV[['a']] # sets the first row of b to a matrix with values 1. R only. -customMV[['b']][[1]] <- matrix(1, nrow = 2, ncol = 2) +customMV[['b']][[1]] <- matrix(1, nrow = 2, ncol = 2) # sets the second row of b. R only. -customMV[['b']][[2]] <- matrix(2, nrow = 2, ncol = 2) +customMV[['b']][[2]] <- matrix(2, nrow = 2, ncol = 2) # make sure the size of inputs is correct # customMV['a', 1] <- 1:10 " @@ -132,13 +132,13 @@ customMV[['b']][[2]] <- matrix(2, nrow = 2, ncol = 2) ``` Currently, only the syntax `customMV["a", 2]` works in the NIMBLE -language, not `customMV[["a"][[2]]`. +language, not `customMV[["a"][[2]]`. We can query and change the number of rows using `getsize` and `resize`, respectively. These work in both R and NIMBLE. Note that we don't specify the variables in this case: all variables in a modelValues object will have the same number of rows. - + ```{r, resize-mv} getsize(customMV) resize(customMV, 3) @@ -156,7 +156,7 @@ names from every scalar element of variables (e.g. "b[1, 1]" ,"b[2, 1]", etc.). The rows of the modelValues will be the rows of the matrix, with any matrices or arrays converted to a vector based on column-major ordering. - + ```{r, as.matrix-mv} as.matrix(customMV, 'a') # convert 'a' as.matrix(customMV) # convert all variables @@ -224,9 +224,9 @@ each simulation in a row of its modelValues. `calcNodesMV` and nodes into the model, and then do their job of calculating or collecting log probabilities (densities), respectively. Each of these returns a numeric vector with the summed log probabilities of the -chosen nodes from each row. `calcNodesMV` will +chosen nodes from each row. `calcNodesMV` will save the log probabilities back into the modelValues object if - `saveLP = TRUE`, a run-time argument. + `saveLP = TRUE`, a run-time argument. Here are some examples: @@ -241,11 +241,11 @@ cCalcManyXDeps <- compileNimble(rCalcManyXDeps, project = simpleModel) cGetLogProbMany <- compileNimble(rGetLogProbMany, project = simpleModel) cSimManyXY$run(m = 5) # simulating 5 times -cCalcManyXDeps$run(saveLP = TRUE) # calculating +cCalcManyXDeps$run(saveLP = TRUE) # calculating cGetLogProbMany$run() # result <- as.matrix(cSimManyXY$mv) # extract simulated values ``` - + ## The nimbleList data structure {#sec:nimbleLists} nimbleLists provide a container for storing different types of objects in NIMBLE, similar to the list data structure in R. Before a nimbleList can be created and used, a *definition*^[The *configuration* for a modelValues object is the same concept as a *definition* here; in a future release of NIMBLE we may make the usage more consistent between modelValues and nimbleLists.] for that nimbleList must be created that provides the names, types, and dimensions of the elements in the nimbleList. nimbleList definitions must be created in R (either in R's global environment or in setup code), but the nimbleList @@ -290,7 +290,7 @@ mynf <- nimbleFunction( returnType(exampleNimListDef()) return(vals) }) - + # pass exampleNimList as argument to mynf mynf(exampleNimList) @@ -321,7 +321,7 @@ Note that definitions for inner, or nested, nimbleLists must be created before t ### Pre-defined nimbleList types {#sec:predef-nimbleLists} -Several types of nimbleLists are predefined in NIMBLE for use with particular kinds of nimbleFunctions. +Several types of nimbleLists are predefined in NIMBLE for use with particular kinds of nimbleFunctions. - `eigenNimbleList` and `svdNimbleList` are the return types of the `eigen` and `svd` functions (Section \@ref(sec:eigen-nimFunctions}). - `waicNimbleList` is the type provided for WAIC output from an MCMC. @@ -331,11 +331,11 @@ Several types of nimbleLists are predefined in NIMBLE for use with particular ki ### Using *eigen* and *svd* in nimbleFunctions {#sec:eigen-nimFunctions} -NIMBLE has two linear algebra functions that return nimbleLists. The `eigen` function takes a symmetic matrix, `x`, as an argument and returns a nimbleList of type `eigenNimbleList`. nimbleLists of type `eigenNimbleList` have two elements: `values`, a vector with the eigenvalues of `x`, and `vectors`, a square matrix with the same dimension as `x` whose columns are the eigenvectors of `x`. The `eigen` function has two additional arguments: `symmetric` and `only.values`. The `symmetric` argument can be used to specify if `x` is a symmetric matrix or not. If `symmetric = FALSE` (the default value), `x` will be checked for symmetry. Eigendecompositions in NIMBLE for symmetric matrices are both faster and more accurate. Additionally, eigendecompostions of non-symmetric matrices can have complex entries, which are not supported by NIMBLE. If a complex entry is detected, NIMBLE will issue a warning and that entry will be set to `NaN`. The `only.values` arument defaults to `FALSE`. If `only.values = TRUE`, the `eigen` function will not calculate the eigenvectors of `x`, leaving the `vectors` nimbleList element empty. This can reduce calculation time if only the eigenvalues of `x` are needed. +NIMBLE has two linear algebra functions that return nimbleLists. The `eigen` function takes a symmetric matrix, `x`, as an argument and returns a nimbleList of type `eigenNimbleList`. nimbleLists of type `eigenNimbleList` have two elements: `values`, a vector with the eigenvalues of `x`, and `vectors`, a square matrix with the same dimension as `x` whose columns are the eigenvectors of `x`. The `eigen` function has two additional arguments: `symmetric` and `only.values`. The `symmetric` argument can be used to specify if `x` is a symmetric matrix or not. If `symmetric = FALSE` (the default value), `x` will be checked for symmetry. Eigendecompositions in NIMBLE for symmetric matrices are both faster and more accurate. Additionally, eigendecompostions of non-symmetric matrices can have complex entries, which are not supported by NIMBLE. If a complex entry is detected, NIMBLE will issue a warning and that entry will be set to `NaN`. The `only.values` argument defaults to `FALSE`. If `only.values = TRUE`, the `eigen` function will not calculate the eigenvectors of `x`, leaving the `vectors` nimbleList element empty. This can reduce calculation time if only the eigenvalues of `x` are needed. + +The `svd` function takes an $n \times p$ matrix `x` as an argument, and returns a nimbleList of type `svdNimbleList`. nimbleLists of type `svdNimbleList` have three elements: `d`, a vector with the singular values of `x`, `u` a matrix with the left singular vectors of `x`, and `v`, a matrix with the right singular vectors of `x`. The `svd` function has an optional argument `vectors` which defaults to a value of `"full"`. The `vectors` argument can be used to specify the number of singular vectors that are returned. If `vectors = "full"`, `v` will be an $n \times n$ matrix and `u` will be an $p \times p$ matrix. If `vectors = "thin"`, `v` will be an$n \times m$ matrix, where $m = \min(n,p)$, and `u` will be an $m \times p$ matrix. If `vectors = "none"`, the `u` and `v` elements of the returned nimbleList will not be populated. -The `svd` function takes an $n \times p$ matrix `x` as an argument, and returns a nimbleList of type `svdNimbleList`. nimbleLists of type `svdNimbleList` have three elements: `d`, a vector with the singular values of `x`, `u` a matrix with the left singular vectors of `x`, and `v`, a matrix with the right singular vectors of `x`. The `svd` function has an optional argument `vectors` which defaults to a value of `"full"`. The `vectors` argument can be used to specify the number of singular vectors that are returned. If `vectors = "full"`, `v` will be an $n \times n$ matrix and `u` will be an $p \times p$ matrix. If `vectors = "thin"`, `v` will be an$n \times m$ matrix, where $m = \min(n,p)$, and `u` will be an $m \times p$ matrix. If `vectors = "none"`, the `u` and `v` elements of the returned nimbleList will not be populated. - nimbleLists created by either `eigen` or `svd` can be returned from a nimbleFunction, using `returnType(eigenNimbleList())` or `returnType(svdNimbleList())` respectively. nimbleLists created by `eigen` and `svd` can also be used within other nimbleLists by specifying the nimbleList element types as `eigenNimbleList()` and `svdNimbleList()`. The below example demonstrates the use of `eigen` and `svd` within a nimbleFunction. @@ -343,13 +343,13 @@ nimbleLists created by either `eigen` or `svd` can be returned from a nimbleFunc eigenListFunctionGenerator <- nimbleFunction( setup = function(){ demoMatrix <- diag(4) + 2 - eigenAndSvdListDef <- nimbleList(demoEigenList = eigenNimbleList(), + eigenAndSvdListDef <- nimbleList(demoEigenList = eigenNimbleList(), demoSvdList = svdNimbleList()) eigenAndSvdList <- eigenAndSvdListDef$new() }, run = function(){ # we will take the eigendecomposition and svd of a symmetric matrix - eigenAndSvdList$demoEigenList <<- eigen(demoMatrix, symmetric = TRUE, + eigenAndSvdList$demoEigenList <<- eigen(demoMatrix, symmetric = TRUE, only.values = TRUE) eigenAndSvdList$demoSvdList <<- svd(demoMatrix, vectors = 'none') returnType(eigenAndSvdListDef()) @@ -363,5 +363,3 @@ outputList$demoSvdList$d ``` The eigenvalues and singular values returned from the above function are the same since the matrix being decomposed was symmetric. However, note that both eigendecompositions and singular value decompositions are numerical procedures, and computed solutions may have slight differences even for a symmetric input matrix. - - diff --git a/UserManual/src/chapter_Laplace.Rmd b/UserManual/src/chapter_Laplace.Rmd index 7ea2ffbf6..9d80174c3 100644 --- a/UserManual/src/chapter_Laplace.Rmd +++ b/UserManual/src/chapter_Laplace.Rmd @@ -4,7 +4,7 @@ ```{r, echo=FALSE} require(nimble) require(nimbleQuad, quietly = TRUE, warn.conflicts = FALSE) -``` +``` # Laplace, AGHQ, and nested approximations {#cha-laplace} @@ -39,11 +39,11 @@ We'll re-introduce the simple Poisson Generalized Linear Mixed Model (GLMM) exam ```{r eval=TRUE} model_code <- nimbleCode({ - # priors + # priors intercept ~ dnorm(0, sd = 100) beta ~ dnorm(0, sd = 100) sigma ~ dhalfflat() - # random effects and data + # random effects and data for(i in 1:10) { # random effects ran_eff[i] ~ dnorm(0, sd = sigma) @@ -129,7 +129,7 @@ If one wants finer grain control over using the approximation, one can use the m ```{r eval=eval} # Get the Laplace approximation for one set of parameter values. -Cglmm_laplace$calcLaplace(c(0, 0, 1)) +Cglmm_laplace$calcLaplace(c(0, 0, 1)) Cglmm_laplace$gr_Laplace(c(0, 0, 1)) # Get the corresponding gradient. MLE <- Cglmm_laplace$findMLE(c(0, 0, 1)) # Find the (approximate) MLE. MLE$par # MLE parameter values @@ -156,9 +156,9 @@ To run a regularized regression that uses the prior as a penalty but excludes th When finding the MLE via Laplace approximation or adaptive Gauss-Hermite quadrature (AGHQ), there are two numerical optimizations: (1) maximizing the joint log-likelihood of random effects and data given a set of parameter values to construct the approximation to the marginal log-likelihood at the given parameter values, and (2) maximizing the approximation to the marginal log-likelihood over the parameter values. Optimization (1) is the "inner" optimization and optimization (2) is the "outer" optimization. -Finding the MLE via Laplace approximation may be sensitive to the optimization methods used, in particular the choice of optimizer for the inner optimization, and the "BFGS" optimizer available through `optim()` may not perform well for inner optimization. +Finding the MLE via Laplace approximation may be sensitive to the optimization methods used, in particular the choice of optimizer for the inner optimization, and the "BFGS" optimizer available through `optim()` may not perform well for inner optimization. -As of version 1.3.0, the default choices for both the inner and outer optimization use R's `nlminb` optimizer. +As of version 1.3.0, the default choices for both the inner and outer optimization use R's `nlminb` optimizer. Users can choose a different optimizer for both of the optimizations. To change the inner or outer optimizers, one can use the `innerOptimMethod` and `outerOptimMethod` elements of the `control` list argument to `buildLaplace`. One can modify various settings that control the behavior of the inner and outer optimizers via `control` as well. See `help(buildLaplace)` for more details. @@ -177,7 +177,7 @@ library(Matrix) ## and wraps the optimizer of interest. nimbleTMBnewton <- function(par, fn, gr, he, lower, upper, control, hessian) { ## Wrap `he` as return value needs to be of class `dsCMatrix`. - he_matrix <- function(p) Matrix(he(p), doDiag = FALSE, sparse = TRUE) + he_matrix <- function(p) Matrix(he(p), doDiag = FALSE, sparse = TRUE) invalid <- function(x) is.null(x) || is.na(x) || is.infinite(x) if(invalid(control$trace)) control$trace <- 1 if(invalid(control$maxit)) control$maxit <- 100 @@ -223,18 +223,18 @@ Here we list some of the similarities and differences from INLA and the ELGM app - For inference on the latent nodes, we provide joint simulation from a multivariate normal mixture over the (hyper)parameter grid points as done in ELGM and also available in INLA. Unlike in INLA, we do not provide univariate latent inference using deterministic nested Laplace approximation. The simulation-based approach may not be as accurate, but it allows for joint inference, including inference on quantities that depend on more than one latent node. - Unlike either ELGM or INLA, latent nodes are not required to have a joint normal distribution, though accuracy may be less when the latent nodes have other distributions. - For latent nodes whose conditional distributions factor into univariate conditionally independent sets, the Laplace approximation is a product of univariate approximations, and one can instead use NIMBLE's AGHQ approximation for higher accuracy. - - We allow the user to choose the grid used for the (hyper)parameters. By default for $d>2$ parameters, we use the CCD grid used by INLA, but one can choose to use the AGHQ grid as used in ELGM or provide one's own grid. - + - We allow the user to choose the grid used for the (hyper)parameters. By default for $d>2$ parameters, we use the CCD grid used by INLA, but one can choose to use the AGHQ grid as used in ELGM or provide one's own grid. + ### Overview of the methodology -NIMBLE's nested approximation consists of several pieces, with different options that allow a user to choose how the approximation is done. We briefly describe the pieces here. For simplicity, we refer to the (hyper)parameters simply as "parameters". +NIMBLE's nested approximation consists of several pieces, with different options that allow a user to choose how the approximation is done. We briefly describe the pieces here. For simplicity, we refer to the (hyper)parameters simply as "parameters". #### Marginalization over the latent nodes -We approximately marginalize (integrate) over the latent nodes to approximate the marginal joint distribution of the parameters. This uses Laplace approximation and is sometimes referred to as the "inner" marginalization. The Laplace approximation is computed using the gradient and Hessian with respect to the latent nodes for a given set of parameter values. +We approximately marginalize (integrate) over the latent nodes to approximate the marginal joint distribution of the parameters. This uses Laplace approximation and is sometimes referred to as the "inner" marginalization. The Laplace approximation is computed using the gradient and Hessian with respect to the latent nodes for a given set of parameter values. -In the case that the conditional distributions of the latent nodes factor into univariate conditionally-independent sets of nodes (conditional on the parameters), this approximation by default uses the product of univariate Laplace approximations and in NIMBLE can be made more accurate by using AGHQ with more than one quadrature point. +In the case that the conditional distributions of the latent nodes factor into univariate conditionally-independent sets of nodes (conditional on the parameters), this approximation by default uses the product of univariate Laplace approximations and in NIMBLE can be made more accurate by using AGHQ with more than one quadrature point. #### Approximating the marginal parameter density on a grid @@ -244,7 +244,7 @@ NIMBLE primarily offers the options of using a CCD grid (as used by INLA) or an For $d <= 2$ NIMBLE defaults to the AGHQ grid and otherwise uses the CCD grid. -NIMBLE also allows users to provide their own grid. +NIMBLE also allows users to provide their own grid. #### Joint inference for latent nodes @@ -284,7 +284,7 @@ Note that the treatment of fixed effects differs between the nested approximatio It is also possible to choose which nodes go in which set (Section \@ref(subsec:nested-node-determ)). One might choose to include fixed effects in the parameter set, akin to how Laplace approximation works. Another situation where one might configure this manually is if random effects in a model are specified in a "non-centered" parameterization such as: $$ -\eta_i = \mu + \sigma b_i +\eta_i = \mu + \sigma b_i $$ $$ b_i \sim \mathcal{N}(0, 1) @@ -292,19 +292,19 @@ $$ In this parameterization, the random effects, $b_i$, do not depend on any parameters because $\mu$ and $\sigma$ are involved directly in the linear predictor, $\eta_i$. This stands in contrast to the common centered parameterization, with $b_i \sim \mathcal{N}(\mu, \sigma)$. Because each $b_i$ does not depend on any parameters, NIMBLE's nested approximation node determination algorithm would by default put the $b_i$ nodes into the parameter set, which in general would result in very slow computation (because there would be a large number of parameters). (And in reality this would generally lead to a failure to build the algorithm because there would be no latent nodes found.) -One further note is that it can be advantageous computationally to have fixed effects in the parameter set. This can sometimes allow NIMBLE to set up a product of lower-dimensional Laplace approximations instead of one higher-dimensional Laplace approximation. However, the benefit of having a product of Laplace approximations trades off with the higher-dimensionality of the parameter vector, which increases computation. How these trade off in a particular modeling context can be worth experimentation by the user. +One further note is that it can be advantageous computationally to have fixed effects in the parameter set. This can sometimes allow NIMBLE to set up a product of lower-dimensional Laplace approximations instead of one higher-dimensional Laplace approximation. However, the benefit of having a product of Laplace approximations trades off with the higher-dimensionality of the parameter vector, which increases computation. How these trade off in a particular modeling context can be worth experimentation by the user. #### Computational bottlenecks -To sample the latent nodes, the algorithm needs to find the weights and the parameters of the multivariate normal approximation to the latent nodes at each point in the grid of parameter values. For a large number of parameters this can be expensive, because Laplace approximation is done at many grid points. For a large number of latent node elements, even a single Laplace approximation at a single set of parameter values can be expensive. Computing the parameters of this approximation involves optimization over a space whose dimension is the number of latent node elements, as well as computation of the Hessian at the maximum. Furthermore, for inference, one then needs to simulate from the high-dimensional multivariate normal. +To sample the latent nodes, the algorithm needs to find the weights and the parameters of the multivariate normal approximation to the latent nodes at each point in the grid of parameter values. For a large number of parameters this can be expensive, because Laplace approximation is done at many grid points. For a large number of latent node elements, even a single Laplace approximation at a single set of parameter values can be expensive. Computing the parameters of this approximation involves optimization over a space whose dimension is the number of latent node elements, as well as computation of the Hessian at the maximum. Furthermore, for inference, one then needs to simulate from the high-dimensional multivariate normal. -To approximate the univariate marginals for the parameters via AGHQ, this requires $d-1$-dimensional AGHQ, which can be expensive, because the number of quadrature points grows as $(d-1)^k$ where $k$ is the number of Gauss-Hermite quadrature grid points in one dimension. Often this would be chosen to be small, such as 3 or 5, but even then the computations increase rapidly in $d$. +To approximate the univariate marginals for the parameters via AGHQ, this requires $d-1$-dimensional AGHQ, which can be expensive, because the number of quadrature points grows as $(d-1)^k$ where $k$ is the number of Gauss-Hermite quadrature grid points in one dimension. Often this would be chosen to be small, such as 3 or 5, but even then the computations increase rapidly in $d$. ### Using NIMBLE's nested approximation Next we'll give example usage of NIMBLE's nested approximation. Further details are available in the usual R help content for the various functions. -The core functions are `buildNestedApprox` and `runNestedApprox`. These are analagous to the "build" and "run" functions for MCMC and for Laplace. +The core functions are `buildNestedApprox` and `runNestedApprox`. These are analogous to the "build" and "run" functions for MCMC and for Laplace. - `buildNestedApprox` sets up the approximation for a model of interest and allows the user to control various aspects of how the approximation is done. It returns an uncompiled nested approximation algorithm. - `runNestedApprox` runs the basic steps of the approximation, giving initial univariate marginal parameter inference (and optionally allowing the user to request latent node or parameter samples). It returns a nested approximation (`nestedApprox`) object. @@ -319,7 +319,7 @@ The core functions are `buildNestedApprox` and `runNestedApprox`. These are anal #### Example use -We'll use the penicillin example from the `faraway` package, which has data on penicillin production as a function of treatment (four levels) and blend (five levels). The treatment is considered as a fixed effect (note the constant, large variance/small precision for `b[i]`) while the blend/block is considered as a random effect. Note that because of the normal likelihood and normal priors for the latent nodes, the conditional distribution for the latent nodes given the data and hyperparameters is a multivariate normal, so the Laplace approximation in this case is exact (and one could also marginalize analytically before carrying out maximum likelihood estimation or MCMC). In many uses of nested approximation, the likelihood is not normal, but the latent node distribution is. +We'll use the penicillin example from the `faraway` package, which has data on penicillin production as a function of treatment (four levels) and blend (five levels). The treatment is considered as a fixed effect (note the constant, large variance/small precision for `b[i]`) while the blend/block is considered as a random effect. Note that because of the normal likelihood and normal priors for the latent nodes, the conditional distribution for the latent nodes given the data and hyperparameters is a multivariate normal, so the Laplace approximation in this case is exact (and one could also marginalize analytically before carrying out maximum likelihood estimation or MCMC). In many uses of nested approximation, the likelihood is not normal, but the latent node distribution is. ```{r} data(penicillin, package="faraway") @@ -362,7 +362,7 @@ Next we run our nested approximation, getting back initial inference on the para ```{r} # To prefer fixed rather than scientific notation for easier viewing. -options(scipen = 2) +options(scipen = 2) result <- runNestedApprox(comp_approx) result ``` @@ -550,4 +550,3 @@ approx <- buildNestedApprox(model, paramNodes = c("beta0", "beta1", "sigma")) ``` Note the messaging indicating that 30 Laplace approximations (one for each conditionally-independent set of nodes containing a ${y_i, b_i}$ pair) were built. (In other, perhaps more common situations, one might have grouping structure such that a single random effect is associated with multiple observations such that the random effect and associated observations are conditionally independent of other random effects and observations given the parameters.) - diff --git a/UserManual/src/chapter_MCMC.Rmd b/UserManual/src/chapter_MCMC.Rmd index 60234d1e7..592a6c9c2 100644 --- a/UserManual/src/chapter_MCMC.Rmd +++ b/UserManual/src/chapter_MCMC.Rmd @@ -2,7 +2,7 @@ ```{r, echo=FALSE} require(nimble) -``` +``` # MCMC {#cha-mcmc} @@ -13,7 +13,7 @@ if(!require(nimble, warn.conflicts = FALSE, quietly = TRUE)) { } require(methods, warn.conflicts = FALSE, quietly = TRUE) # seems to be needed, but why? require(igraph, warn.conflicts = FALSE, quietly = TRUE) # same question -``` +``` NIMBLE provides a variety of paths to creating and executing an MCMC algorithm, which differ greatly in their simplicity of use, and also in the options available and customizability. @@ -23,10 +23,10 @@ The lengthier and more customizable approach to invoking the MCMC engine on a pa 1. (Optional) Create and customize an MCMC configuration for a particular model: - + a. Use `configureMCMC` to create an MCMC configuration (see Section \@ref(sec:mcmc-configuration)). The configuration contains a list of samplers with the node(s) they will sample. - a. (Optional) Customize the MCMC configuration: - + a. (Optional) Customize the MCMC configuration: + i. Add, remove, or re-order the list of samplers (Section \@ref(sec:samplers-provided) and `help(samplers)` in R for details), including adding your own samplers (Section \@ref(sec:user-samplers)); i. Change the tuning parameters or adaptive properties of individual samplers; i. Change the variables to monitor (record for output) and thinning intervals for MCMC samples. @@ -50,7 +50,7 @@ End-to-end examples of MCMC in NIMBLE can be found in Sections \@ref(sec:creatin - @@ -61,10 +61,10 @@ End-to-end examples of MCMC in NIMBLE can be found in Sections \@ref(sec:creatin - - + @@ -127,7 +127,7 @@ mcmc_out <- nimbleMCMC(code = code, data = data, inits = initsFunction, initsList <- initsFunction() Rmodel <- nimbleModel(code, data = data, inits = initsList) -# Using the existing Rmodel object, execute three MCMC chains with +# Using the existing Rmodel object, execute three MCMC chains with # specified burn-in. Return samples, overall log-probability, # summary statistics, and WAIC. mcmc_out <- nimbleMCMC(model = Rmodel, @@ -192,7 +192,7 @@ The behavior of `configureMCMC` can be customized to control how samplers are as #### Sampling posterior predictive nodes {#sec:post-pred-sampling} -A posterior predictive node is a node that is not itself data and has no data nodes in its entire downstream (descendant) dependency network. Such nodes play no role in inference for model parameters but have often been included in BUGS models to accomplish posterior predictive checks and similar calculations. +A posterior predictive node is a node that is not itself data and has no data nodes in its entire downstream (descendant) dependency network. Such nodes play no role in inference for model parameters but have often been included in BUGS models to accomplish posterior predictive checks and similar calculations. As of version 0.13.0, NIMBLE's handling of posterior predictive nodes in MCMC sampling has changed in order to improve MCMC mixing. Samplers for nodes that are not posterior predictive nodes no longer condition on the values of the posterior predictive nodes. This produces a valid MCMC over the posterior distribution marginalizing over the posterior predictive nodes. This MCMC will generally mix better than an MCMC that conditions on the values of posterior predictive nodes, by reducing the dimensionality of the parameter space and removing the dependence between the sampled nodes and the posterior predictive nodes. At the end of each MCMC iteration, the posterior predictive nodes are sampled by `posterior_predictive` sampler(s) based on their conditional distribution(s). @@ -203,14 +203,14 @@ Note that predictive node sampling is done at each iteration. When thinning is u Very basic control of default sampler assignments is provided via two arguments to `configureMCMC`. The `useConjugacy` argument controls whether conjugate samplers are assigned when possible, and the `multivariateNodesAsScalars` argument controls whether scalar elements of multivariate nodes are sampled individually. See `help(configureMCMC)` for usage details. - - + #### Default monitors -The default MCMC configuration includes monitors on all top-level stochastic nodes of the model. Only variables that are monitored will have their samples saved for use outside of the MCMC. MCMC configurations include two sets of monitors, each with different thinning intervals. By default, the second set of monitors (`monitors2`) is empty. +The default MCMC configuration includes monitors on all top-level stochastic nodes of the model. Only variables that are monitored will have their samples saved for use outside of the MCMC. MCMC configurations include two sets of monitors, each with different thinning intervals. By default, the second set of monitors (`monitors2`) is empty. #### Automated parameter blocking @@ -222,7 +222,7 @@ The default configuration may be replaced by one generated from an automated par autoBlockConf <- configureMCMC(Rmodel, autoBlock = TRUE) ``` -Note that this using `autoBlock = TRUE` compiles and runs MCMCs, progressively exploring different sampler assignments, so it takes some time and generates some output. It is most useful for determining effective blocking strategies that can be re-used for later runs. The additional control argument `autoIt` may also be provided to indicate the number of MCMC samples to be used in each trial of the automated blocking procedure (default 20,000). +Note that this using `autoBlock = TRUE` compiles and runs MCMCs, progressively exploring different sampler assignments, so it takes some time and generates some output. It is most useful for determining effective blocking strategies that can be re-used for later runs. The additional control argument `autoIt` may also be provided to indicate the number of MCMC samples to be used in each trial of the automated blocking procedure (default 20,000). ### Customizing the MCMC configuration {#sec:customizing-mcmc-conf} @@ -234,8 +234,8 @@ One can create an MCMC configuration with default samplers on just a particular If a data node is included in `nodes`, *it will be assigned a sampler*. This is the only way in which a default sampler may be placed on a data node and will result in overwriting data values in the node. -#### Creating an empty configuration -If you plan to customize the choice of all samplers, it can be useful to obtain a configuration with no sampler assignments at all. This can be done by any of `nodes = NULL`, `nodes = character()`, or `nodes = list()`. +#### Creating an empty configuration +If you plan to customize the choice of all samplers, it can be useful to obtain a configuration with no sampler assignments at all. This can be done by any of `nodes = NULL`, `nodes = character()`, or `nodes = list()`. #### Overriding the default sampler control list values -The default values of control list elements for all sampling algorithms may be overridden through use of the `control` argument to `configureMCMC`, which should be a named list. +The default values of control list elements for all sampling algorithms may be overridden through use of the `control` argument to `configureMCMC`, which should be a named list. Named elements in the `control` argument will be used for all default samplers and any subsequent sampler added via `addSampler` (see below). For example, the following will create the default MCMC configuration, except all `RW` samplers will have their initial `scale` set to 3, and none of the samplers (`RW`, or otherwise) will be adaptive. ```{r, mcmcconf6, eval=FALSE} @@ -359,7 +359,7 @@ samplerConfList <- mcmcConf$getSamplers() samplerConfList[[1]]$setName("newNameForThisSampler") # change the sampler function of the second sampler, -# assuming existance of a nimbleFunction 'anotherSamplerNF', +# assuming existence of a nimbleFunction 'anotherSamplerNF', # which represents a valid nimbleFunction sampler. samplerConfList[[2]]$setSamplerFunction(anotherSamplerNF) @@ -411,7 +411,7 @@ The variables specified in `monitors` and `monitors2` will be recorded (with th Monitors may be added to the MCMC configuration either in the original call to `configureMCMC` or using the `addMonitors` method: ```{r, addMonitors, eval=FALSE} # Using an argument to configureMCMC -mcmcConf <- configureMCMC(Rmodel, monitors = c("alpha", "beta"), +mcmcConf <- configureMCMC(Rmodel, monitors = c("alpha", "beta"), monitors2 = "x") # Calling a member method of the mcmcconf object @@ -439,7 +439,7 @@ mcmcConf$setThin2(100) The current lists of monitors and thinning intervals may be displayed using the `printMonitors` method. Both sets of monitors (`monitors` and `monitors2`) may be reset to empty character vectors by calling the `resetMonitors` method. The methods `getMonitors` and `getMonitors2` return the currently specified `monitors` and `monitors2` as character vectors. -#### Monitoring model log-probabilities +#### Monitoring model log-probabilities One can record log-probabilities for all the individual nodes in a given variable as discussed next. Alternatively, one can monitor log-probabilities for individual nodes or summed log-probabilities for sets of nodes using NIMBLE's derived quantities system (Section \@ref(sec:derived-quantities). @@ -452,14 +452,14 @@ Cmodel <- compileNimble(Rmodel) Cmcmc <- compileNimble(Rmcmc, project = Rmodel) Cmcmc$run(10000) samples <- as.matrix(Cmcmc$mvSamples) -``` +``` The `samples` matrix will contain both MCMC samples and model log-probabilities. #### Using numeric samples to define a prior distribution -A set of numeric samples, perhaps generated from another MCMC algorithm, can be used to define the prior distribution of model nodes. This is accomplished using the `prior_samples` MCMC sampler. When assigning the `prior_samples` sampler to a `target` node, you must also provide a vector of numeric values (when `target` is a scalar node), or a matrix of values in the multidimenional case. A new value (or rows of values) is selected either sequentially (or randomly) from this numeric vector/matrix, and assigned into the `target` node on each MCMC iteration. See `help(prior_samples)` for more details: +A set of numeric samples, perhaps generated from another MCMC algorithm, can be used to define the prior distribution of model nodes. This is accomplished using the `prior_samples` MCMC sampler. When assigning the `prior_samples` sampler to a `target` node, you must also provide a vector of numeric values (when `target` is a scalar node), or a matrix of values in the multidimensional case. A new value (or rows of values) is selected either sequentially (or randomly) from this numeric vector/matrix, and assigned into the `target` node on each MCMC iteration. See `help(prior_samples)` for more details: ```{r, eval=FALSE} mcmcConf <- configureMCMC(Rmodel) mcmcConf$addSampler(target = 'theta', type = 'prior_samples', samples = rnorm(100)) @@ -469,17 +469,17 @@ mcmcConf$addSampler(target = 'theta', type = 'prior_samples', samples = rnorm(10 As of version 1.4.0, NIMBLE allows one to specify "derived quantities" to calculate or record additional quantities of interest. Potential derived quantities include: - - posterior means, - - posterior variances, + - posterior means, + - posterior variances, - log-probabilities (log-density values) for single or (summed) groups of nodes, - posterior predictive samples for predictive nodes, and - arbitrary calculations by providing a user-defined nimbleFunction that follows the required specification. User-defined derived quantity nimbleFunctions provide a general way for a user to intervene at the end of an MCMC iteration and do something. This would often be recording a value (to be returned at the end of the MCMC) or updating some summary statistic, but the derived quantities system doesn't require this, giving users flexibility. -The operation of every derived quantity function is governed by an `interval` parameter. Each derived quantity function is executed at the end of every `interval` MCMC iterations. The default value for `interval` is given by the `thin` interval of the MCMC. This means that, by default, each derived quantity function will execute after all the MCMC samplers for the MCMC iterations for which samples are recorded. If `interval = 1` for a derived quantity function, then that derived quantity function will execute at the end of every MCMC sampling iteration. +The operation of every derived quantity function is governed by an `interval` parameter. Each derived quantity function is executed at the end of every `interval` MCMC iterations. The default value for `interval` is given by the `thin` interval of the MCMC. This means that, by default, each derived quantity function will execute after all the MCMC samplers for the MCMC iterations for which samples are recorded. If `interval = 1` for a derived quantity function, then that derived quantity function will execute at the end of every MCMC sampling iteration. -Operation of derived quantity functions is entirely post-burnin. Derived quantity functions are not executed during the burnin period. +Operation of derived quantity functions is entirely post-burnin. Derived quantity functions are not executed during the burnin period. #### Specifying derived quantities @@ -507,10 +507,10 @@ The `logProb` derived quantity function calculates and stores the log-probabilit - any sets of nodes (summed), or - all stochastic model nodes (summed). -When the `nodes` argument is a character vector of node names, the log-probability of each node is recorded individually. We'll provide some examples using an old friend, the `pump` model. +When the `nodes` argument is a character vector of node names, the log-probability of each node is recorded individually. We'll provide some examples using an old friend, the `pump` model. ```{r} -pumpCode <- nimbleCode({ +pumpCode <- nimbleCode({ for(i in 1:N) { theta[i] ~ dgamma(alpha, beta) lambda[i] <- theta[i]*t[i] @@ -660,7 +660,7 @@ The predictive function accepts a `nodes` argument, which has the names of the p If nodes is omitted, the default set of nodes will be all (non-data) terminal nodes (i.e., nodes in the model graph that have no children). If `nodes` is specified as `nodes = '.all'`, then all predictive model nodes (including both stochastic and deterministic and terminal and non-terminal) will be simulated and saved. -The `predictive` derived quantity function will automatically determine the set of nodes that need to be simulated (based on the model graph structure) in order to correctly simulate the `nodes` specified. Alternatively, for additional user control in unusual situations, one can specify which nodes should be simulated using the `simNodes` argument. The set of nodes specified in `simNodes` will be simulated, in topological order, prior to saving the values of nodes. If simulation in the exact order specified in the `simNodes` argument is required, one may also provide the argument `sort = FALSE`. +The `predictive` derived quantity function will automatically determine the set of nodes that need to be simulated (based on the model graph structure) in order to correctly simulate the `nodes` specified. Alternatively, for additional user control in unusual situations, one can specify which nodes should be simulated using the `simNodes` argument. The set of nodes specified in `simNodes` will be simulated, in topological order, prior to saving the values of nodes. If simulation in the exact order specified in the `simNodes` argument is required, one may also provide the argument `sort = FALSE`. The automatic determination of the set of nodes to simulate assumes a valid MCMC assignment of samplers to unobserved model nodes. If one removes samplers from the default MCMC configuration or otherwise does non-standard modifications to the set of MCMC samplers, then the automatic determination of nodes that need to be simulated may be incorrect. In this case, both the `nodes` and `simNodes` arguments should be carefully specified to ensure correct calculation. @@ -670,7 +670,7 @@ We introduce a toy model to demonstrate the functionality of the predictive deri toyCode <- nimbleCode({ mu ~ dunif(0, 10) ## parameter y ~ dnorm(mu, 1) ## data - muSq <- mu^2 ## determistic, non-terminal predictive node + muSq <- mu^2 ## deterministic, non-terminal predictive node z ~ dnorm(muSq, 1) ## stochastic, terminal predictive node }) @@ -740,9 +740,9 @@ Some additional comments are provided within this skeleton for a derived quantit derived_name <- nimbleFunction( name = 'derived_name', ## All derived functions must contain (inherit from) the `derived_BASE` class. - contains = derived_BASE, + contains = derived_BASE, setup = function(model, mcmc, interval, control) { - ## + ## ## setup function must have these arguments, and has unrestricted access to: ## - model ## - mcmc @@ -755,13 +755,13 @@ derived_name <- nimbleFunction( ## defaultValue = character()) }, run = function(timesRan = double()) { - ## - ## run fuction has one argument: + ## + ## run function has one argument: ## timesRan, which is the number of times this method has executed during ## the present MCMC chain, ## including the present time. This can be used as an index for a matrix ## of results. - ## + ## }, methods = list( set_interval = function(newInterval = double()) { @@ -770,13 +770,13 @@ derived_name <- nimbleFunction( ## exactly. ## this method is necessary, to allow the default value of 'interval' to ## match the thinning interval of the MCMC. - ## + ## interval <<- newInterval ## }, before_chain = function(niter = double(), nburnin = double(), thin = double(1), chain = double()) { - ## + ## ## (optional) `before_chain` method, ## must have must have these arguments specifications: ## - niter (total number of *post-burnin* MCMC iterations) @@ -784,33 +784,33 @@ derived_name <- nimbleFunction( ## - thin (a length=2 vector, containing the thinning intervals for ## mvSamples and mvSamples2) ## - chain (the chain number of the current MCMC chain being run) - ## + ## }, after_chain = function() { - ## + ## ## (optional) `after_chain` method, ## no arguments and no return value - ## + ## }, get_results = function() { - ## + ## ## REQUIRED `get_results` method, ## which must return a dim=2 numeric array - ## + ## returnType(double(2)) }, get_names = function() { - ## + ## ## REQUIRED `get_names` method, ## which must return a dim=1 character vector - ## + ## returnType(character(1)) }, reset = function() { - ## + ## ## (optional) `reset` method ## no arguments and no return value - ## + ## } ) ) @@ -829,7 +829,7 @@ Once the MCMC configuration object has been created, and customized to one's lik Rmcmc <- buildMCMC(mcmcConf) ``` -`buildMCMC` is a nimbleFunction. The returned object `Rmcmc` is an instance of the nimbleFunction specific to configuration `mcmcConf` (and of course its associated model). +`buildMCMC` is a nimbleFunction. The returned object `Rmcmc` is an instance of the nimbleFunction specific to configuration `mcmcConf` (and of course its associated model). Note that if you would like to be able to calculate the WAIC of the model, you should usually set `enableWAIC = TRUE` as an argument to`configureMCMC` (or to `buildMCMC` if not using `configureMCMC`), or set `nimbleOptions(MCMCenableWAIC = TRUE)`, which will enable WAIC calculations for all subsequently built MCMC functions. For more information on WAIC calculations, including situations in which you can calculate WAIC without having set `enableWAIC = TRUE` see Section \@ref(sec:WAIC) or `help(waic)` in R. @@ -871,7 +871,7 @@ Not doing so can cause slow convergence or even failure of an MCMC algorithm (e. The following cases are particularly important to consider when initializing: - In a model with flat priors (i.e., using `x~dflat()` or `x~dhalfflat()`), NIMBLE cannot generate initial values from those priors. - - In a model with diffuse priors, initializing from the prior can give unreasonable/extreme initial values. + - In a model with diffuse priors, initializing from the prior can give unreasonable/extreme initial values. - In a model with stochastic indices (e.g., `x[idx[i]]` with `idx[i]` unknown), those indices should have (valid) initial values. - In a model with constraints (via `dconstraint`), the model should be initialized such that the constraints are satisfied. - In a model with censoring (via `dinterval`), the initial value(s) of `t` in `dinterval` for the censored node(s) should be consistent with the censoring that is specified. @@ -933,7 +933,7 @@ initsList <- list(list(mu = 1, sigma = 1), list(mu = 2, sigma = 10)) mcmc_out <- runMCMC(Cmcmc, nchains = 2, inits = initsList) -# run ten chains of 100,000 iterations each, using a function to +# run ten chains of 100,000 iterations each, using a function to # generate initial values and a fixed random number seed for each chain. # only return summary statistics from each chain; not all the samples. initsFunction <- function() @@ -965,7 +965,7 @@ This means that one can begin a new run of an existing MCMC without having to re In contrast, when `mcmc$run(niter, reset = FALSE)` is called, the MCMC picks up from where it left off, continuing the previous chain and expanding the output as needed. No values in the model are checked or altered, and sampler functions are not reset to their initial states. (Similarly, if using WAIC, one can use `resetWAIC = FALSE` so that the WAIC calculation is based on all the samples from the expanded set of samples.) -The `run` method also has an optional `resetMV` argument. This argument is only considered when`'reset` is set to `FALSE`. When `mcmc$run(niter, reset = FALSE, resetMV = TRUE)` is called the internal modelValues objects `mvSamples` and `mvSamples2` are each resized to the appropriate length for holding the requested number of samples (`niter/thin`, and `niter/thin2`, respectively) and the MCMC carries on from where it left off. In other words, the previously obtained samples are deleted (e.g. to reduce memory usage) prior to continuing the MCMC. The default value of `resetMV` is `FALSE`. +The `run` method also has an optional `resetMV` argument. This argument is only considered when`'reset` is set to `FALSE`. When `mcmc$run(niter, reset = FALSE, resetMV = TRUE)` is called the internal modelValues objects `mvSamples` and `mvSamples2` are each resized to the appropriate length for holding the requested number of samples (`niter/thin`, and `niter/thin2`, respectively) and the MCMC carries on from where it left off. In other words, the previously obtained samples are deleted (e.g. to reduce memory usage) prior to continuing the MCMC. The default value of `resetMV` is `FALSE`. @@ -977,13 +977,13 @@ If you want to obtain the computation time spent in each sampler, you can set `t ```{r, eval=FALSE} Cmcmc$run(niter, time = TRUE) Cmcmc$getTimes() -``` -will return a vector of the total time spent in each sampler, measured in seconds. - @@ -1061,7 +1061,7 @@ First, by default NIMBLE provides the conditional WAIC, namely the version of WA Second, WAIC relies on a partition of the observations, i.e., 'pointwise' prediction. By default, in NIMBLE 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. If one wants to group data nodes such that the joint density within each group is used, one can provide the `waicControl` list with a `dataGroups` element to `configureMCMC` or `buildMCMC` (when providing a model as the argument to `buildMCMC`). See `help(waic)` for more details. -Note that based on a limited set of simulation experiments in @Hug_Paciorek_2021, our tentative recommendation is that users only use marginal WAIC if also using grouping. +Note that based on a limited set of simulation experiments in @Hug_Paciorek_2021, our tentative recommendation is that users only use marginal WAIC if also using grouping. Marginal WAIC requires using Monte Carlo simulation at each iteration of the MCMC to average over draws for the latent variables. To assess the stability of the marginal WAIC to the number of Monte Carlo iterations, one can examine how the WAIC changes with increasing iterations (up to the full number of iterations specified via the `nItsMarginal` element of the `waicControl` list) based on the `WAIC_partialMC`, `lppd_partialMC`, and `pWAIC_partialMC` elements of the detailed WAIC output. @@ -1083,7 +1083,7 @@ This approach only provides conditional WAIC without any grouping of data nodes ## k-fold cross-validation -The `runCrossValidate` function in NIMBLE performs k-fold cross-validation on a `nimbleModel` fit via MCMC. More information can be found by calling `help(runCrossValidate)`. +The `runCrossValidate` function in NIMBLE performs k-fold cross-validation on a `nimbleModel` fit via MCMC. More information can be found by calling `help(runCrossValidate)`. ## Variable selection using Reversible Jump MCMC {#sec:rjmcmc} @@ -1092,20 +1092,20 @@ A common method for Bayesian variable selection in regression-style problems is A different view of the problem is that the different combinations of coefficients represent models of different dimensions. Reversible Jump MCMC [@Green_1995] (RJMCMC) is a general framework for MCMC simulation in which the dimension of the parameter space can vary between iterates of the Markov chain. The reversible jump sampler can be viewed as an extension of the Metropolis-Hastings algorithm onto more general state spaces. NIMBLE provides an implementation of RJMCMC for variable selection that requires the user to write the largest model of interest and supports use of indicator variables but formally uses RJMCMC sampling for better mixing and efficiency. When a coefficient is not in the model (or its indicator is 0), it will not be sampled, and it will therefore not follow its prior in that case. In technical detail, given two models $M_1$ and $M_2$ of possibly different dimensions, the core idea of RJMCMC is to remove the difference in the dimensions of models $M_1$ and $M_2$ by supplementing the corresponding parameters $\boldsymbol{\theta_1}$ and $\boldsymbol{\theta_2}$ with auxiliary variables $\boldsymbol{u}_{1 \rightarrow 2}$ and $\boldsymbol{u}_{2 \rightarrow 1}$ such that $(\boldsymbol{\theta_1}, \boldsymbol{u}_{1 \rightarrow 2})$ and $(\boldsymbol{\theta_2}, \boldsymbol{u}_{2 \rightarrow 1})$ -are in bijection, $(\boldsymbol{\theta_2}, \boldsymbol{u}_{2 \rightarrow 1}) = \Psi(\boldsymbol{\theta_1}, \boldsymbol{u}_{1 \rightarrow 2})$. The corresponding Metropolis-Hastings acceptance probability is generalized accounting for the proposal density for the auxiliary variables. +are in bijection, $(\boldsymbol{\theta_2}, \boldsymbol{u}_{2 \rightarrow 1}) = \Psi(\boldsymbol{\theta_1}, \boldsymbol{u}_{1 \rightarrow 2})$. The corresponding Metropolis-Hastings acceptance probability is generalized accounting for the proposal density for the auxiliary variables. NIMBLE implements RJMCMC for variable selection using a univariate normal distribution centered on 0 (or some fixed value) as the proposal density for parameters being added to the model. Two ways to set up models for RJMCMC are supported, which differ by whether the inclusion probabilities for each parameter are assumed known or must be written in the model: - If the inclusion probabilities are assumed known, then RJMCMC may be used with regular model code, i.e. model code written without heed to variable selection. - If the inclusion probability is a parameter in the model, perhaps with its own prior, then RJMCMC requires that indicator variables be written in the model code. The indicator variables will be sampled using RJMCMC, but otherwise they are like any other nodes in the model. - + The steps to set up variable selection using RJMCMC are: 1. Write the model with indicator variables if needed. 2. Configure the MCMC as usual with `configureMCMC`. 3. Modify the resulting MCMC configuration object with `configureRJ`. -The `configureRJ` function modifies the MCMC configuration to (1) assign samplers that turn on and off variable in the model and (2) modify the existing samplers for the regression coefficients to use 'toggled' versions. The toggled versions invoke the samplers only when corresponding variable is currently in the model. In the case where indicator variables are included in the model, sampling to turn on and off variables is done using `RJ_indicator` samplers. In the case where indicator variables are not included, sampling is done using `RJ_fixed_prior` samplers. +The `configureRJ` function modifies the MCMC configuration to (1) assign samplers that turn on and off variable in the model and (2) modify the existing samplers for the regression coefficients to use 'toggled' versions. The toggled versions invoke the samplers only when corresponding variable is currently in the model. In the case where indicator variables are included in the model, sampling to turn on and off variables is done using `RJ_indicator` samplers. In the case where indicator variables are not included, sampling is done using `RJ_fixed_prior` samplers. @@ -1121,7 +1121,7 @@ code <- nimbleCode({ beta0 ~ dnorm(0, sd = 100) beta1 ~ dnorm(0, sd = 100) beta2 ~ dnorm(0, sd = 100) - sigma ~ dunif(0, 100) + sigma ~ dunif(0, 100) z1 ~ dbern(psi) ## indicator variable associated with beta1 z2 ~ dbern(psi) ## indicator variable associated with beta2 psi ~ dbeta(1, 1) ## hyperprior on inclusion probability @@ -1142,7 +1142,7 @@ Y <- rnorm(N, 1.5 + 2 * x1, sd = 1) ## build the model and configure default MCMC RJexampleModel <- nimbleModel(code, constants = list(N = N), - data = list(Y = Y, x1 = x1, x2 = x2), + data = list(Y = Y, x1 = x1, x2 = x2), inits = list(beta0 = 0, beta1 = 0, beta2 = 0, sigma = sd(Y), z2 = 1, z1 = 1, psi = 0.5)) @@ -1158,13 +1158,13 @@ At this point we can modify the MCMC configuration object, `RJexampleConf`, to u ```{r indicators: configureRJ, results = "hide"} -configureRJ(conf = RJexampleConf, +configureRJ(conf = RJexampleConf, targetNodes = c("beta1", "beta2"), indicatorNodes = c('z1', 'z2'), control = list(mean = c(0, 0), scale = 2)) ``` -The `targetNodes` argument gives the coefficients (nodes) for which we want to do variable selection. The `indicatorNodes` gives the corresponding indicator nodes, ordered to match `targetNodes`. The `control` list gives the means and scale (standard deviation) for normal reversible jump proposals for `targetNodes`. The means will typically be 0 (by default `mean = 0` and `scale = 1`), but they could be anything. An optional `control` element called `fixedValue` can be provided in the non-indicator setting; this gives the value taken by nodes in `targetNodes` when they are out of the model (by default, `fixedValue` is 0). +The `targetNodes` argument gives the coefficients (nodes) for which we want to do variable selection. The `indicatorNodes` gives the corresponding indicator nodes, ordered to match `targetNodes`. The `control` list gives the means and scale (standard deviation) for normal reversible jump proposals for `targetNodes`. The means will typically be 0 (by default `mean = 0` and `scale = 1`), but they could be anything. An optional `control` element called `fixedValue` can be provided in the non-indicator setting; this gives the value taken by nodes in `targetNodes` when they are out of the model (by default, `fixedValue` is 0). All `control` elements can be scalars or vectors. A scalar values will be used for all `targetNodes`. A vector value must be of equal length as `targetNodes` and will be used in order. To use RJMCMC on a vector of coefficients with a corresponding vector of indicator variables, simply provide the variable names for `targetNodes` and `indicatorNodes`. For example, `targetNodes = "beta"` is equivalent to `targetNodes = c("beta[1]", "beta[2]")` and `indicatorNodes = "z"` is equivalent to `indicatorNodes = c("z[1]", "z[2]")`. Expansion of variable names into a vector of node names will occur as described in Section \@ref(sec:arbitr-coll-nodes). When using this method, *both* arguments must be provided as variable names to be expanded. @@ -1174,7 +1174,7 @@ Next we can see the result of `configureRJ` by printing the modified list of sam RJexampleConf$printSamplers() ``` -An `RJ_indicator` sampler was assigned to each of `z[1]` and `z[2]` in place of the `binary` sampler, while the samplers for `beta[1]`and `beta[2]` have been changed to `RJ_toggled` samplers. The latter samplers contain the original samplers, in this case `conjugate_dnorm_dnorm` samplers, and use them only when the corresponding indicator variable is equal to $1$ (i.e., when the coefficient is in the model). +An `RJ_indicator` sampler was assigned to each of `z[1]` and `z[2]` in place of the `binary` sampler, while the samplers for `beta[1]`and `beta[2]` have been changed to `RJ_toggled` samplers. The latter samplers contain the original samplers, in this case `conjugate_dnorm_dnorm` samplers, and use them only when the corresponding indicator variable is equal to $1$ (i.e., when the coefficient is in the model). Notice that the order of the samplers has changed, since `configureRJ` calls `removeSampler` for nodes in `targetNodes` and `indicatorNodes`, and subsequently `addSampler`, which appends the sampler to the end of current sampler list. Order can be modified by using `setSamplers`. @@ -1200,7 +1200,7 @@ code <- nimbleCode({ ## build the model RJexampleModel2 <- nimbleModel(code, constants = list(N = N), - data = list(Y = Y, x1 = x1, x2 = x2), + data = list(Y = Y, x1 = x1, x2 = x2), inits = list(beta0 = 0, beta1 = 0, beta2 = 0, sigma = sd(Y))) @@ -1212,11 +1212,11 @@ RJexampleConf2 <- configureMCMC(RJexampleModel2) RJexampleConf2$printSamplers() ``` -In this case, since there are no indicator variables, we need to pass to `configureRJ` the prior inclusion probabilities for each node in `targetNodes`, by specifying either one common value or a vector of values for the argument `priorProb`. This case does not allow for a stochastic prior. +In this case, since there are no indicator variables, we need to pass to `configureRJ` the prior inclusion probabilities for each node in `targetNodes`, by specifying either one common value or a vector of values for the argument `priorProb`. This case does not allow for a stochastic prior. ```{r no indicators: configureRJ, results = "hide"} -configureRJ(conf = RJexampleConf2, +configureRJ(conf = RJexampleConf2, targetNodes = c("beta1", "beta2"), priorProb = 0.5, control = list(mean = 0, scale = 2, fixedValue = c(1.5, 0))) @@ -1227,7 +1227,7 @@ configureRJ(conf = RJexampleConf2, RJexampleConf2$printSamplers() ``` -Since there are no indicator variables, the `RJ_fixed_prior` sampler is assigned directly to each of `beta[1]` and `beta[2]` along with the `RJ_toggled` sampler. The former sets a coefficient to its `fixedValue` when it is out of the model. The latter invokes the regular sampler for the coefficient only if it is in the model at a given iteration. +Since there are no indicator variables, the `RJ_fixed_prior` sampler is assigned directly to each of `beta[1]` and `beta[2]` along with the `RJ_toggled` sampler. The former sets a coefficient to its `fixedValue` when it is out of the model. The latter invokes the regular sampler for the coefficient only if it is in the model at a given iteration. If `fixedValue` is given when using `indicatorNodes` the values provided in `fixedValue` are ignored. However the same behavior can be obtained in this situation, using a different model specification. For example, the model in \@ref(sec:rjmcmc-indicator) can be modified to have `beta1` equal to $1.5$ when not in the model as follows: @@ -1259,7 +1259,7 @@ In the remainder of this section, we provide additional explanation of conjugate ### Conjugate ('Gibbs') samplers -By default, `configureMCMC()` and `buildMCMC()` will assign conjugate samplers to all nodes satisfying a conjugate relationship, unless the option `useConjugacy = FALSE` is specified. +By default, `configureMCMC()` and `buildMCMC()` will assign conjugate samplers to all nodes satisfying a conjugate relationship, unless the option `useConjugacy = FALSE` is specified. The current release of NIMBLE supports conjugate sampling of the relationships listed in Table \@ref(tab:conjugacy)^[NIMBLE's internal definitions of these relationships can be viewed with `nimble:::conjugacyRelationshipsInputList`.]. @@ -1269,7 +1269,7 @@ The current release of NIMBLE supports conjugate sampling of the relationships l Conjugate sampler functions may (optionally) dynamically check that their own posterior likelihood calculations are correct. If incorrect, a warning is issued. However, this functionality will roughly double the run-time required for conjugate sampling. By default, this option is disabled in NIMBLE. This option may be enabled with `nimbleOptions(verifyConjugatePosteriors = TRUE)`. -If one wants information about conjugate node relationships for other purposes, they can be obtained using the `checkConjugacy` method on a model. This returns a named list describing all conjugate relationships. The `checkConjugacy` method can also accept a character vector argument specifying a subset of node names to check for conjugacy. +If one wants information about conjugate node relationships for other purposes, they can be obtained using the `checkConjugacy` method on a model. This returns a named list describing all conjugate relationships. The `checkConjugacy` method can also accept a character vector argument specifying a subset of node names to check for conjugacy. @@ -1279,7 +1279,7 @@ If one wants information about conjugate node relationships for other purposes, - + @@ -1291,7 +1291,7 @@ If one wants information about conjugate node relationships for other purposes, - + @@ -1330,11 +1330,11 @@ The `nimbleHMC` package provides two distinct HMC samplers, which implement two The samplers provided in `nimbleHMC` can be used in `nimble`'s general MCMC system, and may be used in combination with other samplers provided with `nimble` operating on different (or overlapping) parts of the model. For convenience, `nimbleHMC` provides several different ways to set up and run HMC which allow different degrees of control: -- `nimbleHMC` is analogous to `nimbleMCMC`: It accpets model `code`, `data`, `inits` and `constants` arguments, and performs all steps from building the model to running one or more MCMC chains. Use the `nimbleHMC` function if you want one call to apply `NUTS` sampling to all continuous-valued model dimensions, and `nimble`'s default samplers for all discrete-valued dimensions. +- `nimbleHMC` is analogous to `nimbleMCMC`: It accepts model `code`, `data`, `inits` and `constants` arguments, and performs all steps from building the model to running one or more MCMC chains. Use the `nimbleHMC` function if you want one call to apply `NUTS` sampling to all continuous-valued model dimensions, and `nimble`'s default samplers for all discrete-valued dimensions. - `buildHMC` is analogous to `buildMCMC`: It accepts a model object as an argument, and builds an MCMC algorithm which applies `NUTS` sampling to all continuous-valued model dimensions, and `nimble`'s default samplers for all discrete-valued dimensions. Use `buildHMC` if you have created a model object, and want to automatically build the default HMC algorithm for your model. - `configureHMC` is analogous to `configureMCMC`: It accepts a model object as an argument, and sets up a default MCMC configuration containing one `NUTS` sampler operating on all continuous-valued model dimensions, and `nimble`'s default samplers for all discrete-valued dimensions. You can also specify which `nodes` should have samplers assigned. You can then modify the MCMC configuration (use `addSampler` and other methods) to further customize it before building and running the MCMC. Use `configureHMC` if you are familiar with MCMC configuration objects, and may want to further customize your MCMC using HMC or other samplers or monitors before building the MCMC. -- `addHMC` manages use of `addSampler` to specifically add an HMC sampler for a chosen set of nodes to an exisiting MCMC configuration. It optionally removes other samplers operating on the specified nodes. Use `addHMC` to add an HMC sampler operating one or more target nodes to an exisiting MCMC configuration object. -- For the most fine-grained control, you can use the `addSampler` method speciying the the argument `type="NUTS"` or `type="NUTS_classic"` to add HMC sampler(s) to an MCMC configuration object. +- `addHMC` manages use of `addSampler` to specifically add an HMC sampler for a chosen set of nodes to an existing MCMC configuration. It optionally removes other samplers operating on the specified nodes. Use `addHMC` to add an HMC sampler operating one or more target nodes to an existing MCMC configuration object. +- For the most fine-grained control, you can use the `addSampler` method specifying the argument `type="NUTS"` or `type="NUTS_classic"` to add HMC sampler(s) to an MCMC configuration object. As of version 1.1.0, we now support the following new functionality for AD in models: @@ -1353,11 +1353,11 @@ We'll use a Poisson Generalized Linear Mixed Model (GLMM) as a simple example mo ```{r eval=TRUE} model_code <- nimbleCode({ - # priors + # priors intercept ~ dnorm(0, sd = 100) beta ~ dnorm(0, sd = 100) sigma ~ dunif(0, 10) - # random effects and data + # random effects and data for(i in 1:10) { # random effects ran_eff[i] ~ dnorm(0, sd = sigma) @@ -1422,10 +1422,10 @@ The `nimbleSMC` package provides samplers that perform particle MCMC, primarily ### Customized log-likelihood evaluations: *RW_llFunction sampler* - + Sometimes it is useful to control the log-likelihood calculations used for an MCMC updater instead of simply using the model. For example, one could use a sampler with a log-likelihood that analytically (or numerically) integrates over latent model nodes. Or one could use a sampler with a log-likelihood that comes from a stochastic approximation such as a particle filter (see below), allowing composition of a particle MCMC (PMCMC) algorithm [@Andrieu_Doucet_Holenstein_2010]. The `RW_llFunction` sampler handles this by using a Metropolis-Hastings algorithm with a normal proposal distribution and a user-provided log-likelihood function. To allow compiled execution, the log-likelihood function must be provided as a specialized instance of a nimbleFunction. The log-likelihood function may use the same model as the MCMC as a setup argument (as does the example below), but if so the state of the model should be unchanged during execution of the function (or you must understand the implications otherwise). -The `RW_llFunction` sampler can be customized using the `control` list argument to set the initial proposal distribution scale and the adaptive properties for the Metropolis-Hastings sampling. In addition, the `control` list argument must contain a named `llFunction` element. This is the specialized nimbleFunction that calculates the log-likelihood; it must accept no arguments and return a scalar double number. The return value must be the total log-likelihood of all stochastic dependents of the `target` nodes -- and, if `includesTarget = TRUE`, of the target node(s) themselves -- or whatever surrogate is being used for the total log-likelihood. This is a required `control` list element with no default. See `help(samplers)` for details. +The `RW_llFunction` sampler can be customized using the `control` list argument to set the initial proposal distribution scale and the adaptive properties for the Metropolis-Hastings sampling. In addition, the `control` list argument must contain a named `llFunction` element. This is the specialized nimbleFunction that calculates the log-likelihood; it must accept no arguments and return a scalar double number. The return value must be the total log-likelihood of all stochastic dependents of the `target` nodes -- and, if `includesTarget = TRUE`, of the target node(s) themselves -- or whatever surrogate is being used for the total log-likelihood. This is a required `control` list element with no default. See `help(samplers)` for details. Here is a complete example: @@ -1462,7 +1462,7 @@ Rmcmc <- buildMCMC(mcmcConf) Note that we need to return `ll[1]` and not just `ll` because there are no scalar variables in compiled models. Hence `y` and other variables, and therefore `ll`, are of dimension 1 (i.e., vectors / one-dimensional arrays), so we need to specify the first element in order to have the return type be a scalar. - @@ -1483,7 +1483,7 @@ litters_code <- nimbleCode({ b[i] ~ dgamma(1, .001) for (j in 1:N) { r[i,j] ~ dbin(p[i,j], n[i,j]) - p[i,j] ~ dbeta(a[i], b[i]) + p[i,j] ~ dbeta(a[i], b[i]) } mu[i] <- a[i] / (a[i] + b[i]) theta[i] <- 1 / (a[i] + b[i]) @@ -1494,12 +1494,12 @@ litters_code <- nimbleCode({ constants <- list(G = 2, N = 16, n = matrix(c(13, 12, 12, 11, 9, 10, 9, 9, 8, 11, 8, 10, 13, - 10, 12, 9, 10, 9, 10, 5, 9, 9, 13, 7, 5, 10, 7, 6, + 10, 12, 9, 10, 9, 10, 5, 9, 9, 13, 7, 5, 10, 7, 6, 10, 10, 10, 7), nrow = 2)) # list specifying model data data <- list(r = matrix(c(13, 12, 12, 11, 9, 10, 9, 9, 8, 10, 8, 9, 12, 9, - 11, 8, 9, 8, 9, 4, 8, 7, 11, 4, 4, 5 , 5, 3, 7, 3, + 11, 8, 9, 8, 9, 4, 8, 7, 11, 4, 4, 5 , 5, 3, 7, 3, 7, 0), nrow = 2)) # list specifying initial values @@ -1601,5 +1601,3 @@ It is possible to run multiple chains in parallel using standard R parallelizati Thus, in your parallel loop or lapply-style statement, you should run `nimbleModel` and all subsequent calls to create and compile the model and MCMC algorithm within the parallelized block of code, once for each MCMC chain being run in parallel. For a worked example, please see [the parallelization example on our webpage](https://r-nimble.org/nimbleExamples/parallelizing_NIMBLE.html). We also now provide a [wrapper function, `nimbleParallelMCMC`](https://github.com/nimble-dev/nimble-demos/blob/master/parallelization/parallelizingNimble.R), that makes it easier to run chains in parallel. - - diff --git a/UserManual/src/chapter_MoreIntroduction.Rmd b/UserManual/src/chapter_MoreIntroduction.Rmd index e8846abc0..4a139895d 100644 --- a/UserManual/src/chapter_MoreIntroduction.Rmd +++ b/UserManual/src/chapter_MoreIntroduction.Rmd @@ -5,7 +5,7 @@ # More introduction {#cha-more-introduction} Now that we have shown a brief example, we will introduce more about -the concepts and design of NIMBLE. +the concepts and design of NIMBLE. One of the most important concepts behind NIMBLE is to allow a combination of high-level processing in R and low-level processing in @@ -22,7 +22,7 @@ We adopted the BUGS language, and we have extended it to make it more flexible. The BUGS language became widely used in WinBUGS, then in OpenBUGS and JAGS. These systems all provide automatically-generated MCMC algorithms, but we have adopted only the language for describing -models, not their systems for generating MCMCs. +models, not their systems for generating MCMCs. NIMBLE extends BUGS by: @@ -35,14 +35,14 @@ distributions and use them in BUGS models; By supporting new functions and distributions, NIMBLE makes BUGS an extensible language, which is a major departure from previous -packages that implement BUGS. +packages that implement BUGS. We adopted BUGS because it has been so successful, with over 30,000 users by the time they stopped counting [@Lunn_Spiegelhalter_Thomas_Best_2009]. Many papers and books provide BUGS code as a way to document their statistical models. We describe NIMBLE's version of BUGS later. The web sites for WinBUGS, -OpenBUGS and JAGS provide other useful documntation on writing models +OpenBUGS and JAGS provide other useful documentation on writing models in BUGS. For the most part, if you have BUGS code, you can try NIMBLE. @@ -159,13 +159,13 @@ In Version `r nimbleVersion`, the NIMBLE algorithm library includes: algorithm. 1. A variety of basic functions that can be used as programming tools for larger algorithms. These include: - + a. A likelihood function for arbitrary parts of any model. a. Functions to simulate one or many sets of values for arbitrary parts of any model. a. Functions to calculate the summed log probability (density) for one or many sets of values for arbitrary parts of any model along with stochastic dependencies in the model structure. - + - + More about the NIMBLE algorithm library is in Chapter \@ref(cha-algos-provided). diff --git a/UserManual/src/chapter_OtherAlgorithms.Rmd b/UserManual/src/chapter_OtherAlgorithms.Rmd index a89460d8e..b9103b396 100644 --- a/UserManual/src/chapter_OtherAlgorithms.Rmd +++ b/UserManual/src/chapter_OtherAlgorithms.Rmd @@ -4,7 +4,7 @@ ```{r, echo=FALSE} require(nimble) require(nimbleSMC, quietly = TRUE, warn.conflicts = FALSE) -``` +``` # Particle Filters, PMCMC, and MCEM {#cha-algos-provided} @@ -16,7 +16,7 @@ The NIMBLE algorithm library includes a growing set of non-MCMC algorithms, incl - nested approximation of Bayesian models using INLA-like methods. In this chapter we discuss sequential Monte Carlo methods and MCEM and in the next chapter we discuss methods related to Laplace approximation. - + ## Particle filters / sequential Monte Carlo and iterated filtering As of Version 0.10.0 of NIMBLE, all of NIMBLE's sequential Monte Carlo/particle filtering functionality lives in the `nimbleSMC` package, described in @Michaud_etal_2021. Please load this package before trying to use these algorithms. @@ -24,13 +24,13 @@ As of Version 0.10.0 of NIMBLE, all of NIMBLE's sequential Monte Carlo/particle ### Filtering algorithms NIMBLE includes algorithms for four different types of sequential Monte Carlo (also known as particle filters), which can be used to sample from the latent states and approximate the log likelihood of a state-space model. These include the bootstrap filter, the auxiliary particle filter, the Liu-West filter, and the ensemble Kalman filter. The iterated filtering version 2 (IF2) is a related method for maximum-likelihood estimation. Each of these is built with the eponymous functions `buildBootstrapFilter`, `buildAuxiliaryFilter`, `buildLiuWestFilter`, `buildEnsembleKF`, and `buildIteratedFilter2`. Each method requires setup arguments `model` and `nodes`; the latter should be a character vector specifying latent model nodes. In addition, each method can be customized using a `control` list argument. Details on the control options and specifics of the algorithms can be found in the help pages for the functions. - + Once built, each filter can be run by specifying the number of particles. Each filter has a modelValues object named `mvEWSamples` that is populated with equally-weighted samples from the posterior distribution of the latent states (and in the case of the Liu-West filter, the posterior distribution of the top level parameters as well) as the filter is run. The bootstrap, auxiliary, and Liu-West filters, as well as the IF2 method, also have another modelValues object, `mvWSamples`. This has unequally-weighted samples from the posterior distribution of the latent states, along with weights for each particle. In addition, the bootstrap and auxiliary particle filters return estimates of the log-likelihood of the given state-space model. - - We first create a linear state-space model to use as an example for our particle filter algorithms. - + + We first create a linear state-space model to use as an example for our particle filter algorithms. + ```{r, particle_Filter_Chunk, results = "hide"} -# Building a simple linear state-space model. +# Building a simple linear state-space model. # x is latent space, y is observed data timeModelCode <- nimbleCode({ x[1] ~ dnorm(mu_0, 1) @@ -39,7 +39,7 @@ timeModelCode <- nimbleCode({ x[i] ~ dnorm(x[i-1] * a + b, 1) y[i] ~ dnorm(x[i] * c, 1) } - + a ~ dunif(0, 1) b ~ dnorm(0, 1) c ~ dnorm(1,1) @@ -57,7 +57,7 @@ for(i in 2:t){ } # build the model -rTimeModel <- nimbleModel(timeModelCode, constants = list(t = t), +rTimeModel <- nimbleModel(timeModelCode, constants = list(t = t), data <- list(y = y), check = FALSE ) # Set parameter values and compile the model @@ -70,22 +70,22 @@ cTimeModel <- compileNimble(rTimeModel) ``` #### Bootstrap filter - + Here is an example of building and running the bootstrap filter. - + ```{r, boot_Filter_Chunk} # Build bootstrap filter -rBootF <- buildBootstrapFilter(rTimeModel, "x", - control = list(thresh = 0.8, saveAll = TRUE, +rBootF <- buildBootstrapFilter(rTimeModel, "x", + control = list(thresh = 0.8, saveAll = TRUE, smoothing = FALSE)) -# Compile filter +# Compile filter cBootF <- compileNimble(rBootF,project = rTimeModel) # Set number of particles parNum <- 5000 # Run bootstrap filter, which returns estimate of model log-likelihood bootLLEst <- cBootF$run(parNum) -# The bootstrap filter can also return an estimate of the effective +# The bootstrap filter can also return an estimate of the effective # sample size (ESS) at each time point bootESS <- cBootF$returnESS() ``` @@ -93,38 +93,38 @@ bootESS <- cBootF$returnESS() #### Auxiliary particle filter Next, we provide an example of building and running the auxiliary particle filter. Note that a filter cannot be built on a model that already has a filter specialized to it, so we create a new copy of our state space model first. - + ```{r, aux_Filter_Chunk, results = "hide"} # Copy our state-space model for use with the auxiliary filter auxTimeModel <- rTimeModel$newModel(replicate = TRUE) compileNimble(auxTimeModel) # Build auxiliary filter -rAuxF <- buildAuxiliaryFilter(auxTimeModel, "x", +rAuxF <- buildAuxiliaryFilter(auxTimeModel, "x", control = list(thresh = 0.5, saveAll = TRUE)) -# Compile filter +# Compile filter cAuxF <- compileNimble(rAuxF,project = auxTimeModel) # Run auxiliary filter, which returns estimate of model log-likelihood auxLLEst <- cAuxF$run(parNum) -# The auxiliary filter can also return an estimate of the effective +# The auxiliary filter can also return an estimate of the effective # sample size (ESS) at each time point auxESS <- cAuxF$returnESS() ``` #### Liu and West filter -Now we give an example of building and running the Liu and West filter, which can sample from the posterior distribution of top-level parameters as well as latent states. Note that the Liu-West filter ofen performs poorly and is provided primarily for didactic purposes. The Liu and West filter accepts an additional `params` argument, specifying the top-level parameters to be sampled. +Now we give an example of building and running the Liu and West filter, which can sample from the posterior distribution of top-level parameters as well as latent states. Note that the Liu-West filter ofen performs poorly and is provided primarily for didactic purposes. The Liu and West filter accepts an additional `params` argument, specifying the top-level parameters to be sampled. ```{r, lw_Filter_Chunk, results = "hide"} # Copy model LWTimeModel <- rTimeModel$newModel(replicate = TRUE) compileNimble(LWTimeModel) -# Build Liu-West filter, also +# Build Liu-West filter, also # specifying which top level parameters to estimate rLWF <- buildLiuWestFilter(LWTimeModel, "x", params = c("a", "b", "c"), - control = list(saveAll = FALSE)) -# Compile filter + control = list(saveAll = FALSE)) +# Compile filter cLWF <- compileNimble(rLWF,project = LWTimeModel) # Run Liu-West filter cLWF$run(parNum) @@ -132,7 +132,7 @@ cLWF$run(parNum) #### Ensemble Kalman filter -Next we give an example of building and running the ensemble Kalman filter, which can sample from the posterior distribution of latent states. +Next we give an example of building and running the ensemble Kalman filter, which can sample from the posterior distribution of latent states. ```{r, ENKF_Filter_Chunk, results = "hide"} # Copy model @@ -140,14 +140,14 @@ ENKFTimeModel <- rTimeModel$newModel(replicate = TRUE) compileNimble(ENKFTimeModel) # Build and compile ensemble Kalman filter rENKF <- buildEnsembleKF(ENKFTimeModel, "x", - control = list(saveAll = FALSE)) + control = list(saveAll = FALSE)) cENKF <- compileNimble(rENKF,project = ENKFTimeModel) # Run ensemble Kalman filter cENKF$run(parNum) ``` - + Once each filter has been run, we can extract samples from the posterior distribution of our latent states as follows: - + ```{r, particle_Filter_Samples, eval=FALSE} # Equally-weighted samples (available from all filters) bootEWSamp <- as.matrix(cBootF$mvEWSamples) # alternative: as.list @@ -155,14 +155,14 @@ auxEWSamp <- as.matrix(cAuxF$mvEWSamples) LWFEWSamp <- as.matrix(cLWF$mvEWSamples) ENKFEWSamp <- as.matrix(cENKF$mvEWSamples) -# Unequally-weighted samples, along with weights (available +# Unequally-weighted samples, along with weights (available # from bootstrap, auxiliary, and Liu and West filters) bootWSamp <- as.matrix(cBootF$mvWSamples, "x") bootWts <- as.matrix(cBootF$mvWSamples, "wts") auxWSamp <- as.matrix(xAuxF$mvWSamples, "x") auxWts <- as.matrix(cAuxF$mvWSamples, "wts") -# Liu and West filter also returns samples +# Liu and West filter also returns samples # from posterior distribution of top-level parameters: aEWSamp <- as.matrix(cLWF$mvEWSamples, "a") ``` @@ -179,7 +179,7 @@ library(FKF) flowCode <- nimbleCode({ for(t in 1:n) y[t] ~ dnorm(x[t], sd = sigmaMeasurements) - x[1] ~ dnorm(x0, sd = sigmaInnovations) + x[1] ~ dnorm(x0, sd = sigmaInnovations) for(t in 2:n) x[t] ~ dnorm((t-1==28)*meanShift1899 + x[t-1], sd = sigmaInnovations) logSigmaInnovations ~ dnorm(0, sd = 100) @@ -245,7 +245,7 @@ fkfResult <- fkf(HHt = matrix(exp(2*est[1])), fkfResult$logLik ``` - + ### Particle MCMC (PMCMC) {#sec:particle-mcmc} Particle MCMC (PMCMC) is a method that uses MCMC for top-level model @@ -287,10 +287,10 @@ timeConf$addSampler(target = c("a", "b", "c", "mu_0"), type = "RW_PF_block", control = list(propCov= diag(4), adaptScaleOnly = FALSE, latents = "x", pfOptimizeNparticles = TRUE)) ``` - + The `type = "RW_PF"` and `type = "RW_PF_block"` samplers default to using a bootstrap filter. The adaptation control -parameters `adaptive`, `adaptInterval`, and +parameters `adaptive`, `adaptInterval`, and `adaptScaleOnly` work in the same way as for an `RW` and `RW*block` samplers. However, it is not clear if the same approach to adaptation works well for PMCMC, so one should consider turning off adaptation and using @@ -304,7 +304,7 @@ list entry. The `pfType` entry can be set either to user-defined nimbleFunction that returns a likelihood approximation. Any user-defined filtering nimbleFunction named in the `pfType` -control list entry must satsify the following: +control list entry must satisfy the following: 1. The nimbleFunction must be the result of a call to @@ -321,7 +321,7 @@ control list entry must satsify the following: `control` list can be used to pass in any additional information or arguments that the custom filter may require. 1. The `nimbleFunction` must have a `run` function that - accepts a single integer arugment (the number of particles to use), + accepts a single integer argument (the number of particles to use), and returns a scalar double (the log-likelihood estimate). 1. The `nimbleFunction` must define, in setup code, a `modelValues` object named `mvEWSamples` that is used to @@ -344,13 +344,13 @@ control list entry must satsify the following: conditional distributions cannot be directly sampled from) at each iteration. MCEM can be slow but is also a general workhorse method that is applicable to a wide range of problems. - + `buildMCEM` provides an (optionally) ascent-based MCEM algorithm based on @Caffo_etal_2005. `buildMCEM` was re-written for nimble version 1.2.0 and is not compatible with previous versions. For the E (expectation) step, `buildMCEM` uses an MCMC for the latent states given the data and current parameters. For the M (maximization) step, `buildMCEM` uses `optim`. - + The ascent-based feature uses an estimate of the standard error of the move in one iteration to determine if that move is clearly uphill or is swamped by noise due to using a Monte Carlo (MCMC) sample. In the latter case, the MCMC @@ -359,7 +359,7 @@ control list entry must satsify the following: confidently less than a tolerance level taken to indicate that the algorithm is at the MLE. The associated tuning parameters, along with details on control over each step, are described at `help(buildMCEM)`. - + The maximization part of each MCEM iteration will use gradients from automatic differentiation if they are supported for the model and `buildDerivs=TRUE` when `nimbleModel` was called to create the model. If any parameters have @@ -369,12 +369,12 @@ control list entry must satsify the following: Additionally, the MCEM algorithm can provide an estimate of the asymptotic covariance matrix of the parameters. See \@ref(sec:estimate-mcem-cov). - + We will revisit the *pump* example to illustrate the use of NIMBLE's MCEM algorithm. ```{r, echo=FALSE} -pumpCode <- nimbleCode({ +pumpCode <- nimbleCode({ for (i in 1:N){ theta[i] ~ dgamma(alpha,beta) lambda[i] <- theta[i]*t[i] @@ -392,7 +392,7 @@ 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)) -``` +``` ```{r, build-MCEM, results = "hide", eval = runMCEMs} @@ -415,7 +415,7 @@ model provided cannot be part of another MCMC sampler. ```{r, run-MCEM, eval = runMCEMs} pumpMLE <- CpumpMCEM$findMLE(initM = 1000) pumpMLE$par -``` +``` For this model, we can check the result by analytically integrating over the latent nodes (possible for this model but often not feasible), which gives @@ -441,6 +441,4 @@ Continuing the above example, here is the covariance matrix: ```{r, run-MCEM-cov, eval = runMCEMs} pumpCov <- CpumpMCEM$vcov(pumpMLE$par) pumpCov -``` - - +``` diff --git a/UserManual/src/chapter_RCfunctions.Rmd b/UserManual/src/chapter_RCfunctions.Rmd index a6ad9bc81..619d6f168 100644 --- a/UserManual/src/chapter_RCfunctions.Rmd +++ b/UserManual/src/chapter_RCfunctions.Rmd @@ -5,11 +5,11 @@ ```{r, echo=FALSE} require(nimble) -``` +``` ```{r, echo=FALSE} read_chunk(file.path('chunks', 'programmingWithModels_chunks.R')) # one can put code chunks here if one wants -``` +``` ## Introduction to simple nimbleFunctions {#sec:RC-intro} @@ -31,13 +31,13 @@ extension of BUGS (Chapter \@ref(cha-user-defined)). least squares estimation of linear regression parameters^[Of course, in general, explicitly calculating the inverse is not the recommended numerical recipe for least squares.]: - + ```{r, nf-RCfun, message=FALSE} -``` +``` In this example, we fit a linear model for 100 random response values (`y`) to four columns of randomly generated explanatory variables (`X`). We ran the -nimbleFunction +nimbleFunction `solveLeastSquares` uncompiled, natively in R, allowing testing and debugging (Section \@ref(sec:debugging)). Then we compiled it and showed that the compiled version does the same thing, but faster^[On the machine this is @@ -73,7 +73,7 @@ The next sections cover each of these topics in turn. Often, R help pages are available for NIMBLE's versions of R functions using the prefix 'nim' and capitalizing the next letter. For -example, help on NIMBLE's version of `numeric` can be found by `help(nimNumeric)`. In some cases help is found directly using the name of the function as it appears in R. +example, help on NIMBLE's version of `numeric` can be found by `help(nimNumeric)`. In some cases help is found directly using the name of the function as it appears in R. ### Basic operations @@ -87,22 +87,22 @@ In addition to the above functions, we provide functions `any_na()` and `any_nan Other R functions with numeric arguments and return value can be called during compiled execution by wrapping them as a `nimbleRcall` (see Section \@ref(sec:calling-R-code)). -Next we cover some of these functions in more detail. +Next we cover some of these functions in more detail. #### *numeric*, *integer*, *logical*, *matrix* and *array* `numeric`, `integer`, or `logical` will create a 1-dimensional vector of -floating-point (or 'double' [precision]), integer, or logical values, respectively. The `length` argument specifies the vector -length (default 0), and the `value` argument specifies the initial +floating-point (or 'double' [precision]), integer, or logical values, respectively. The `length` argument specifies the vector +length (default 0), and the `value` argument specifies the initial value either as a scalar (used for all vector elements, with default - 0) or a vector. If a vector and its length is not equal to `length`, the remaining values will be filled by R-style recycling if the `recycle` argument is `TRUE` (which is the default). The `init` argument specifies + 0) or a vector. If a vector and its length is not equal to `length`, the remaining values will be filled by R-style recycling if the `recycle` argument is `TRUE` (which is the default). The `init` argument specifies whether or not to initialize the elements in compiled code (default `TRUE`). If first use of the variable does not rely on initial values, using `init = FALSE` will yield slightly more efficient performance. -`matrix` creates a 2-dimensional matrix object of either floating-point (if -`type = "double"`, the default), integer (if `type = "integer"`), or logical (if `type = "logical"`) values. -As in R, `nrow` and `ncol` arguments specify the number of rows and columns, respectively. +`matrix` creates a 2-dimensional matrix object of either floating-point (if +`type = "double"`, the default), integer (if `type = "integer"`), or logical (if `type = "logical"`) values. +As in R, `nrow` and `ncol` arguments specify the number of rows and columns, respectively. The `value` and `init` arguments are used in the same way as for `numeric`, `integer`, and `logical`. @@ -127,7 +127,7 @@ for changing the size of a numeric object, are given in Section \@ref(sec:chang- `dim` function for matrices or arrays and like R's `length` function for vectors. In other words, regardless of whether the number of dimensions is 1 or more, it returns a vector - of the sizes. + of the sizes. #### Deprecated creation of non-scalar objects using *declare* @@ -181,7 +181,7 @@ from nimbleFunctions (and in BUGS model code), but care is needed in the syntax, - Names of the distributions generally (but not always) match those of R, which sometimes differ from BUGS. See the list below. - - Supported parameterizations are also indicated in the list below. + - Supported parameterizations are also indicated in the list below. - For multivariate distributions (multivariate normal, Dirichlet, and Wishart), 'r' functions only return one random draw at a time, and the first argument must always be 1. - R's recycling rule (re-use of an argument as needed to accommodate @@ -195,7 +195,7 @@ As in R (and nimbleFunctions), arguments are matched by order or by name (`log`, `log.p`, `lower.tail`) can be used and have the same defaults. User-defined distributions for BUGS (Chapter \@ref(cha-user-defined)) can also be used from -nimbleFunctions. +nimbleFunctions. For standard distributions, we rely on R's regular help pages (e.g., `help(dgamma)`. For distributions unique to NIMBLE @@ -211,7 +211,7 @@ Supported distributions, listed by their 'd' function, include: - `dbeta(x, shape1, shape2, log)` - `dchisq(x, df, log)` - `ddexp(x, location, rate, log)` - - `ddexp(x, location, scale, log)` + - `ddexp(x, location, scale, log)` - `dexp(x, rate, log)` - `dexp_nimble(x, rate, log)` - `dexp_nimble(x, scale, log)` @@ -235,7 +235,7 @@ Supported distributions, listed by their 'd' function, include: - `dinvwish_chol(x, cholesky, df, scale_param, log)` -In the last four, `cholesky` stands for Cholesky decomposition of the relevant matrix; +In the last four, `cholesky` stands for Cholesky decomposition of the relevant matrix; `prec_param` is a logical indicating whether the Cholesky is of a precision matrix (TRUE) or covariance matrix (FALSE)^[For the multivariate t, these are more properly termed the 'inverse scale' and 'scale' matrices]; and `scale_param` is a logical indicating whether the Cholesky is of a scale matrix (TRUE) or an inverse scale matrix (FALSE). @@ -262,7 +262,7 @@ argument that will be displayed in an error message. `print`, or equivalently `nimPrint`, prints an arbitrary set of outputs in order and adds a newline character at the end. `cat` or `nimCat` is identical, -except without a newline at the end. +except without a newline at the end. ### Checking for user interrupts: *checkInterrupt* {#sec:check-user-interr} @@ -290,7 +290,7 @@ another one with the additional parameters. Finally, `nimOptim` requires a simple R `list`. To define the `control` parameter, create a default value with the `nimOptimDefaultControl` function, and set any desired fields. -See `help(nimOptim)` for details, including example usage. +See `help(nimOptim)` for details, including example usage. ### Integration: *integrate* and *nimIntegrate* @@ -422,7 +422,7 @@ When a variable is first created by assignment, its type is determined automatically by that assignment. For example, if `x` has not appeared before, then ```{r, eval=FALSE} x <- A %*% B # assume A and B are double matrices or vectors -``` +``` will create `x` to be a double matrix of the correct size (determined during execution). @@ -431,7 +431,7 @@ will create `x` to be a double matrix of the correct size Because NIMBLE is statically typed, you cannot use the same variable name for two objects of different types (including objects of different dimensions). -Suppose we have (implicitly) created `x` as a double matrix. +Suppose we have (implicitly) created `x` as a double matrix. If `x` is used subsequently, it can only be used as a double matrix. This is true even if it is assigned a new value, which will again set its size automatically but cannot @@ -441,7 +441,7 @@ change its type. x <- A %*% B # assume A and B are double matrices or vectors x <- nimMatrix(0, nrow = 5, ncol = 2) # OK: 'x' is still a double matrix x <- rnorm(10) # NOT OK: 'x' is a double vector -``` +``` #### When a numeric object needs to be created before being used {#sec:when-numeric-object} @@ -456,8 +456,8 @@ contents for specific indices. ```{r, exampleDeclare, eval=FALSE} x <- numeric(10) -for(i in 1:10) - x[i] <- foo(y[i]) +for(i in 1:10) + x[i] <- foo(y[i]) ``` #### How NIMBLE handles `NA` @@ -475,19 +475,19 @@ These issues arise because NIMBLE uses R's encoding of NA values in C(++) but us #### Changing the sizes of existing objects: *setSize* {#sec:chang-sizes-exist} `setSize` changes the size of an object, preserving its contents in -column-major order. +column-major order. ```{r, exampleNumeric, eval=FALSE} # Example of creating and resizing a floating-point vector -# myNumericVector will be of length 10, with all elements initialized to 2 -myNumericVector <- numeric(10, value = 2) +# myNumericVector will be of length 10, with all elements initialized to 2 +myNumericVector <- numeric(10, value = 2) # resize this numeric vector to be length 20; last 10 elements will be 0 setSize(myNumericVector, 20) ``` ```{r, exampleMatrix, eval=FALSE} # Example of creating a 1-by-10 matrix with values 1:10 and resizing it -myMatrix <- matrix(1:10, nrow = 1, ncol = 10) +myMatrix <- matrix(1:10, nrow = 1, ncol = 10) # resize this matrix to be a 10-by-10 matrix setSize(myMatrix, c(10, 10)) # The first column will have the 1:10 @@ -509,11 +509,11 @@ myfun <- nimbleFunction( d <- a + b # double vector f <- c(i) # integer vector }) -``` +``` In the line that creates `b`, the index range `i:i` is not evaluated until run time. Even though `i:i` will always evaluate -to simpy `i`, the compiler does not determine that. Since there is +to simply `i`, the compiler does not determine that. Since there is a vector index range provided, the result of `randomValues[i:i]` is determined to be a vector. The following line then creates `d` as a vector, because a vector plus a scalar returns a vector. Another @@ -525,7 +525,7 @@ illustrated in the last line. Consider the following code: ```{r, } myfun <- nimbleFunction( - run = function() { + run = function() { A <- matrix(value = rnorm(9), nrow = 3) B <- rnorm(3) Cmatrix <- A %*% B # double matrix, one column @@ -533,7 +533,7 @@ myfun <- nimbleFunction( Cmatrix <- (A %*% B)[,1] # error, vector assigned to matrix Cmatrix[,1] <- (A %*% B)[,1] # ok, if Cmatrix is large enough }) -``` +``` This creates a matrix `A`, a vector `B`, and matrix-multiplies them. The vector `B` is automatically treated as a one-column matrix in matrix algebra computations. The result of matrix multiplication is always a matrix, but @@ -570,11 +570,11 @@ unless it is known at compile time that `M1` has 1 column, in which case - `v1 %*% v2` promotes `v1` to a 1-row matrix and `v2` to a 1-column matrix, so the returned values is a 1x1 matrix with the inner product of `v1` and `v2`. If you want the inner product as a scalar, use -`inprod(v1, v2)`. +`inprod(v1, v2)`. - `asRow(v1)` explicitly promotes `v1` to a 1-row matrix. Therefore `v1 %*% asRow(v2)` gives the outer product of `v1` and `v2`. - - `asCol(v1)` explicitly promotes `v1` to a 1-column matrix. + - `asCol(v1)` explicitly promotes `v1` to a 1-column matrix. - The default promotion for a vector is to a 1-column matrix. Therefore, `v1 %*% t(v2)` is equivalent to `v1 %*% asRow(v2)` . - When indexing, dimensions with scalar indices will be dropped. @@ -585,8 +585,8 @@ product of `v1` and `v2`. If you want the inner product as a scalar, use it must already be correctly sized for the result. For example, `Y[5:10, 20:30] <- x` will not work -- and could crash your R session with a segmentation fault -- if Y is not - already at least 10x30 in size. This can be done by `setSize(Y, c(10, 30))`. See Section \@ref(sec:chang-sizes-exist) for more details. Note that non-indexed assignment to `Y`, such as `Y <- x`, will automatically set `Y` to the necessary size. - + already at least 10x30 in size. This can be done by `setSize(Y, c(10, 30))`. See Section \@ref(sec:chang-sizes-exist) for more details. Note that non-indexed assignment to `Y`, such as `Y <- x`, will automatically set `Y` to the necessary size. + Here are some examples to illustrate the above points, assuming `M2` is @@ -626,15 +626,15 @@ As illustrated in the example in Section \@ref(sec:RC-intro), the syntax for a type declaration is: `type(nDim, sizes)` - -where `type` is `double`, `integer`, `logical` or `character`. + +where `type` is `double`, `integer`, `logical` or `character`. (In more general nimbleFunction programming, a type can also be a nimbleList type, discussed in Section \@ref(sec:nimbleLists).) For example `run = function(x = double(1)) { ... `} sets the single argument of the run function to be a vector of numeric values of unknown size. For `type(nDim, sizes)`, `nDim` is the number of dimensions, with 0 indicating scalar and omission of `nDim` defaulting to a scalar. -`sizes` is an optional vector of fixed, known sizes. +`sizes` is an optional vector of fixed, known sizes. For example, `double(2, c(4, 5))` declares a $4 \times 5$ matrix. If sizes are omitted, they will be set either by assignment or by `setSize`. @@ -645,11 +645,11 @@ provided. For example, to provide 1.2 as a default: myfun <- nimbleFunction( run = function(x = double(0, default = 1.2)) { }) -``` +``` Functions with return values must have their return type explicitly declared using `returnType`, which can occur anywhere in the run code. For example `returnType(integer(2))` declares the return type to be a matrix of integers. A return type of `void()` means there is no return value, which is the default if no `returnType` statement -is included. +is included. Note that because all values in models are stored as doubles and because of some limitations in NIMBLE's automatic casting, non-scalar return values of user-defined distributions must be doubles. @@ -658,18 +658,18 @@ Note that because all values in models are stored as doubles and because of some Uncompiled nimbleFunctions pass arguments like R does, by copy. If `x` is passed as an argument to function `foo`, and `foo` modifies `x` internally, it is modifying its copy of `x`, not -the original `x` that was passed to it. +the original `x` that was passed to it. Compiled nimbleFunctions pass arguments to other compiled nimbleFunctions by reference (or pointer). This is very different. Now if `foo` modifies `x` internally, it is modifying the same `x` that was passed to it. This allows much faster execution but is obviously a fundamentally different behavior. - + Uncompiled execution of nimbleFunctions is primarily intended for debugging. However, debugging of how nimbleFunctions interact via arguments requires testing the compiled versions. - + ## Calling external compiled code {#sec:calling-external-code} If you have a function in your own compiled C or C++ code and an appropriate header file, you can generate a nimbleFunction that wraps access to that function, which can then be used in other nimbleFunctions. See `help(nimbleExternalCall)` for an example. This also contains an example of using an externally compiled function in the BUGS code of a model. diff --git a/UserManual/src/chapter_UsingModels.Rmd b/UserManual/src/chapter_UsingModels.Rmd index 4893e7e3a..8a5a08278 100644 --- a/UserManual/src/chapter_UsingModels.Rmd +++ b/UserManual/src/chapter_UsingModels.Rmd @@ -6,14 +6,14 @@ require(nimble) ```{r, echo=FALSE} read_chunk(file.path('chunks', 'usingBugsModels_chunks.R')) # one can put code chunks here if one wants -``` +``` # Working with NIMBLE models {#cha-using-models} Here we describe how one can get information about NIMBLE models and carry out operations on a model. While all of this functionality can -be used from R, its primary use occurs when writing nimbleFunctions (see Chapter \@ref(cha-progr-with-models)). +be used from R, its primary use occurs when writing nimbleFunctions (see Chapter \@ref(cha-progr-with-models)). Information about node types, distributions, and dimensions can be used to determine algorithm behavior in *setup* code of nimbleFunctions. Information about node or variable values or the parameter and bound values of a node would generally be used for algorithm @@ -30,18 +30,18 @@ Section \@ref(sec:nodes-and-variables) defines what we mean by variables and nod One can determine the variables in a model using `getVarNames` and the nodes in a model using `getNodeNames`, with optional arguments allowing you to select only certain types of nodes. We illustrate here with the pump model from Chapter \@ref(cha-lightning-intro). ```{r, getVarAndNodeNamesPump} -``` +``` -You can see one lifted node (see next section), `lifted_d1_over_beta`, involved in a reparameterization to NIMBLE's canonical parameterization of the gamma distribution for the `theta` nodes. +You can see one lifted node (see next section), `lifted_d1_over_beta`, involved in a reparameterization to NIMBLE's canonical parameterization of the gamma distribution for the `theta` nodes. We can determine the set of nodes contained in one or more nodes or variables using `expandNodeNames`, illustrated here for an example with multivariate nodes. The -`returnScalarComponents` argument also allows us to return all of the scalar elements of multivariate nodes. +`returnScalarComponents` argument also allows us to return all of the scalar elements of multivariate nodes. ```{r, expandNodeNames, message=FALSE} -``` +``` -As discussed in Section \@ref(sec:cdisdata), you can determine whether a node is flagged as data using `isData`. +As discussed in Section \@ref(sec:cdisdata), you can determine whether a node is flagged as data using `isData`. ### Understanding lifted nodes {#sec:introduced-nodes} @@ -73,7 +73,7 @@ generated. These include: uses `sd` as a parameter in the normal distribution. If many nodes use the same `tau`, only one new `sd` node will be created, so the computation `1/sqrt(tau)` will not be repeated - redundantly. + redundantly. @@ -109,7 +109,7 @@ information about nodes or variables. ```{r, distAPI-example, message=FALSE} code <- nimbleCode({ - for(i in 1:4) + for(i in 1:4) y[i] ~ dnorm(mu, sd = sigma) mu ~ T(dnorm(0, 5), -20, 20) sigma ~ dunif(0, 10) @@ -122,18 +122,18 @@ m$isDiscrete(c('y', 'mu', 'sigma')) m$isDeterm('mu') m$getDimension('mu') m$getDimension('mu', includeParams = TRUE) -``` +``` Note that any variables provided to these functions are expanded into their constituent node names, so the length of results may not be the same length as the input vector of node and variable names. However the order of the results should be preserved relative to the order of -the inputs, once the expansion is accounted for. +the inputs, once the expansion is accounted for. ### Getting information about a distribution One can also get generic information about a distribution based on the -name of the distribution using the function `getDistributionInfo`. +name of the distribution using the function `getDistributionInfo`. In particular, one can determine whether a distribution was provided by the user (`isUserDefined`), whether a distribution provides CDF and quantile functions (`pqDefined`), whether a distribution is a @@ -152,7 +152,7 @@ distribution followed by that node. The parameter does not have to be one of the parameters used when the node was declared. Alternative parameterization values can also be obtained. See Section \@ref(subsec:distributions) for available parameterizations. -(These can also be seen in `nimble:::distributionsInputList`.) +(These can also be seen in `nimble:::distributionsInputList`.) Here is an example: ```{r, getParamExample, message=FALSE} @@ -163,7 +163,7 @@ gammaModel <- nimbleModel( }), data = list(x = 2.4), inits = list(a = 1.2)) gammaModel$getParam('x', 'scale') gammaModel$getParam('x', 'rate') -``` +``` `getParam` is part of the NIMBLE language, so it can be used in run code of nimbleFunctions. @@ -172,14 +172,14 @@ gammaModel$getParam('x', 'rate') The function `getBound` provides access to the lower and upper bounds of the distribution for a node. In most cases these bounds will be fixed based on the distribution, but for the uniform distribution -the bounds are the parameters of the distribution, and when truncation is used (Section \@ref(subsec:trunc)), the bounds will be determined +the bounds are the parameters of the distribution, and when truncation is used (Section \@ref(subsec:trunc)), the bounds will be determined by the truncation. Like the functions described in the previous section, `getBound` can be used as global function taking a model as the first argument, or it can be used as a model member function. The next two arguments must be the name of one (stochastic) node and either `"lower"` or `"upper"` indicating whether the lower or upper bound is desired. For multivariate nodes the bound is a scalar that is the bound for all elements of the node, as we -do not handle truncation for multivariate nodes. +do not handle truncation for multivariate nodes. Here is an example: ```{r, getBoundExample, message=FALSE} @@ -196,7 +196,7 @@ exampleModel$b <- 3 exampleModel$calculate(exampleModel$getDependencies('b')) getBound(exampleModel, 'a', 'upper') exampleModel$getBound('b','lower') -``` +``` `getBound` is part of the NIMBLE language, so it can be used in run code of nimbleFunctions. In fact, we anticipate that most use of `getBound` will be for algorithms, such as for the reflection @@ -230,10 +230,10 @@ deterministic node. - **calculateDiff** `calculateDiff` is identical to `calculate`, but it returns the new log probability value minus the one that was previously stored. This is useful when one wants to - change the value or values of node(s) in the model (e.g., by setting a value or + change the value or values of node(s) in the model (e.g., by setting a value or `simulate`) and then determine the change in the log probability, such as needed for a - Metropolis-Hastings acceptance probability. + Metropolis-Hastings acceptance probability. Each of these functions is accessed as a member function of a model @@ -246,15 +246,15 @@ an example using `simulate`. #### Example: simulating arbitrary collections of nodes {#sec:arbitr-coll-nodes} -```{r, calcSimGLPdemos} -``` +```{r, calcSimGLPdemos} +``` The example illustrates a number of features: 1. `simulate(model, nodes)` is equivalent to `model$simulate(nodes)`. You can use either, but the latter is - encouraged and the former may be deprecated inthe future. + encouraged and the former may be deprecated in the future. 1. Inputs like `"y[1:3]"` are automatically expanded into `c("y[1]", "y[2]", "y[3]")`. In fact, simply `"y"` will be expanded into all nodes within `y`. @@ -279,17 +279,17 @@ except that they return a value (described above) and they have no `includeData` `simNodes`, `calcNodes` and `getLogProbNodes` are basic nimbleFunctions that simulate, calculate, or get the log probabilities (densities), respectively, of the same set of nodes each time they - are called. Each of these -takes a model and a character string of node names + are called. Each of these +takes a model and a character string of node names as inputs. If `nodes` is left blank, then all the nodes of the model - are used. - - For `simNodes`, the nodes provided will be topologically sorted to + are used. + + For `simNodes`, the nodes provided will be topologically sorted to simulate in the correct order. For `calcNodes` and `getLogProbNodes`, the nodes will be sorted and dependent nodes will be included. Recall that the calculations must be up to date (from a `calculate` call) for `getLogProbNodes` to return the values you are probably looking for. - + ```{r, Basic_Utils_Algs} simpleModelCode <- nimbleCode({ for(i in 1:4){ @@ -304,7 +304,7 @@ simpleModel <- nimbleModel(simpleModelCode, check = FALSE) cSimpleModel <- compileNimble(simpleModel) # simulates all the x's and y's -rSimXY <- simNodes(simpleModel, nodes = c('x', 'y') ) +rSimXY <- simNodes(simpleModel, nodes = c('x', 'y') ) # calls calculate on x and its dependents (y, but not z) rCalcXDep <- calcNodes(simpleModel, nodes = 'x') @@ -339,18 +339,18 @@ cGetLogProbXDep$run() cCalcXDep$run() ``` - + ### Accessing log probabilities via *logProb* variables {#sec:access-log-prob} -For each variable that contains at least one stochastic node, NIMBLE generates a model variable with the prefix 'logProb_'. In general users will not need to access these logProb variables directly but rather will use `getLogProb`. However, knowing they exist can be useful, in part because these variables can be monitored in an MCMC. +For each variable that contains at least one stochastic node, NIMBLE generates a model variable with the prefix 'logProb_'. In general users will not need to access these logProb variables directly but rather will use `getLogProb`. However, knowing they exist can be useful, in part because these variables can be monitored in an MCMC. When the stochastic node is scalar, the `logProb` variable will have the same size. For example: ```{r, usingModelLogProbs} -``` +``` Creation of `logProb` variables for stochastic multivariate nodes is trickier, because they can represent an arbitrary block of a larger @@ -363,11 +363,9 @@ example, in ```{r, eval = FALSE} for(i in 1:10) x[i,] ~ dmnorm(mu[], prec[,]) -``` +``` `x` may be $10 \times 20$ (dimensions must be provided), but `logProb_x` will be $10 \times 1$. For the most part you do not need to worry about how NIMBLE is storing the log probability values, because you can -always get them using `getLogProb`. - - +always get them using `getLogProb`. diff --git a/UserManual/src/chapter_WritingModels.Rmd b/UserManual/src/chapter_WritingModels.Rmd index bbd2e374d..6d3499dd8 100644 --- a/UserManual/src/chapter_WritingModels.Rmd +++ b/UserManual/src/chapter_WritingModels.Rmd @@ -4,7 +4,7 @@ ```{r, echo=FALSE} require(nimble) -``` +``` # Writing models in NIMBLE's dialect of BUGS {#cha-writing-models} @@ -56,8 +56,8 @@ NIMBLE extends the BUGS language in the following ways: 1. Link functions can be used in stochastic, as well as deterministic, declarations.^[But beware of the possibility of needing to set values for 'lifted' nodes created by NIMBLE.] 1. Data values can be reset, and which parts of a model are flagged as data can be changed, allowing one model to be used for different data sets without rebuilding the model each time. 1. One can use stochastic/dynamic indexes, i.e., indexes can be other nodes or functions of other nodes. For a given dimension of a node being indexed, if the index is not constant, it must be a scalar value. So expressions such as `mu[k[i], 3]` or `mu[k[i], 1:3]` or `mu[k[i], j[i]]` are allowed, but not `mu[k[i]:(k[i]+1)]`^[In some cases, one can write a nimbleFunction to achieve the desired result, such as replacing `sum(mu[start[i]:end[i]])` with a nimbleFunction that takes `mu`, `start[i]`, and `end[i]` as arguments and does the needed summation.] Nested dynamic indexes such as `mu[k[j[i]]]` are also allowed. - - + + ### Not-supported features of BUGS and JAGS {#sec:not-yet-supported} The following are not supported. @@ -98,19 +98,19 @@ The machine in this case is a graphical model^[Technically, a *directed acyclic graph*]. A *node* (sometimes called a *vertex*) holds one value, which may be a scalar or a vector. *Edges* define the relationships between nodes. A huge variety -of statistical models can be thought of as graphs. +of statistical models can be thought of as graphs. -Here is the code to define and create a simple linear regression model with four observations. +Here is the code to define and create a simple linear regression model with four observations. ```{r, chunk-WMinit, echo = FALSE} # source the code read_chunk(file.path('chunks', 'writingModels_chunks.R')) # one can put code chunks here if one wants -``` +``` ```{r, linearRegressionGraph, fig.height=10, fig.width=20, fig.cap="Graph of a linear regression model"} -``` +``` -The graph representing the model is +The graph representing the model is shown in Figure \@ref(fig:linearRegressionGraph). Each observation, `y[i]`, is a node whose edges say that it follows a normal distribution depending on a predicted value, `predicted.y[i]`, and @@ -122,7 +122,7 @@ which are each nodes. This graph is created from the following BUGS code: ```{r, linearRegressionCode, eval=FALSE} -``` +``` In this code, stochastic relationships are declared with '$\sim$' and deterministic relationships are declared with '`<-`'. For @@ -140,7 +140,7 @@ the everything gets drawn in the end. Available distributions, default and alte An equivalent graph can be created by this BUGS code: ```{r, linearRegressionAltCode, eval=FALSE} -``` +``` In this case, the `predicted.y[i]` nodes in Figure \@ref(fig:linearRegressionGraph) will be created automatically by @@ -157,30 +157,30 @@ NIMBLE's version of BUGS. ```{r, didacticnimbleCode, eval=FALSE} { # 1. normal distribution with BUGS parameter order - x ~ dnorm(a + b * c, tau) + x ~ dnorm(a + b * c, tau) # 2. normal distribution with a named parameter - y ~ dnorm(a + b * c, sd = sigma) + y ~ dnorm(a + b * c, sd = sigma) # 3. For-loop and nested indexing for(i in 1:N) { for(j in 1:M[i]) { - z[i,j] ~ dexp(r[ blockID[i] ]) + z[i,j] ~ dexp(r[ blockID[i] ]) } } # 4. multivariate distribution with arbitrary indexing - for(i in 1:3) + for(i in 1:3) mvx[8:10, i] ~ dmnorm(mvMean[3:5], cov = mvCov[1:3, 1:3, i]) # 5. User-provided distribution - w ~ dMyDistribution(hello = x, world = y) + w ~ dMyDistribution(hello = x, world = y) # 6. Simple deterministic node d1 <- a + b # 7. Vector deterministic node with matrix multiplication - d2[] <- A[ , ] %*% mvMean[1:5] + d2[] <- A[ , ] %*% mvMean[1:5] # 8. Deterministic node with user-provided function - d3 <- foo(x, hooray = y) + d3 <- foo(x, hooray = y) } -``` +``` -When a variable appears only on the right-hand side, it can be provided via `constants` (in which case it can never be changed) or via `data` or `inits`, as discussed in Chapter \@ref(cha-building-models). +When a variable appears only on the right-hand side, it can be provided via `constants` (in which case it can never be changed) or via `data` or `inits`, as discussed in Chapter \@ref(cha-building-models). Notes on the comment-numbered lines are: @@ -189,7 +189,7 @@ Notes on the comment-numbered lines are: 1. `y` follows a normal distribution with the same mean as `x` but a named standard deviation parameter instead of a precision parameter (sd = 1/sqrt(precision)). 1. `z[i, j]` follows an exponential distribution with parameter `r[ blockID[i] ]`. This shows how for-loops can be used for indexing of variables containing - multiple nodes. Variables that define for-loop indices (`N` and `M`) must also be provided as constants. + multiple nodes. Variables that define for-loop indices (`N` and `M`) must also be provided as constants. 1. The arbitrary block `mvx[8:10, i]` follows a multivariate normal distribution, with a named covariance matrix instead of BUGS' default of a precision matrix. As in R, curly braces for for-loop @@ -203,7 +203,7 @@ Notes on the comment-numbered lines are: function. See Chapter \@ref(cha-user-defined). -#### More about indexing {#sec:indexing} +#### More about indexing {#sec:indexing} Examples of allowed indexing include: @@ -216,8 +216,8 @@ Examples of allowed indexing include: - `x[k[i]+1]` \# a dynamic (and computed) index - `x[k[j[i]]]` \# nested dynamic indexes - `x[k[i], 1:3]` \# nested indexing of rows or columns - - + + NIMBLE does not allow multivariate nodes to be used without square brackets, which is an incompatibility with JAGS. Therefore a statement like `xbar <- mean(x)` in JAGS must be converted to `xbar <- mean(x[])` (if `x` is a vector) or `xbar <- @@ -225,7 +225,7 @@ square brackets, which is an incompatibility with JAGS. Therefore a statement l explained in later chapters, square brackets with blank indices are not necessary for multivariate objects.]. Section \@ref(sec:provide-dimensions) discusses how to provide NIMBLE with dimensions of `x` when needed. -One cannot provide a vector of indices that are not constant. For example, `x[start[i]:end[i]]` is allowed only if `start` and `end` are provided as constants. Also, one cannot provide a vector, or expression evaluating to a vector (apart from use of `:`) as an index. For example, `x[inds]` is not allowed. Often one can write a nimbleFunction to achieve the desired result, such as definining a nimbleFunction that takes `x`, `start[i]`, and `end[i]` or takes `x` and `inds` as arguments and does the subsetting (possibly in combination with some other calculation). Note that `x[L:U]` and `x[inds]` are allowed in nimbleFunction code. +One cannot provide a vector of indices that are not constant. For example, `x[start[i]:end[i]]` is allowed only if `start` and `end` are provided as constants. Also, one cannot provide a vector, or expression evaluating to a vector (apart from use of `:`) as an index. For example, `x[inds]` is not allowed. Often one can write a nimbleFunction to achieve the desired result, such as defining a nimbleFunction that takes `x`, `start[i]`, and `end[i]` or takes `x` and `inds` as arguments and does the subsetting (possibly in combination with some other calculation). Note that `x[L:U]` and `x[inds]` are allowed in nimbleFunction code. Generally NIMBLE supports R-like linear algebra expressions and attempts to follow the same rules as R about dimensions (although in some cases this is not possible). For @@ -249,7 +249,7 @@ this would be created with a for loop: logY[i] <- log(Y[i]) } } -``` +``` Since NIMBLE supports R-like algebraic expressions, an alternative in NIMBLE's dialect of BUGS is to use a vectorized declaration like this: @@ -257,7 +257,7 @@ NIMBLE's dialect of BUGS is to use a vectorized declaration like this: { logY[1:10] <- log(Y[1:10]) } -``` +``` There is an important difference between the models that are created by the @@ -312,7 +312,7 @@ Here is some example code illustrating how to use the distribution in a model. code <- nimbleCode({ Ustar[1:p,1:p] ~ dlkj_corr_cholesky(1.3, p) U[1:p,1:p] <- uppertri_mult_diag(Ustar[1:p, 1:p], sds[1:p]) - for(i in 1:n) + for(i in 1:n) y[i, 1:p] ~ dmnorm(mu[1:p], cholesky = U[1:p, 1:p], prec_param = 0) }) ``` @@ -366,9 +366,9 @@ as an example: 1. A *canonical* parameterization is used directly for computations^[Usually this is the parameterization in the `Rmath` header of R's C implementation of distributions.]. For - `gamma`, this is (shape, scale). + `gamma`, this is (shape, scale). 1. The BUGS parameterization is the one defined in the - original BUGS language. In general this is the parameterization for which conjugate MCMC samplers can be executed most efficiently. For `dgamma`, this is (shape, rate). + original BUGS language. In general this is the parameterization for which conjugate MCMC samplers can be executed most efficiently. For `dgamma`, this is (shape, rate). 1. An *alternative* parameterization is one that must be converted into the *canonical* parameterization. For `dgamma`, NIMBLE provides both (shape, rate) and (mean, sd) parameterization @@ -385,7 +385,7 @@ parameter names are not given, they are taken in order, so that (shape, rate) is the default. This happens to match R's order of parameters, but it need not. If names are given, they can be given in any order. NIMBLE knows that rate is an alternative to scale and that -(mean, sd) are an alternative to (shape, scale or rate). +(mean, sd) are an alternative to (shape, scale or rate). ```{r, child = 'tables/parameterizationTableLong.md'} ``` @@ -401,23 +401,23 @@ would result in very inefficient AD calculations. Note that in this case, NIMBLE automatically determines the inverse of the provided matrix (which could be either the precision or the covariance), so any model calculations that need that inverse can simply use `getParam` at no additional -computational expense. +computational expense. - -In addition, NIMBLE supports alternative distribution names, known as aliases, as in JAGS, as specified in Table \@ref(tab:densityAliases). +In addition, NIMBLE supports alternative distribution names, known as aliases, as in JAGS, as specified in Table \@ref(tab:densityAliases). ```{r, child = 'tables/densityAliasesTable.md'} ``` -We plan to, but do not currently, include the following distributions as part of core NIMBLE: beta-binomial, Dirichlet-multinomial, F, Pareto, or forms of the multivariate t other than the standard one provided. +We plan to, but do not currently, include the following distributions as part of core NIMBLE: beta-binomial, Dirichlet-multinomial, F, Pareto, or forms of the multivariate t other than the standard one provided. @@ -425,30 +425,30 @@ We plan to, but do not currently, include the following distributions as part of ### Available BUGS language functions {#subsec:BUGS-lang-fxns} Tables \@ref(tab:functions)-\@ref(tab:functionsMatrix) show the available operators and functions. - +--> Support for more general R expressions is covered in Chapter \@ref(cha-RCfunctions) about programming -with nimbleFunctions. +with nimbleFunctions. For the most part NIMBLE supports the functions used in BUGS and JAGS, with exceptions indicated in the table. Additional functions provided by NIMBLE are also listed. Note that we provide distribution functions for use in calculations, namely the 'p', 'q', and 'd' functions. - See Section \@ref(sec:nimble-dist-funs) for details on the syntax for using distribution functions as functions in deterministic calculations, as only some parameterizations are allowed and the names of some distributions differ from those used to define stochastic nodes in a model. - - - - - - ```{r, child = 'tables/functionTableLong.md'} @@ -461,7 +461,7 @@ system("cat tables/functionTableMatrixLong.md tables/eigenSvdTableMatrixElements ```{r echo=FALSE, include=FALSE} system("rm -rf tmp.md") ``` - + See Section \@ref(sec:user-functions) to learn how to use nimbleFunctions to write new functions for use in BUGS code. @@ -472,25 +472,25 @@ NIMBLE allows the link functions listed in Table \@ref(tab:links). ```{r, child = 'tables/linksTable.md'} ``` - + Link functions are specified as functions applied to a node on the left hand side of a BUGS expression. To handle link functions in deterministic declarations, NIMBLE converts the declaration into an equivalent inverse declaration. For example, `log(y) <- x` is converted into `y <- exp(x)`. In other words, the link function is -just a simple variant for conceptual clarity. +just a simple variant for conceptual clarity. To handle link functions in a stochastic declaration, NIMBLE does some processing that inserts an additional node into the model. For example, the declaration `logit(p[i]) ~ dnorm(mu[i],1)`, is equivalent -to the following two declarations: +to the following two declarations: - `logit_p[i] ~ dnorm(mu[i], 1)`, - `p[i] <- expit(logit_p[i])` -where `expit` is the inverse of `logit`. +where `expit` is the inverse of `logit`. -Note that NIMBLE does not provide an automatic way of initializing the additional node (`logit_p[i]` in this case), which is a parent node of the explicit node (`p[i]`), without explicitly referring to the additional node by the name that NIMBLE generates. +Note that NIMBLE does not provide an automatic way of initializing the additional node (`logit_p[i]` in this case), which is a parent node of the explicit node (`p[i]`), without explicitly referring to the additional node by the name that NIMBLE generates. ### Truncation, censoring, and constraints {#subsec:trunc} @@ -498,15 +498,15 @@ NIMBLE provides three ways to declare boundaries on the value of a variable, eac #### Truncation -Either of the following forms, +Either of the following forms, ``` x ~ dnorm(0, sd = 10) T(0, a) x ~ T(dnorm(0, sd = 10), 0, a) ``` - + declares that `x` follows a normal distribution between 0 and - `a` (inclusive of 0 and `a`). Either boundary may be omitted or may be another node, such as `a` in this example. The first form is compatible with JAGS, but in NIMBLE it can only be used when reading code from a text file. When writing model code in R, the second version must be used. + `a` (inclusive of 0 and `a`). Either boundary may be omitted or may be another node, such as `a` in this example. The first form is compatible with JAGS, but in NIMBLE it can only be used when reading code from a text file. When writing model code in R, the second version must be used. Truncation means the possible values of `x` are limited a priori, hence the probability density of `x` must be normalized^[NIMBLE uses the CDF and inverse CDF (quantile) functions of a distribution to do this; in some cases if one uses truncation to include only the extreme tail of a distribution, numerical difficulties can arise.]. In this example it would be the normal probability density divided by its integral from 0 to `a`. Like JAGS, NIMBLE also provides `I` as a synonym for `T` to accommodate older BUGS code, but `T` is preferred because it disambiguates multiple usages of `I` in BUGS. @@ -526,11 +526,11 @@ nodes will be sampled in an MCMC.) The vector for `c` should give the censoring times corresponding to censored entries and a value above the observed times for uncensored entries (e.g., `Inf`). The values for `c` can be provided via `constants`, `inits`, or `data`. We recommend providing -via `constants` for values that are fixed and via `inits` for values that are +via `constants` for values that are fixed and via `inits` for values that are to be estimated or might be changed for some reason. -Left-censored observations would be specified by setting `censored[i]` to 0 -and `t[i]` to `NA`. +Left-censored observations would be specified by setting `censored[i]` to 0 +and `t[i]` to `NA`. Important: The value of `t[i]` corresponding to each censored data point should be given a valid initial value (via `inits`) that is consistent with `censored[i]`. Otherwise, NIMBLE will initialize @@ -544,25 +544,25 @@ The `dinterval` is not really a distribution but rather a trick: in the above ex `censored[i] = M` if `c[i, M] < t[i]` -(The `i` index is provided only for consistency with the previous example.) The most common uses of `dinterval` will be for left- and right-censored data, in which case `c[i,]` will be a single value (and typically given as simply `c[i]`), and for interval-censored data, in which case `c[i,]` will be a vector of two values. - Nodes following a `dinterval` distribution should normally be set -as `data` with known values. Otherwise, the node may be simulated during initialization in some algorithms (e.g., MCMC) and thereby establish a permanent, perhaps unintended, constraint. +as `data` with known values. Otherwise, the node may be simulated during initialization in some algorithms (e.g., MCMC) and thereby establish a permanent, perhaps unintended, constraint. -Censoring differs from truncation because censoring an observation involves bounds on a random variable that could have taken any value, while in truncation we know a priori that a datum could not have occurred outside the truncation range. +Censoring differs from truncation because censoring an observation involves bounds on a random variable that could have taken any value, while in truncation we know a priori that a datum could not have occurred outside the truncation range. #### Constraints and ordering NIMBLE provides a more general way to enforce constraints using `dconstraint(cond)`. For example, we could specify that the sum of `mu1` and `mu2` must be positive like this: ```{r, dconstraint-example, eval=FALSE} -mu1 ~ dnorm(0, 1) -mu2 ~ dnorm(0, 1) +mu1 ~ dnorm(0, 1) +mu2 ~ dnorm(0, 1) constraint_data ~ dconstraint( mu1 + mu2 > 0 ) -``` +``` with `constraint_data` set (as `data`) to 1. This is equivalent to a half-normal distribution on the half-plane $\mu_1 + \mu_2 > 0$. Nodes following `dconstraint` should be provided as data for the same reason of avoiding unintended initialization described above for `dinterval`. @@ -570,8 +570,8 @@ equivalent to a half-normal distribution on the half-plane $\mu_1 + Important: the model should be initialized so that the constraints are satisfied. Otherwise, NIMBLE will initialize from the prior, which may not produce a valid initial value and may cause algorithms (particularly) MCMC to fail. -Formally, `dconstraint(cond)` is a probability distribution on $\left\{ 0, 1 \right\}$ such that $P(1) = 1$ if `cond` is `TRUE` and $P(0) = 1$ if `cond` is `FALSE`. - Of course, in many cases, parameterizing the model so that the @@ -585,10 +585,10 @@ algorithms that are not tailored to the constraints may fail or be inefficient. ##### Ordering -To specify an ordering of parameters, such as $\alpha_1 <= \alpha_2 <= \alpha_3$ one can use `dconstraint` as follows: +To specify an ordering of parameters, such as $\alpha_1 <= \alpha_2 <= \alpha_3$ one can use `dconstraint` as follows: ```{r, ordering-example, eval=FALSE} constraint_data ~ dconstraint( alpha1 <= alpha2 & alpha2 <= alpha3 ) -``` +``` Note that unlike in BUGS, one cannot specify prior ordering using syntax such as @@ -598,8 +598,8 @@ alpha[2] ~ dnorm(0, 1) I(alpha[1], alpha[3]) alpha[3] ~ dnorm(0, 1) I(alpha[2], ) ``` -as this does not represent a directed acyclic graph. - Also note that specifying prior ordering using `T(,)` can result in possibly unexpected results. For example: @@ -677,4 +677,4 @@ mod <- nimbleModel(code, constants = chickConst, data = chickData) mod$getCode() ``` -As discussed in Section \@ref(sec:user-macros), NIMBLE users can also create new macros to use in models. \ No newline at end of file +As discussed in Section \@ref(sec:user-macros), NIMBLE users can also create new macros to use in models. diff --git a/UserManual/src/chapter_WritingNimbleFunctions.Rmd b/UserManual/src/chapter_WritingNimbleFunctions.Rmd index 282aa87b9..fc52a400a 100644 --- a/UserManual/src/chapter_WritingNimbleFunctions.Rmd +++ b/UserManual/src/chapter_WritingNimbleFunctions.Rmd @@ -4,16 +4,16 @@ ```{r, echo=FALSE} require(nimble) -``` +``` # Writing nimbleFunctions to interact with models {#cha-progr-with-models} ```{r, echo=FALSE} -read_chunk(file.path('chunks', 'programmingWithModels_chunks.R')) -``` +read_chunk(file.path('chunks', 'programmingWithModels_chunks.R')) +``` + - ## Overview {#sec:writ-nimble-funct} When you write an R function, you say what the input arguments are, @@ -42,18 +42,18 @@ are two kinds of code and two steps of execution: arguments, the return value, and in some circumstances for local variables. There are two kinds of *run* code: - + a. There is always a primary function, given as the argument `run`^[This can be omitted if you don't need it.]. b. There can optionally be other functions, or 'methods' in the language of object-oriented programming, that share the same objects created by the *setup* function. - + Here is a small example to fix ideas: ```{r, nf-intro} -``` +``` The call to the R function called `nimbleFunction` returns a function, similarly to defining a function in R. That function, @@ -84,33 +84,33 @@ structure is determined, and then the `run` code can use that information without repeatedly, redundantly recomputing it. A nimbleFunction can be called repeatedly (one can think of it as a generator), each time returning a specialized -nimbleFunction. - +nimbleFunction. + Readers familiar with object-oriented programming may find it useful to think in terms of class definitions and objects. `nimbleFunction` creates a class definition. Each specialized nimbleFunction is one object in the class. The setup arguments are used to define member data in -the object. +the object. ## Using and compiling nimbleFunctions {#sec:using-comp-nimbl} To compile the nimbleFunction, together with its model, we use `compileNimble`: ```{r, nf-compiling} -``` +``` These have been initialized with the values from their uncompiled versions and can be used in the same way: ```{r, nf-using} -``` +``` ## Writing setup code ### Useful tools for setup functions The setup function is typically used to determine information on nodes -in a model, set up modelValues or nimbleList objects, set up (nested) nimbleFunctions +in a model, set up modelValues or nimbleList objects, set up (nested) nimbleFunctions or nimbleFunctionLists, and set up any persistent numeric objects. For example, the setup code of an MCMC nimbleFunction creates the nimbleFunctionList of sampler nimbleFunctions. The values of numeric @@ -125,9 +125,9 @@ Some of the useful tools and objects to create in setup functions include: - **modelValues objects** These are discussed in Sections \@ref(sec:modelValues-struct) and \@ref(sec:access-model-modelv). - **nimbleList objects** New instances of `nimbleList`s can be created from a nimbleList definition in either setup or run code. See Section \@ref(sec:nimbleLists) for more information. - - **specializations of other nimbleFunctions** A useful NIMBLE programming technique is to have one nimbleFunction contain other @@ -137,7 +137,7 @@ Some of the useful tools and objects to create in setup functions include: nimbleFunctions (Section \@ref(sec:virt-nimbl-nimbl)). If one wants a nimbleFunction that does get specialized but has -empty setup code, use `setup = function() {}` or `setup = TRUE`. +empty setup code, use `setup = function() {}` or `setup = TRUE`. ### Accessing and modifying numeric values from setup {#sec:access-modify-numer} @@ -146,10 +146,10 @@ modified^[Actually, they can be, but only for uncompiled nimbleFunctions.], numeric values and modelValues can be, as illustrated by extending the example from above. ```{r, nf-modifyValueToAdd} -``` +``` ### Determining numeric types in nimbleFunctions - + For numeric variables from the `setup` function that appear in the `run` function or other member functions (or are declared in `setupOutputs`), the @@ -159,7 +159,7 @@ For numeric variables from the `setup` function that example if `X` is created as a matrix (two-dimensional double) in one specialization but as a vector (one-dimensional double) in another, there will be a problem during compilation. The sizes may differ in each specialization. - + Treatment of vectors of length one presents special challenges because they could be treated as scalars or vectors. Currently they are treated as scalars. If you want a vector, ensure that the length is @@ -181,13 +181,13 @@ example `X` and `Y`, simply include ```{r, eval = FALSE} setupOutputs(X, Y) -``` +``` -anywhere in the setup code. +anywhere in the setup code. ## Writing run code {#sec:nimble-lang-comp} -In Chapter \@ref(cha-RCfunctions) we described the functionality of the NIMBLE language that could be used in run code without setup code (typically in cases where no models or modelValues are needed). Next we explain the additional features that allow use of models and modelValues in the run code. +In Chapter \@ref(cha-RCfunctions) we described the functionality of the NIMBLE language that could be used in run code without setup code (typically in cases where no models or modelValues are needed). Next we explain the additional features that allow use of models and modelValues in the run code. ### Driving models: *calculate*, *calculateDiff*, *simulate*, *getLogProb* {#sec:driv-models:-calc} @@ -197,7 +197,7 @@ for `getLogProb`, it is usually important for the `nodes` vector to be sorted in topological order. Model member functions such as `getDependencies` and `expandNodeNames` will -always return topoligically sorted node names. +always return topologically sorted node names. Most R-like indexing of a node vector is allowed within the argument to `calculate`, `calculateDiff`, `simulate`, and `getLogProb`. For example, all of the following are allowed: @@ -209,7 +209,7 @@ myModel$calculate(nodes[1:3]) myModel$calculate(nodes[c(1,3)]) myModel$calculate(nodes[2:i]) myModel$calculate(nodes[ values(model, nodes) + 0.1 < x ]) -``` +``` Note that in the last line of code, one must have that the length of `nodes` is equal to that of `values(model, nodes)`, which means that all the nodes in `nodes` @@ -222,7 +222,7 @@ can only be indexed within a call to `calculate`, `calculateDiff`, ### Getting and setting variable and node values #### Using indexing with nodes - + Here is an example that illustrates getting and setting of nodes, subsets of nodes, or variables. Note the following: @@ -231,7 +231,7 @@ subsets of nodes, or variables. Note the following: - In `model[[v]]`, `v` can only be a single node or variable name, not a vector of multiple nodes nor an element of such a vector (`model[[ nodes[i] ]]` does not work). The node itself may be a vector, matrix or array node. - In fact, `v` can be a node-name-like character string, even if it is not actually a node in the model. See example 4 in the code below. - One can also use `model$varName`, with the caveat that `varName` must be a variable name. This usage would only make sense for a nimbleFunction written for models known to have a specific variable. (Note that if `a` is a scalar node in `model`, then `model[['a']]` will be a scalar but `model$a` will be a vector of length 1. - - one should use the `<<-` global assignment operator to assign into model nodes. + - one should use the `<<-` global assignment operator to assign into model nodes. Note that NIMBLE does not allow variables to change dimensions. Model nodes are the same, and indeed are more restricted because they can't change sizes. In addition, NIMBLE distinguishes between scalars and vectors of length 1. These rules, and ways to handle them correctly, are illustrated in the following code as well as in Section \@ref(sec:how-nimble-handles). @@ -253,7 +253,7 @@ nfGen <- nimbleFunction( # node1 and node2 would typically be setup arguments, so they could # have different values for different models. We are assigning values # here so the example is clearer. - node1 <- 'sigma' # a scalar node + node1 <- 'sigma' # a scalar node node2 <- 'y[1:5]' # a vector node notReallyANode <- 'y[2:4]' # y[2:4] allowed even though not a node! }, @@ -264,19 +264,19 @@ nfGen <- nimbleFunction( tmp3 <- model[[notReallyANode]] # 4. tmp3 will be a vector tmp4 <- model$y[3:4] # 5. hard-coded access to a model variable # 6. node1 is scalar so can be assigned a scalar: - model[[node1]] <<- runif(1) - model[[node2]][1] <<- runif(1) - # 7. an element of node2 can be assigned a scalar + model[[node1]] <<- runif(1) + model[[node2]][1] <<- runif(1) + # 7. an element of node2 can be assigned a scalar model[[node2]] <<- runif(length(model[[node2]])) # 8. a vector can be assigned to the vector node2 - model[[node2]][1:3] <<- vals[1:3] + model[[node2]][1:3] <<- vals[1:3] # elements of node2 can be indexed as needed returnType(double(1)) - out <- model[[node2]] # we can return a vector + out <- model[[node2]] # we can return a vector return(out) } ) - + Rnf <- nfGen(m) Cnf <- compileNimble(Rnf) Cnf$run(rnorm(10)) @@ -300,7 +300,7 @@ probabilities. This can be done using `values`: P <- values(model, nodes) # or put values from a vector into a set of model nodes values(model, nodes) <- P -``` +``` where the first line would assign the collection of values from `nodes` into `P`, and the second would do the inverse. In both cases, values @@ -311,7 +311,7 @@ column-wise order. vector in other expressions, e.g., ```{r, eval=FALSE} Y <- A %*% values(model, nodes) + b -``` +``` One can also index elements of nodes in the argument to values, in the same manner as discussed for `calculate` and related functions in Section \@ref(sec:driv-models:-calc). @@ -324,18 +324,18 @@ For example: ```{r, eval=FALSE} # c(rnorm(1)) creates vector of length one: -values(model, nodes[1]) <- c(rnorm(1)) +values(model, nodes[1]) <- c(rnorm(1)) # won't compile because rnorm(1) is a scalar -# values(model, nodes[1]) <- rnorm(1) +# values(model, nodes[1]) <- rnorm(1) out <- values(model, nodes[1]) # out is a vector out2 <- values(model, nodes[1])[1] # out2 is a scalar -``` +``` ### Getting parameter values and node bounds -Sections \@ref(sec:getParam)-\@ref(sec:getBound) describe how to get the parameter values for a node and the range (bounds) of possible values for the node using `getParam` and `getBound`. Both of these can be used in run code. +Sections \@ref(sec:getParam)-\@ref(sec:getBound) describe how to get the parameter values for a node and the range (bounds) of possible values for the node using `getParam` and `getBound`. Both of these can be used in run code. ### Using modelValues objects {#sec:access-model-modelv} @@ -346,7 +346,7 @@ designed to easily save values from a model object during the running of a nimbleFunction. A modelValues object used in run code must always exist in the setup code, either by passing it in as a setup argument or creating it in the setup code. - + To illustrate this, we will create a nimbleFunction for computing importance weights for importance sampling. This function will use two modelValues objects. `propModelValues` will contain a set of @@ -354,35 +354,35 @@ values simulated from the importance sampling distribution and a field `propLL` for their log probabilities (densities). `savedWeights` will contain the difference in log probability (density) between the model and the -`propLL` value provided for each set of values. +`propLL` value provided for each set of values. ```{r, mv-setup-code} ``` - + The simplest way to pass values back and forth between models and modelValues inside of a nimbleFunction is with `copy`, which has the synonym `nimCopy`. See `help(nimCopy)` for argument details. - - + Alternatively, the values may be accessed via indexing of individual rows, using the notation `mv[var, i]`, where `mv` is a modelValues object, `var` is a variable name (not a node name), and `i` is a row number. Likewise, the `getsize` and `resize` functions can be used as discussed in Section \@ref(sec:modelValues-struct). However the function `as.matrix` does not work in run code. - + Here is a run function to use these modelValues: - + ```{r, mv-run-time} ``` @@ -391,7 +391,7 @@ using `$`, which is shown in more detail below. In fact, since modelValues, like most NIMBLE objects, are reference class objects, one can get a reference to them before the function is executed and then use that reference afterwards. - + ```{r, mv-compilation-example} ``` @@ -408,12 +408,12 @@ statements would be valid: ```{r, eval = FALSE} model[["x[2:8, ]"]][2:4, 1:3] %*% Z -``` +``` if Z is a vector or matrix, and ```{r, eval = FALSE} C[6:10] <- mv[v, i][1:5, k] + B -``` +``` if B is a vector or matrix. The NIMBLE language allows scalars, but models defined from BUGS code @@ -432,7 +432,7 @@ setup code in just the same ways as the run function. In fact, the run function is just a default main method name. Any method can then call another method. ```{r, usingMemberFunctions} -``` +``` ### Using other nimbleFunctions {#sec:using-other-nimbl} @@ -440,10 +440,10 @@ One nimbleFunction can use another nimbleFunction that was passed to it as a setup argument or was created in the setup function. This can be an effective way to program. When a nimbleFunction needs to access a setup variable or method of another nimbleFunction, use -`$`. +`$`. ```{r, owningMemberFunctions} -``` +``` ### Virtual nimbleFunctions and nimbleFunctionLists {#sec:virt-nimbl-nimbl} @@ -464,7 +464,7 @@ simple, single-level inheritance. Here is how it works: ```{r, nimbleFunctionLists} -``` +``` One can also use `seq_along` with nimbleFunctionLists (and only with nimbleFunctionLists). As in R, `seq_along(myFunList)` is equivalent to @@ -479,7 +479,7 @@ Virtual nimbleFunctions cannot define setup values to be inherited. NIMBLE provides limited uses of character objects in run code. Character vectors created in setup code will be available in run code, but the only thing you can really do with them is -include them in a `print` or `stop` statement. +include them in a `print` or `stop` statement. Note that character vectors of model node and variable names are processed during compilation. For example, in `model[[node]]`, `node` @@ -498,7 +498,7 @@ Before the introduction of nimbleLists in Version 0.6-4, NIMBLE did not explicit they are compiled using `setupOutputs`. Here is an example: ```{r, dataStructures} -``` +``` You'll notice that: @@ -516,7 +516,7 @@ interfaces by setting the `buildInterfacesForCompiledNestedNimbleFunctions` option via `nimbleOptions` to TRUE. If we had left that option FALSE (its default value), we could still get to the values of interest using -`valueInCompiledNimbleFunction(CmyDataNF, 'X')` +`valueInCompiledNimbleFunction(CmyDataNF, 'X')` - We need to take care that at the time of compilation, the `X`, `Y` and `Z` values contain doubles via `as.numeric` so that they are not compiled as integer objects. @@ -526,34 +526,34 @@ default value), we could still get to the values of interest using ## Example: writing user-defined samplers to extend NIMBLE's MCMC engine {#sec:user-samplers} -One important use of nimbleFunctions is to write additional samplers that can be used in NIMBLE's MCMC engine. This allows a user to write a custom sampler for one or more nodes in a model, as well as for programmers to provide general samplers for use in addition to the library of samplers provided with NIMBLE. +One important use of nimbleFunctions is to write additional samplers that can be used in NIMBLE's MCMC engine. This allows a user to write a custom sampler for one or more nodes in a model, as well as for programmers to provide general samplers for use in addition to the library of samplers provided with NIMBLE. The following code illustrates how a NIMBLE developer would implement and use a Metropolis-Hastings random walk sampler with fixed proposal standard deviation. ```{r, custom-sampler, eval=FALSE} my_RW <- nimbleFunction( - + contains = sampler_BASE, - + setup = function(model, mvSaved, target, control) { # proposal standard deviation scale <- if(!is.null(control$scale)) control$scale else 1 calcNodes <- model$getDependencies(target) }, - + run = function() { # initial model logProb - model_lp_initial <- getLogProb(model, calcNodes) + model_lp_initial <- getLogProb(model, calcNodes) # generate proposal - proposal <- rnorm(1, model[[target]], scale) + proposal <- rnorm(1, model[[target]], scale) # store proposal into model - model[[target]] <<- proposal + model[[target]] <<- proposal # proposal model logProb model_lp_proposed <- model$calculate(calcNodes) - + # log-Metropolis-Hastings ratio log_MH_ratio <- model_lp_proposed - model_lp_initial - + # Metropolis-Hastings step: determine whether or # not to accept the newly proposed value u <- runif(1, 0, 1) @@ -561,24 +561,24 @@ my_RW <- nimbleFunction( else jump <- FALSE # keep the model and mvSaved objects consistent - if(jump) copy(from = model, to = mvSaved, row = 1, + if(jump) copy(from = model, to = mvSaved, row = 1, nodes = calcNodes, logProb = TRUE) else copy(from = mvSaved, to = model, row = 1, nodes = calcNodes, logProb = TRUE) }, - + methods = list( reset = function () {} ) ) -``` +``` -The name of this sampler function, for the purposes of using it in an -MCMC algorithm, is `my_RW`. Thus, this sampler can be added -to an exisiting MCMC configuration object `conf` using: +The name of this sampler function, for the purposes of using it in an +MCMC algorithm, is `my_RW`. Thus, this sampler can be added +to an existing MCMC configuration object `conf` using: ```{r, custom-sampler-add, eval=FALSE} mcmcConf$addSampler(target = 'x', type = 'my_RW', control = list(scale = 0.1)) -``` +``` To be used within the MCMC engine, sampler functions definitions must adhere exactly to the following: @@ -590,7 +590,7 @@ adhere exactly to the following: no return value. Further, after execution it must leave the `mvSaved` modelValues object as an up-to-date copy of the values and logProb values in the model object. - - The nimbleFunction must have a member method called `reset`, which takes no arguments + - The nimbleFunction must have a member method called `reset`, which takes no arguments and has no return value. @@ -621,7 +621,7 @@ As of version 0.13.0, NIMBLE's handling of posterior predictive nodes in MCMC sa - turning off posterior predictive nodes as dependencies in the use of `getDependencies` when building an MCMC, and - moving all posterior predictive samplers to be last in the order of samplers used in an MCMC iteration. - + The first change calls for careful consideration when writing new samplers. It is done by making `buildMCMC` set a NIMBLE system option (described below) that tells `getDependencies` to ignore posterior predictive nodes (defined as nodes that are themselves not data and have no data nodes in their entire downstream (descendant) dependency network) before it builds the samplers (i.e., runs the setup code of each sampler). That way, the setup code of all samplers that use `getDependencies` will automatically comply with the new system. User-defined samplers that determine dependencies using `getDependencies` should automatically work in the new system (although it is worth checking). However if a user-defined sampler manually specifies dependencies, and if those dependencies include posterior predictive nodes, the MCMC results could be incorrect. Whether incorrect results occur could depend on other parts of the model structure connected to node(s) sampled by the user-defined sampler (i.e., the target node(s)) and/or posterior predictive nodes. Therefore, it is strongly recommended that all user-defined samplers rely on `getDependencies` to determine which nodes depend on the target node(s). @@ -629,7 +629,7 @@ User-defined samplers that determine dependencies using `getDependencies` should There are several new NIMBLE system options available to take control of the behavior just described. - `getDependenciesIncludesPredictiveNodes` determines whether posterior predictive nodes are included in results of `getDependencies`. This defaults to `TRUE`, so that outside of building samplers, the behavior of `getDependencies` should be identical to previous versions of NIMBLE. - - `MCMCusePredictiveDependenciesInCalculations` gives the value to which + - `MCMCusePredictiveDependenciesInCalculations` gives the value to which `getDependenciesIncludesPredictiveNodes` will be set while building samplers. This defaults to `FALSE`. - `MCMCorderPosteriorPredictiveSamplersLast` determines whether posterior predictive samplers are moved to the end of the sampler list. This default to `TRUE`. Behavior prior to version 0.13.0 corresponds to `FALSE`. @@ -638,7 +638,7 @@ Here are some examples of how to use these options: - If you want to have MCMC samplers condition on posterior predictive nodes, do `nimbleOptions(MCMCusePredictiveDependenciesInCalculations = TRUE)`. - If you want to get identical MCMC samplers, including their ordering, as in versions prior to 0.13.0, do `nimbleOptions(MCMCusePredictiveDependenciesInCalculations = TRUE)` and `nimbleOptions(MCMCorderPosteriorPredictiveSamplersLast = FALSE)`. - If you want to experiment with the behavior of `getDependencies` when ignoring posterior predictive nodes, do `nimbleOptions(getDependenciesIncludesPredictiveNodes = FALSE)`. - + ## Copying nimbleFunctions (and NIMBLE models) NIMBLE relies heavily on R's reference class system. When models, @@ -657,7 +657,7 @@ the model member function `newModel` will create a new copy of the model from the same model definition (Section \@ref(sub:multiple-instances)). This new model can then be used with newly instantiated nimbleFunctions. - @@ -672,7 +672,7 @@ tell you during compilation that the second set of `foo` instances cannot be built from the previous compiled version. A solution is to re-define `foo` from the beginning -- i.e. call `nimbleFunction` again -- and then proceed with building and compiling the algorithm -for `model2`. +for `model2`. ## Debugging nimbleFunctions {#sec:debugging} @@ -681,7 +681,7 @@ of each nimbleFunction is for debugging. One can call `debug` on nimbleFunction methods (in particular the main *run* method, e.g., `debug(mynf$run`) and then step through the code in R using R's debugger. One can also insert `browser` calls into run code and then run the -nimbleFunction from R. +nimbleFunction from R. In contrast, directly debugging a compiled nimbleFunction is difficult, although those familiar with running R through a debugger @@ -747,7 +747,7 @@ controlled via `nimbleOptions`. As noted above, the option `buildInterfacesForCompiledNestedNimbleFunctions` defaults to FALSE, which means NIMBLE will not build full interfaces to compiled -nimbleFunctions that ony appear within other nimbleFunctions. If you +nimbleFunctions that only appear within other nimbleFunctions. If you want access to all such nimbleFunctions, use the option `buildInterfacesForCompiledNestedNimbleFunctions = TRUE`. This will use more memory but can be useful for debugging. @@ -758,6 +758,5 @@ contents of an uncompiled nimbleFunction object after it has been compiled in an effort to free some memory. We expect to be able to keep making NIMBLE more efficient -- faster execution and lower memory use -- in the future. - -