gaussint {excursions}

Sequential estimation of Gaussian integrals

Description

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
)

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

Decides the number of threads the program can use. Set to 0 for using the maximum number of threads allowed by the system (default).

seed

The random seed to use (optional).

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.

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.

Author(s)

David Bolin [e-mail]

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.

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