Introduction

A problem that arises in many statistics applications is that of computing the multivariate normal distribution function

$$F({\bf a, b}) = \frac{1}{\sqrt{\vert\Sigma\vert (2\pi)^m}} ... ...ma^{-1} \mbox{\boldmath $\theta$}} d\mbox{\boldmath $\theta$},$$

where $\mbox{\boldmath$\theta$}= (\theta_1, \theta_2, ..., \theta_m)^t$ and $\Sigma$ is an $m \times m$ symmetric positive definite covariance matrix. If, for some i, $a_i$ is $-\infty$ and $b_i$ is $\infty$, an appropriate transformation allows the ith variable to be integrated explicitly and reduces the number of variables in the problem, so, for each i, at least one of $a_i$ and $b_i$ is assumed to be finite.

There is reliable and efficient software available for the cases $m = 1$ and $m = 2$ (for $m = 2$ see [DON 73], [DW 90] and [CW 91]), so assume $m > 2$. An algorithm by Schervish [SCHRV 84] is implemented for $m < 8$. This algorithm uses a locally adaptive numerical integration strategy based on Simpson's rule. It also provides an error bound for $F$. Tests reported in Section 5 of this article, however, show that this algorithm can require very long computation times for $m > 4$.

An obvious approach to computing F is to try currently available numerical integration software (e.g Monte-Carlo methods or subregion adaptive methods) directly. However, the infinite integration limits need to be handled either by using some type of transformation to a finite region or, by using carefully selected cutoff values. In either case, any numerical integration method that does not take account of the peaked character of the integrand can waste significant amounts of computational effort for many problems.

The purpose of this article is to show that a relatively straightforward sequence of transformations can be used to put the problem into a form that allows reliable and efficient computation of $F$ using either simple Monte-Carlo or subregion adaptive numerical integration algorithms. These transformations will be described in the Section 2. Section 3 gives a complete description of a simple Monte-Carlo algorithm and Section 4 provides an illustrative example. Test results comparing the Schervish algorithm with Monte-Carlo and subregion adaptive algorithms that use the transformations will be presented in the Section 5.