R: Multivariate normal distribution

dmnorm {mnormt}

R Documentation

Multivariate normal distribution

Description

The probability density function, the distribution function and random number generation for the multivariate normal (Gaussian) probability distribution

Usage

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

Arguments

`x`	for `dmnorm`, this is either a vector of length `d` or a matrix with `d` columns, where `d=ncol(varcov)`, giving the coordinates of the point(s) where the density must be evaluated; for `pmnorm`, only a vector of length `d` is allowed, and `d` cannot exceed 20
`mean`	a numeric vector representing the expected value of the distribution; it must be of length `d`, as defined above. For `dmnorm` it can be a matrix; in this case its dimensions must match those of `x`
`varcov`	a 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)
`log`	a logical value (default value is `FALSE`); if `TRUE`, the logarithm of the density is computed
`...`	parameters passed to `sadmvn`, among `maxpts`, `abseps`, `releps`
`n`	the number of random numbers 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 function pmnorm works by making a suitable call to sadmvn if d>2, or to biv.nt.prob if d=2, or to pnorm if d=1. Function sadmvn is an interface to a Fortran-77 routine with the same name written by Alan Genz, and available from his web page; this makes uses of some auxiliary functions whose authors are documented in the Fortran code. The routine uses an adaptive integration method.

Value

dmnorm returns a vector of density values (possibly log-transformed); pmnorm and sadmvn return a single probability with attributes giving details on the achieved accuracy; rmnorm returns a matrix of n rows of random vectors

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 an higher value of maxpts

Author(s)

Fortran code of SADMVN and most auxiliary functions by Alan Genz, some additional auxiliary functions by people referred to within his program. Porting to R and additional R code by Adelchi Azzalini

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 code available at http://www.math.wsu.edu/math/faculty/genz/software/fort77/mvn.f

Examples

x <- seq(-2,4,length=21)
y <- 2*x+10
z <- x+cos(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)
p1 <- pmnorm(c(2,11,3), mu, Sigma)
p2 <- pmnorm(c(2,11,3), mu, Sigma, maxpts=10000, abseps=1e-10)
x <- rmnorm(10, mu, Sigma)
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)

[Package mnormt version 1.4-7 Index]