mnorm | R Documentation |
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)
x |
either a vector of length |
mean |
either a vector of length |
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, |
sqrt |
if not |
log |
a logical value (default value is |
... |
arguments passed to |
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 |
upper |
a numeric vector of upper integration limits
of the density function; must be of maximal length |
maxpts |
the maximum number of function evaluations
(default value: |
abseps |
absolute error tolerance (default value: |
releps |
relative error tolerance (default value: |
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.
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
.
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
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.
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
dnorm
, dmt
,
biv.nt.prob
, ptriv.nt
,
plot_fxy
for plotting 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) # p1 <- pnorm(0, 1, 3) p2 <- pmnorm(0, 1, 3^2)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.