# mnorm: The multivariate normal distribution In mnormt: The Multivariate Normal and t Distributions, and Their Truncated Versions

## Description

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

## Usage

 ```1 2 3 4``` ```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, where `d=ncol(varcov)`, representing the coordinates of the point(s) where the density must be evaluated. `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`.

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`. If `d>2`, function `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. If `d=2`, a call to `sadmvn` is translated into one to `biv.nt.prob`; if `d=1`, `pnorm` is used.

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`.

## 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`

## Author(s)

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 (for `dmnormt`, `rmnormt`, etc.) 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.: Fortran77 code downloaded in 2005 and again in 2007 from his web-page, whose URL as of 2020-04-28 is http://www.math.wsu.edu/faculty/genz/software/software.html

`dnorm`, `dmt`, `biv.nt.prob`, `ptriv.nt`

## Examples

 ``` 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) ```

### Example output

``` 3.273202e-01 3.273202e-01 3.273202e-01 0.000000e+00 1.444171e-09
```

mnormt documentation built on Sept. 1, 2020, 5:09 p.m.