mvtnorm (version 0.3-5)

pmvnorm: Multivariate Normal Distribution

Description

Computes the distribution function of the multivariate normal distribution for arbitrary limits and correlation matrices based on algorithms by Genz and Bretz.

Usage

pmvnorm(lower, upper, mean, corr, maxpts = 25000, abseps =0.001, releps = 0)

Arguments

lower

the vector of lower limits of length n.

upper

the vector of upper limits of length n.

mean

the mean vector of length n.

corr

the correlation matrix of dimension n.

maxpts

maximum number of function values as integer.

abseps

absolute error tolerance as double.

releps

relative error tolerance as double.

Value

• A list with the following components:
• value
  estimated integral value.
• error
  estimated absolute error.
• msg
  status messages.

Details

This program involves the computation of multivariate normal-probabilities with arbitrary correlation matrices. It involves both the computation of singular and nonsingular probabilities. The methodology is described in Genz (1992, 1993).

Note that both -Inf and +Inf may be specified in lower and upper. For more details see pmvt.

The mvn case is treated as a special case of pmvt with df=0.

Univariate problems are passed to pnorm.

Multivariate normal density and random numbers are available using dmvnorm (available in Mvnorm) and mvrnorm.

References

Genz, A. (1992). Numerical computation of multivariate normal probabilities. Journal of Computational and Graphical Statistics, 1, 141--150

Genz, A. (1993). Comparison of methods for the computation of multivariate normal probabilities. Computing Science and Statistics, 25, 400--405

Examples

n <- 5
mean <- rep(0, 5)
lower <- rep(-1, 5)
upper <- rep(3, 5)
corr <- diag(5)
corr[lower.tri(corr)] <- 0.5
corr[upper.tri(corr)] <- 0.5
prob <- pmvnorm(lower, upper, mean, corr)
print(prob)

stopifnot(pmvnorm(lower=-Inf, upper=3, mean=0, corr=1)$value == pnorm(3))

a <- pmvnorm(lower=rep(-Inf,2),upper=c(.3,.5),mean=c(2,4),diag(2))$value

stopifnot(round(a,16) == round(prod(pnorm(c(.3,.5),c(2,4))),16))

a <- pmvnorm(lower=rep(-Inf,3),upper=c(.3,.5,1),mean=c(2,4,1),diag(3))$value

stopifnot(round(a,16) == round(prod(pnorm(c(.3,.5,1),c(2,4,1))),16))

# Example from R News paper (original by Genz, 1992):

m <- 3
sigma <- diag(3)
sigma[2,1] <- 3/5
sigma[3,1] <- 1/3
sigma[3,2] <- 11/15
pmvnorm(lower=rep(-Inf, m), upper=c(1,4,2), mean=rep(0, m), sigma)