Sequential estimation of Gaussian integrals

gaussint is used for calculating $n$-dimensional Gaussian integrals $$\int_a^b \frac{|Q|^{1/2}}{(2\pi)^{n/2}} \exp(-\frac1{2}(x-\mu)^{T}Q(x-\mu)) dx$$ A limit value $lim$ can be used to stop the integration if the sequential estimate goes below the limit, which can result in substantial computational savings in cases when one only is interested in testing if the integral is above the limit value. The integral is calculated sequentially, and estimates for all subintegrals are also returned.

Usage

gaussint(
  mu,
  Q.chol,
  Q,
  a,
  b,
  lim = 0,
  n.iter = 10000,
  ind,
  use.reordering = c("natural", "sparsity", "limits"),
  max.size,
  max.threads = 0,
  seed,
  tol = NULL,
  tol.level = NULL,
  size.tol = NULL
)

Arguments

mu

Expectation vector for the Gaussian distribution.

Q.chol

The Cholesky factor of the precision matrix (optional).

Q

Precision matrix for the Gaussian distribution. If Q is supplied but not Q.chol, the cholesky factor is computed before integrating.

a

Lower limit in integral.

b

Upper limit in integral.

lim

If this argument is used, the integration is stopped and 0 is returned if the estimated value goes below $lim$.

n.iter

Number or iterations in the MC sampler that is used for approximating probabilities. The default value is 10000.

ind

Indices of the nodes that should be analyzed (optional).

use.reordering

Determines what reordering to use:

"natural": No reordering is performed.
"sparsity": Reorder for sparsity in the cholesky factor (MMD reordering is used).
"limits": Reorder by moving all nodes with a=-Inf and b=Inf first and then reordering for sparsity (CAMD reordering is used).

max.size

The largest number of sub-integrals to compute. Default is the total dimension of the distribution.

max.threads

The number of threads that the program can use. The default, 0, uses the default number of threads of OpenMP, which can be set with the environment variable OMP_NUM_THREADS. The number of threads is at most OMP_THREAD_LIMIT, and is one if the package was built without OpenMP.

seed

The random seed to use (optional).

tol

Target for the estimated error (optional). If tol is given, the number of iterations is chosen adaptively: the integral is first estimated with 1000 iterations, and iterations are then added until the estimated error is at most tol, using at most n.iter iterations in total. By default, n.iter iterations are always used.

tol.level

The probability level where the error is controlled if tol is given (optional). The error is then controlled for the sub-integral estimates where the estimates go below tol.level, which is useful when one is interested in where the sub-integrals pass a probability level. By default, the error of P is controlled.

size.tol

Target for the estimated error of the number of sub-integrals before the estimates go below tol.level, relative to the number (optional). If size.tol is given, and tol is not, the number of iterations is chosen adaptively as for tol, but with this target, see the details. It requires tol.level.

Value

A list with elements

P: Value of the integral.
E: Estimated error of the P estimate.
Pv: A vector with the estimates of all sub-integrals.
Ev: A vector with the estimated errors of the Pv estimates.
n.iter: The number of iterations that were used.

Details

The function uses sequential importance sampling to estimate the Gaussian integral, and returns all computed sub-integrals. This means that if, for example, the function is used to compute $P(x>0)$ for an n-dimensional Gaussian variable $x$, then all integrals $P(x_1>0,\ldots,x_i>0)$ for $i=1,\ldots,n$ are computed.

If one is only interested in whether $P(x>0)>\alpha$ or not, then one can stop the integration as soon as $P(x_1>0,\ldots,x_i>0)<\alpha$. This can save a lot of computation time if $P(x_1>0,\ldots,x_i>0)< \alpha$ for $i$ much smaller than $n$. This limit value is specified by the lim argument.

Which reordering to use depends on what the purpose of the calculation is and what the integration limits are. However, in general the limits reordering is typically most appropriate since this combines sparisty (which improves accuracy and reduces computational cost) with automatic handling of dimensions with limits a=-Inf and b=Inf, which do not affect the probability but affect the computation time if they are not handled separately.

The estimated error is the standard error of the estimate, which decreases as one over the square root of the number of iterations. With tol, the iterations are added in batches, where each batch is sized to reach tol from the estimated error so far, and the estimates are averages over all batches. The number of iterations that is needed for a given accuracy depends strongly on the problem, and tol can therefore save a lot of computation time compared to a fixed number of iterations.

With size.tol, the target is instead the estimated error of the number of sub-integrals, from the last dimension, whose estimates are at least tol.level, relative to the number. This is the size of the excursion set in excursions(). The error is estimated by the error of the estimate where it goes below tol.level, divided by how fast the estimates decrease there. The target is at least half a sub-integral, since where the estimates pass tol.level is uncertain for any number of iterations. The batches of iterations after the first have at most 10000 iterations, which bounds the memory that is needed.

References

Bolin, D. and Lindgren, F. (2015) Excursion and contour uncertainty regions for latent Gaussian models, JRSS-series B, vol 77, no 1, pp 85-106.

Bolin, D. and Lindgren, F. (2018), Calculating Probabilistic Excursion Sets and Related Quantities Using excursions, Journal of Statistical Software, vol 86, no 1, pp 1-20.

Author

David Bolin davidbolin@gmail.com

Examples

## Create mean and a tridiagonal precision matrix
n <- 11
mu.x <- seq(-5, 5, length = n)
Q.x <- Matrix(toeplitz(c(1, -0.1, rep(0, n - 2))))
## Calculate the probability that the variable is between mu-3 and mu+3
prob <- gaussint(mu = mu.x, Q = Q.x, a = mu.x - 3, b = mu.x + 3, max.threads = 2)
prob$P
#> [1] 0.968005