New Numerical Algorithm for Multivariate Normal Probabilities in Package mvtnorm

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 \(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 (Genz and Bretz 1999) 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 (Miwa et al. 2003). The probability function in a one-sided problem is \[\begin{aligned} P_{m}(\mu, \bf{R}) & = & \mathbb{P}\{X_i \geq 0; \, 1\leq i\leq m\} \\ & = & \int_0^\infty \dots \int_0^\infty \phi_m(\mathbf{x}; \mu, \mathbf{R}) \, dx_1 \dots dx_m \end{aligned}\] where \(\mu =(\mu_1,...,\mu_m)^\top\) is the mean and \(\mathbf{R}= (\rho_{ij})\) the correlation matrix of \(m\) multivariate normal distributed random variables \(X_1, \dots, X_m \sim \mathcal{N}_m(\mu, \mathbf{R})\). The function \(\phi_m\) denotes the density function of the \(m\)-dimensional normal distribution with mean \(\mu\) and correlation matrix \(\mathbf{R}\).

The distribution function for \(\mathbf{c}= (c_1, \dots, c_m)\) can be expressed as: \[\begin{aligned} F_m(\mathbf{c}) & = & \mathbb{P}\{X_i \leq c_i; \, 1\leq i\leq m \} \\ & = & \mathbb{P}\{-X_i \geq - c_i; \, 1\leq i\leq m \} \\ & = & P_m(-\mu+\mathbf{c},\mathbf{R}). \end{aligned}\] The \(m\)-dimensional non-centered orthant one-sided probability can be calculated from at most \((m-1)!\) 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 \times m!\), where \(G\) is the number of grid points for integration.

The two-sided probability \[\begin{aligned} F_{m}(\mathbf{d}, \mathbf{c}) = \mathbb{P}\{d_i \leq X_i \leq c_i; \, 1\leq i\leq m\} \end{aligned}\] can be calculated from \(2^m\) \(m\)-dimensional one-sided probabilities, which have the same mean and correlation matrix. The total order of this two-sided problem is \(G \times m! \times 2^m\).

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

Table 1: Value of probabilities with tri-diagonal correlation coefficients, \(\rho_{i,i\pm1}=\rho, 1\leq i \leq m\) and \(\rho_{i,j}=0,\forall |i-j|>1\). \(\rho = 2^{-1}\) or \(\rho= -2^{-1}\).

Algorithm	\(m=5\)		\(m=10\)
	\(\rho=\frac{1}{2}\)	\(\rho=-\frac{1}{2}\)	\(\rho=\frac{1}{2}\)	\(\rho=-\frac{1}{2}\)
Genz & Bretz \((\varepsilon=10^{-4})\)	0.08468833	0.001385620	0.008863600	\(2.376316 \times 10^{-8}\)
Genz & Bretz \((\varepsilon=10^{-5})\)	0.08472561	0.001390769	0.008863877	2.319286 \(\times\) \(10^{-8}\)
Genz & Bretz \((\varepsilon=10^{-6})\)	0.08472682	0.001388424	0.008862195	\(2.671923 \times 10^{-8}\)
Miwa \((G=128)\)	0.08472222	0.001388889	0.008863235	\(2.505205 \times 10^{-8}\)
Exact.	0.08472222	0.001388889	0.008863236	\(2.505211 \times 10^{-8}\)

The following example shows how to calculate the probability \[\begin{aligned} p & = & F_{m}(\mathbf{d}, \mathbf{c}) \\ & = & \mathbb{P}\{-1<X_1<1,-4<X_2<4,-2<X_3<2\}. \end{aligned}\] with mean \(\mu = (0,0,0)^\top\) and correlation matrix \[\begin{aligned} \mathbf{R}= \left( \begin{array}{ccc} 1 & 1/4 & 1/5 \\ 1/4 & 1 & 1/3 \\ 1/5 & 1/3 & 1 \end{array} \right) \end{aligned}\] 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 \(\mathbf{R}\).

2 Accuracy and time consumption

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\(\circledR\) Pentium\(\circledR\) processor with 2.8 GHz.

Probabilities with tri-diagonal correlation matrix

The exact value of \(P_m(\mu, \mathbf{R})\) is known if \(\mathbf{R}\) has some special structure. For example, when \(\mathbf{R}\) is a \(m\)-dimensional tri-diagonal correlation matrix with correlation coefficients \[\begin{aligned} \rho_{i,j} = \left\{ \begin{array}{cc} -2^{-1} & j = i \pm 1 \\ 0 & |i-j|>1 \end{array} \right. 1 \leq i \leq m \end{aligned}\] the value of \(P_m(\mathbf{0}, \mathbf{R})\) is \(((1+m)!)^{-1}\) (Miwa et al. 2003). The accuracy of Miwa algorithm (\(G = 128\) grid points) and the Genz & Bretz algorithm (with absolute error tolerance \(\varepsilon=10^{-4},10^{-5},10^{-6}\)) 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 \(10^{-5}\) is hard to achieve while Miwa’s algorithm with grid points \(G=128\) achieves error tolerance smaller than \(10^{-7}\).

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 \(\mathbf{R}\) is the correlation matrix with

\[\begin{aligned} \rho_{i,j} = 2^{-1}, \, i\neq j, \, 1 \leq i \leq m\\ \end{aligned}\]

the value of \(P_m(\mathbf{0}, \mathbf{R})\) is known to be \((1+m)^{-1}\) (Miwa et al. 2003). 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 (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 \(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, \(\rho_{i,j}=2^{-1}, i\neq j\), \(1 \leq i \leq m\).

Algorithm	m=5		m=9
	\(\rho=\frac{1}{2}\)	sec.	\(\rho=\frac{1}{2}\)	sec.
Genz & Bretz \((\varepsilon=10^{-4})\)	0.1666398	0.029	0.09998728	0.231
Genz & Bretz \((\varepsilon=10^{-5})\)	0.1666719	0.132	0.09998277	0.403
Genz & Bretz \((\varepsilon=10^{-6})\)	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 \((\varepsilon=10^{-4})\)
	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.