R Graphical Manual

Browse All

Last data update: 2014.03.03

R: Sequential estimation of Gaussian integrals

gaussint

R Documentation

Sequential estimation of Gaussian integrals

Description

gaussint is used for calculating Gaussian integrals

int_a^b frac{|Q|^{1/2}}{(2π)^{n/2}} exp(-frac1{2}(x-μ)^{T}Q(x-μ)) 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,
         LDL=FALSE)

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).
`LDL`	Use LDL factorisations? This is useful for analysis of problems with positive semidefinite precisions.

Value

A list:

`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 davidbolin@gmail.com

References

Bolin, D. and Lindgren, F. (2013) Excursion and contour uncertainty regions for latent Gaussian models, Journal of the Royal Statistical Society: Series B (in press).

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 process 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