The probability density function, the distribution function and random number generation for the multivariate normal (Gaussian) distribution
1 2 3 4
either a vector of length
either a vector of length
a symmetric positive-definite matrix representing the
variance-covariance matrix of the distribution;
a vector of length 1 is also allowed (in this case,
a logical value (default value is
parameters passed to
the number of random vectors to be generated.
a numeric vector of lower integration limits of
the density function; must be of maximal length
a numeric vector of upper integration limits
of the density function; must be of maximal length
the maximum number of function evaluations
absolute error tolerance (default value:
relative error tolerance (default value:
pmnorm works by making a suitable call to
d>2, or to
sadmvn is an interface to a Fortran-77 routine with
the same name written by Alan Genz, available from his web page,
which works using an adaptive integration method.
This Fortran-77 routine makes uses of some auxiliary functions whose authors
are documented in the code.
sqrt=NULL (default value), the working of
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
If, after setting the same seed value to
two calls to
rmnorm are made with the same arguments except that one
n1 vectors and the other
n2 vectors, with
n1<n2, then the
n1 vectors of the first call coincide with the
n2 vectors of the second call.
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
status of the probability
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
Fortran code of
SADMVN and most auxiliary functions by Alan Genz,
some additional auxiliary functions by people referred to within his
program. Interface to R and additional R code by Adelchi Azzalini
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 code available at http://www.math.wsu.edu/math/faculty/genz/software/fort77/mvn.f
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28
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) # p1 <- pnorm(0, 1, 3) p2 <- pmnorm(0, 1, 3^2)
 3.273202e-01 3.273202e-01 3.273202e-01 0.000000e+00 1.444171e-09
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.