The ‘New Numerical Algorithm for Multivariate Normal Probabilities in Package mvtnorm’ article from the 2009-1 issue.
(Miwa et al. 2003) proposed a numerical algorithm for evaluating multivariate normal probabilities. Starting with version 0.9-0 of the mvtnorm package (Hothorn et al. 2001; Genz et al. 2008), this algorithm is available to the R community. We give a brief introduction to Miwa’s procedure and compare it to a quasi-randomized Monte-Carlo procedure proposed by (Genz and Bretz 1999), which has been available through mvtnorm for some years now, both with respect to computing time and accuracy.
For low-dimensional problems, Miwa’s procedure is a numerical method
that it has advantages in achieving higher accuracy with less time
consumption compared to the Monte-Carlo method. For large dimension
An algorithm for calculating any non-centered orthant probability of a
non-singular multivariate normal distribution is described by
(Miwa et al. 2003). The probability function in a one-sided problem is
The distribution function for
The two-sided probability
A new algorithm
argument to pmvnorm()
and qmvnorm()
has been
introduced in mvtnorm
version 0.9-0 to switch between two algorithms: GenzBretz()
is the
default and triggers the use of the above mentioned quasi-randomized
Monte-Carlo procedure by (Genz and Bretz 1999). Alternatively,
algorithm = Miwa()
applies the procedure discussed here. Both
functions can be used to specify hyper-parameters of the algorithm. For
Miwa()
, the argument steps
defines the number of grid points
Algorithm | ||||
---|---|---|---|---|
Genz & Bretz |
0.08468833 | 0.001385620 | 0.008863600 | |
Genz & Bretz |
0.08472561 | 0.001390769 | 0.008863877 | 2.319286 |
Genz & Bretz |
0.08472682 | 0.001388424 | 0.008862195 | |
Miwa |
0.08472222 | 0.001388889 | 0.008863235 | |
Exact. | 0.08472222 | 0.001388889 | 0.008863236 |
The following example shows how to calculate the probability
> library("mvtnorm")
> m <- 3
> S <- diag(m)
> S[2, 1] <- S[1, 2] <- 1 / 4
> S[3, 1] <- S[3, 1] <- 1 / 5
> S[3, 2] <- S[3, 2] <- 1 / 3
> pmvnorm(lower = -c(1,4,2),
+ upper = c(1,4,2),
+ mean=rep(0, m), corr = S,
+ algorithm = Miwa())
[1] 0.6536804
attr(,"error")
[1] NA
attr(,"msg")
[1] "Normal Completion"
The upper limit and lower limit of the integral region are passed by the
vectors upper
and lower
. The mean vector and correlation matrix are
given by the vector mean
and the matrix corr
. From the result, we
know that
In this section, we compare the accuracy and time consumption of the R
implementation of the algorithm of (Miwa et al. 2003) with the default
procedure for approximating multivariate normal probabilities in
mvtnorm by (Genz and Bretz 1999).
The experiments were performed using an Intel
The exact value of
Both algorithms are linear in this simplest case and very fast (<0.01 second), so the time consumption is not discussed here.
When
the value of
We have implemented an R interface to the procedure of (Miwa et al. 2003) in
the R package mvtnorm.
For small dimensions, it is an alternative to quasi-randomized
Monte-Carlo procedures, which are computed by default. However, Miwa’s
algorithm has some disadvantages. When the dimension
The authors would like to thank the reviewers for their comments that help improve the manuscript and the package.
Algorithm | m=5 | m=9 | ||
---|---|---|---|---|
sec. | sec. | |||
Genz & Bretz |
0.1666398 | 0.029 | 0.09998728 | 0.231 |
Genz & Bretz |
0.1666719 | 0.132 | 0.09998277 | 0.403 |
Genz & Bretz |
0.1666686 | 0.133 | 0.09999726 | 0.431 |
Miwa |
0.1666667 | 0.021 | 0.09999995 | 89.921 |
Exact. | 0.1666667 | 0.10000000 |
Dimension | Miwa |
Genz & Bretz |
||
---|---|---|---|---|
One-sided | Two-sided | One-sided | Two-sided | |
0.021 | 0.441 | 0.029 | 0.085 | |
0.089 | 8.731 | 0.089 | 0.149 | |
0.599 | 156.01 | 0.083 | 0.255 | |
9.956 | 4hours | 0.138 | 0.233 | |
89.921 | - | 0.231 | 0.392 |
This article is converted from a Legacy LaTeX article using the texor package. The pdf version is the official version. To report a problem with the html, refer to CONTRIBUTE on the R Journal homepage.
Text and figures are licensed under Creative Commons Attribution CC BY 4.0. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".
For attribution, please cite this work as
Mi, et al., "New Numerical Algorithm for Multivariate Normal Probabilities in Package mvtnorm", The R Journal, 2009
BibTeX citation
@article{RJ-2009-001, author = {Mi, Xuefei and Miwa, Tetsuhisa and Hothorn, Torsten}, title = {New Numerical Algorithm for Multivariate Normal Probabilities in Package mvtnorm}, journal = {The R Journal}, year = {2009}, note = {https://rjournal.github.io/}, volume = {1}, issue = {1}, issn = {2073-4859}, pages = {37-39} }