| pmvnorm | R Documentation |
Computes the distribution function of the multivariate normal distribution for arbitrary limits and correlation matrices.
pmvnorm(lower=-Inf, upper=Inf, mean=rep(0, length(lower)),
corr=NULL, sigma=NULL, algorithm = GenzBretz(), keepAttr=TRUE,
seed = NULL, ...)
lower |
the vector of lower limits of length n. |
upper |
the vector of upper limits of length n. |
mean |
the mean vector of length n. |
corr |
the correlation matrix of dimension n. |
sigma |
the covariance matrix of dimension n less than 1000. Either |
algorithm |
an object of class |
keepAttr |
|
seed |
an object specifying if and how the random number generator
should be initialized, see |
... |
additional parameters (currently given to |
This program involves the computation of
multivariate normal probabilities with arbitrary correlation matrices.
It involves both the computation of singular and nonsingular
probabilities. The implemented methodology is described in
\bibcitetmvtnorm::numerical-:1992 and \bibcitetmvtnorm::comparison:1993 for algorithm
GenzBretz, in \bibcitetmvtnorm::Miwa+Hayter+Kuriki:2003 for algorithm
Miwa, useful up to dimension 20, and \bibcitetmvtnorm::Genz:2004
for the TVPACK algorithm, which covers 2- and 3-dimensional problems
for semi-infinite integration regions.
Note the default algorithm GenzBretz is randomized and hence slightly depends on
.Random.seed and that both -Inf and +Inf may
be specified in lower and upper. For more details see
pmvt.
The multivariate normal
case is treated as a special case of pmvt with df=0 and
univariate problems are passed to pnorm.
The multivariate normal density and random deviates are available using
dmvnorm and rmvnorm.
pmvnorm is based on original implementations by Alan Genz, Frank
Bretz, and Tetsuhisa Miwa developed for computing accurate approximations to
the normal integral. Users interested in computing log-likelihoods involving
such normal probabilities should consider function lpmvnorm,
which is more flexible and efficient for this task and comes with the
ability to evaluate score functions.
An overview is available from \bibcitetmvtnorm::Genz_Bretz_2009.
The evaluated distribution function is returned, if keepAttr is true, with attributes
error |
estimated absolute error |
msg |
status message(s). |
algorithm |
a |
*
qmvnorm for quantiles and lpmvnorm for
log-likelihoods.
n <- 5
mean <- rep(0, 5)
lower <- rep(-1, 5)
upper <- rep(3, 5)
corr <- diag(5)
corr[lower.tri(corr)] <- 0.5
corr[upper.tri(corr)] <- 0.5
prob <- pmvnorm(lower, upper, mean, corr)
print(prob)
stopifnot(pmvnorm(lower=-Inf, upper=3, mean=0, sigma=1) == pnorm(3))
a <- pmvnorm(lower=-Inf,upper=c(.3,.5),mean=c(2,4),diag(2))
stopifnot(round(a,16) == round(prod(pnorm(c(.3,.5),c(2,4))),16))
a <- pmvnorm(lower=-Inf,upper=c(.3,.5,1),mean=c(2,4,1),diag(3))
stopifnot(round(a,16) == round(prod(pnorm(c(.3,.5,1),c(2,4,1))),16))
# Example from R News paper (original by Genz, 1992):
m <- 3
sigma <- diag(3)
sigma[2,1] <- 3/5
sigma[3,1] <- 1/3
sigma[3,2] <- 11/15
pmvnorm(lower=rep(-Inf, m), upper=c(1,4,2), mean=rep(0, m), corr=sigma)
# Correlation and Covariance
a <- pmvnorm(lower=-Inf, upper=c(2,2), sigma = diag(2)*2)
b <- pmvnorm(lower=-Inf, upper=c(2,2)/sqrt(2), corr=diag(2))
stopifnot(all.equal(round(a,5) , round(b, 5)))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.