Last chance! 50% off unlimited learning
Sale ends in
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, ...)The evaluated distribution function is returned, if keepAttr is true, with attributes
estimated absolute error
status message(s).
a character string with class(algorithm).
the vector of lower limits of length n.
the vector of upper limits of length n.
the mean vector of length n.
the correlation matrix of dimension n.
the covariance matrix of dimension n less than 1000. Either corr or
sigma can be specified. If sigma is given, the
problem is standardized internally. If corr is given,
it is assumed that appropriate standardization was performed
by the user. If neither corr nor
sigma is given, the identity matrix is used
for sigma.
an object of class GenzBretz,
Miwa or TVPACK
specifying both the algorithm to be used as well as
the associated hyper parameters.
logical indicating if
attributes such as error and msg should be
attached to the return value. The default, TRUE is back compatible.
an object specifying if and how the random number generator
should be initialized, see simulate.
additional parameters (currently given to GenzBretz for
backward compatibility issues).
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 Genz (1992, 1993) (for algorithm GenzBretz), in Miwa et al. (2003) for algorithm Miwa (useful up to dimension 20) and 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.
Genz, A. (1992). Numerical computation of multivariate normal probabilities. Journal of Computational and Graphical Statistics, 1, 141--150.
Genz, A. (1993). Comparison of methods for the computation of multivariate normal probabilities. Computing Science and Statistics, 25, 400--405.
Genz, A. (2004), Numerical computation of rectangular bivariate and trivariate normal and t-probabilities, Statistics and Computing, 14, 251--260.
Genz, A. and Bretz, F. (2009), Computation of Multivariate Normal and t Probabilities. Lecture Notes in Statistics, Vol. 195. Springer-Verlag, Heidelberg.
Miwa, T., Hayter J. and Kuriki, S. (2003). The evaluation of general non-centred orthant probabilities. Journal of the Royal Statistical Society, Ser. B, 65, 223--234.
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)))
Run the code above in your browser using DataLab