New Numerical Algorithm for Multivariate Normal Probabilities in Package mvtnorm

Abstract:

The ‘New Numerical Algorithm for Multivariate Normal Probabilities in Package mvtnorm’ article from the 2009-1 issue.

Cite PDF Tweet

Published

June 1, 2009

Citation

Mi, et al., 2009

Volume

Pages

1/1

37 - 39


proposed a numerical algorithm for evaluating multivariate normal probabilities. Starting with version 0.9-0 of the mvtnorm package , 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 , 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 m, Miwa’s procedure could get more accurate results, but the time consumption is huge compared to Genz & Bretz’ procedure. It is only applicable to problems with dimension smaller than 20, whereas the Monte-Carlo procedure by can be used to evaluate 1000-dimensional normal distributions. At the end of this article, a suggestion is given for choosing a suitable algorithm in different situations.

1 Introduction

An algorithm for calculating any non-centered orthant probability of a non-singular multivariate normal distribution is described by . The probability function in a one-sided problem is Pm(μ,R)=P{Xi0;1im}=00ϕm(x;μ,R)dx1dxm where μ=(μ1,...,μm) is the mean and R=(ρij) the correlation matrix of m multivariate normal distributed random variables X1,,XmNm(μ,R). The function ϕm denotes the density function of the m-dimensional normal distribution with mean μ and correlation matrix R.

The distribution function for c=(c1,,cm) can be expressed as: Fm(c)=P{Xici;1im}=P{Xici;1im}=Pm(μ+c,R). The m-dimensional non-centered orthant one-sided probability can be calculated from at most (m1)! non-centered probabilities with positive definite tri-diagonal correlation matrix. The algorithm for calculating such probabilities is a recursive linear integration procedure. The total order of a one-sided problem is G×m!, where G is the number of grid points for integration.

The two-sided probability Fm(d,c)=P{diXici;1im} can be calculated from 2m m-dimensional one-sided probabilities, which have the same mean and correlation matrix. The total order of this two-sided problem is G×m!×2m.

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 . 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 G to be evaluated.

Table 1: Value of probabilities with tri-diagonal correlation coefficients, ρi,i±1=ρ,1im and ρi,j=0,|ij|>1. ρ=21 or ρ=21.
Algorithm m=5 m=10
ρ=12 ρ=12 ρ=12 ρ=12
Genz & Bretz (ε=104) 0.08468833 0.001385620 0.008863600 2.376316×108
Genz & Bretz (ε=105) 0.08472561 0.001390769 0.008863877 2.319286 × 108
Genz & Bretz (ε=106) 0.08472682 0.001388424 0.008862195 2.671923×108
Miwa (G=128) 0.08472222 0.001388889 0.008863235 2.505205×108
Exact. 0.08472222 0.001388889 0.008863236 2.505211×108

The following example shows how to calculate the probability p=Fm(d,c)=P{1<X1<1,4<X2<4,2<X3<2}. with mean μ=(0,0,0) and correlation matrix R=(11/41/51/411/31/51/31) by using the following R code:

> 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 p=0.6536804 with given correlation matrix R.

2 Accuracy and time consumption

In this section, we compare the accuracy and time consumption of the R implementation of the algorithm of with the default procedure for approximating multivariate normal probabilities in mvtnorm by . The experiments were performed using an Intel® Pentium® processor with 2.8 GHz.

Probabilities with tri-diagonal correlation matrix

The exact value of Pm(μ,R) is known if R has some special structure. For example, when R is a m-dimensional tri-diagonal correlation matrix with correlation coefficients ρi,j={21j=i±10|ij|>11im the value of Pm(0,R) is ((1+m)!)1 . The accuracy of Miwa algorithm (G=128 grid points) and the Genz & Bretz algorithm (with absolute error tolerance ε=104,105,106) for probabilities with tri-diagonal correlation matrix are compared in Table 1. In each calculation, we have results with small variance. The values which do not hold the tolerance error are marked with bold characters and are underlined in the tables. When the dimension is larger than five, Genz & Bretz’ algorithm with error tolerance smaller than 105 is hard to achieve while Miwa’s algorithm with grid points G=128 achieves error tolerance smaller than 107.

Both algorithms are linear in this simplest case and very fast (<0.01 second), so the time consumption is not discussed here.

Centered orthant probabilities

When R is the correlation matrix with

ρi,j=21,ij,1im

the value of Pm(0,R) is known to be (1+m)1 . Accuracy and time consumption of Miwa’s algorithm and Genz & Bretz’ algorithm for this situation are compared in Table 2. As a numerical algorithm, Miwa’s algorithm still has better tolerance error. However, the time consumption of Miwa’s algorithm increases non-linearly when the dimension of the orthant probabilities increases. A detail time consumption analysis for both methods is given in Table 3. Miwa’s algorithm is much slower than Genz & Bretz’ algorithm in calculating two-sided orthant probability when the dimension m is larger than seven.

3 Conclusion

We have implemented an R interface to the procedure of 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 m increases, the time consumption of Miwa’s algorithm increases dramatically. Moreover, it can’t be applied to singular problems which are common in multiple testing problems, for example.

4 Acknowledgements

The authors would like to thank the reviewers for their comments that help improve the manuscript and the package.

Table 2: Accuracy and time consumption of centered orthant probabilities with correlation coefficients, ρi,j=21,ij, 1im.
Algorithm m=5 m=9
ρ=12 sec. ρ=12 sec.
Genz & Bretz (ε=104) 0.1666398 0.029 0.09998728 0.231
Genz & Bretz (ε=105) 0.1666719 0.132 0.09998277 0.403
Genz & Bretz (ε=106) 0.1666686 0.133 0.09999726 0.431
Miwa (G=128) 0.1666667 0.021 0.09999995 89.921
Exact. 0.1666667 0.10000000
Table 3: Time consumption of centered orthant probabilities (measured in seconds).
Dimension Miwa (G=128) Genz & Bretz (ε=104)
One-sided Two-sided One-sided Two-sided
m=5 0.021 0.441 0.029 0.085
m=6 0.089 8.731 0.089 0.149
m=7 0.599 156.01 0.083 0.255
m=8 9.956 4hours 0.138 0.233
m=9 89.921 - 0.231 0.392




CRAN packages used

mvtnorm

CRAN Task Views implied by cited packages

Distributions, Finance

Note

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.

Footnotes

    References

    A. Genz and F. Bretz. Numerical computation of multivariate t-probabilities with application to power calculation of multiple contrasts. Journal of Statistical Computation and Simulation, 63: 361–378, 1999.
    A. Genz, F. Bretz, T. Hothorn, T. Miwa, X. Mi, F. Leisch and F. Scheipl. Mvtnorm: Multivariate normal and t distribution. 2008. URL http://CRAN.R-project.org. R package version 0.9-0.
    T. Hothorn, F. Bretz and A. Genz. On multivariate t and Gauß probabilities in R. R News, 1(2): 27–29, 2001.
    T. Miwa, A. J. Hayter and S. Kuriki. The evaluation of general non-centred orthant probabilities. Journal of The Royal Statistical Society Series B–Statistical Methodology, 65: 223–234, 2003.

    Reuse

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

    Citation

    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}
    }