The multivariate normal distribution

The probability density function, the distribution function and random number generation for a d-dimensional multivariate normal (Gaussian) random variable.

dmnorm(x, mean = rep(0, d), varcov, log = FALSE) 
pmnorm(x, mean = rep(0, d), varcov, ...) 
rmnorm(n = 1, mean = rep(0, d), varcov, sqrt=NULL) 
sadmvn(lower, upper, mean, varcov, maxpts = 2000*d, abseps = 1e-06, releps = 0)

Arguments

x: either a vector of length d or a matrix with d columns representing the coordinates of the point(s) where the density must be evaluated; see also ‘Details’ for restrictions on d.
mean: either a vector of length d, representing the mean value, or (except for rmnorm) a matrix whose rows represent different mean vectors; in the matrix case, only allowed for dmnorm and pmnorm, its dimensions must match those of x.
varcov: a symmetric positive-definite matrix representing the variance-covariance matrix of the distribution; a vector of length 1 is also allowed (in this case, d=1 is set).
sqrt: if not NULL (default value is NULL), a square root of the intended varcov matrix; see ‘Details’ for a full description.
log: a logical value (default value is FALSE); if TRUE, the logarithm of the density is computed.
...: arguments passed to sadmvn, among maxpts, abseps, releps.
n: the number of (pseudo) random vectors to be generated.
lower: a numeric vector of lower integration limits of the density function; must be of maximal length 20; +Inf and -Inf entries are allowed.
upper: a numeric vector of upper integration limits of the density function; must be of maximal length 20; +Inf and -Inf entries are allowed.
maxpts: the maximum number of function evaluations (default value: 2000*d).
abseps: absolute error tolerance (default value: 1e-6).
releps: relative error tolerance (default value: 0).

Details

The dimension d cannot exceed 20 for pmnorm and sadmvn. If this threshold is exceeded, NA is returned.

The function pmnorm works by making a suitable call to sadmvn if d>3, or to ptriv.nt if d=3, or to biv.nt.prob if d=2, or to pnorm if d=1. The R functions sadmvn, ptriv.nt and biv.nt.prob are, in essence, interfaces to underlying Fortran 77 routines by Alan Genz; see the references below. These routines use adaptive numerical quadrature and other non-random type techniques.

If sqrt=NULL (default value), the working of rmnorm involves computation of a square root of varcov via the Cholesky decomposition. If a non-NULL value of sqrt is supplied, it is assumed that it represents a matrix, \(R\) say, such that \(R' R\) represents the required variance-covariance matrix of the distribution; in this case, the argument varcov is ignored. This mechanism is intended primarily for use in a sequence of calls to rmnorm, all sampling from a distribution with fixed variance matrix; a suitable matrix sqrt can then be computed only once beforehand, avoiding that the same operation is repeated multiple times along the sequence of calls; see the examples below. Another use of sqrt is to supply a different form of square root of the variance-covariance matrix, in place of the Cholesky factor.

For efficiency reasons, rmnorm does not perform checks on the supplied arguments.

If, after setting the same seed value to set.seed, two calls to rmnorm are made with the same arguments except that one generates n1 vectors and the other n2 vectors, with n1<n2, then the n1 vectors of the first call coincide with the initial n2 vectors of the second call.

Value

dmnorm returns a vector of density values (possibly log-transformed); pmnorm returns a vector of probabilities, possibly with attributes on the accuracy in case x is a vector; sadmvn return a single probability with attributes giving details on the achieved accuracy; rmnorm returns a matrix of n rows of random vectors, or a vector in case n=1 or d=1.

References

Genz, A. (1992). Numerical Computation of multivariate normal probabilities. J. Computational and Graphical Statist., 1, 141-149.

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

Genz, A.: Fortran 77 code downloaded in 2005 and again in 2007 from his web-page, whose URL as of 2020-04-28 was https://www.math.wsu.edu/faculty/genz/software/software.html

Author

FORTRAN 77 code of SADMVN, package mvtdstpack.f, package tvpack and most auxiliary functions by Alan Genz; some additional auxiliary functions by people referred to within his programs; interface to R and additional R code (for dmnormt, rmnormt, etc.) by Adelchi Azzalini.

Note

The attributes error and status of the probability returned by pmnorm and sadmvn indicate whether the function had a normal termination, achieving the required accuracy. If this is not the case, re-run the function with a higher value of maxpts

Examples

x <- seq(-2, 4, length=21)
y <- cos(2*x) + 10
z <- x + sin(3*y) 
mu <- c(1,12,2)
Sigma <- matrix(c(1,2,0,2,5,0.5,0,0.5,3), 3, 3)
f <- dmnorm(cbind(x,y,z), mu, Sigma)
f0 <- dmnorm(mu, mu, Sigma)
p1 <- pmnorm(c(2,11,3), mu, Sigma)
p2 <- pmnorm(c(2,11,3), mu, Sigma, maxpts=10000, abseps=1e-10)
p <- pmnorm(cbind(x,y,z), mu, Sigma)
#
set.seed(123)
x1 <- rmnorm(5, mu, Sigma)
set.seed(123)
x2 <- rmnorm(5, mu, sqrt=chol(Sigma)) # x1=x2
eig <- eigen(Sigma, symmetric = TRUE)
R <- t(eig$vectors %*% diag(sqrt(eig$values)))
for(i in 1:50) x <- rmnorm(5, mu, sqrt=R)
#
p <- sadmvn(lower=c(2,11,3), upper=rep(Inf,3), mu, Sigma) # upper tail
#
p0 <- pmnorm(c(2,11), mu[1:2], Sigma[1:2,1:2])
p1 <- biv.nt.prob(0, lower=rep(-Inf,2), upper=c(2, 11), mu[1:2], Sigma[1:2,1:2])
p2 <- sadmvn(lower=rep(-Inf,2), upper=c(2, 11), mu[1:2], Sigma[1:2,1:2]) 
c(p0, p1, p2, p0-p1, p0-p2)
#> [1] 0.3273202 0.3273202 0.3273202 0.0000000 0.0000000
#
p1 <- pnorm(0, 1, 3)
p2 <- pmnorm(0, 1, 3^2)