Introduction

A common problem in many statistics computations is computing the multivariate normal distribution function

$$P({\bf b}) = \frac{1}{\sqrt{\vert\Sigma\vert(2\pi)^m}} \int_... ...{b_m} e^{- \frac{1}{2}{\bf x}^t \Sigma^{-1} {\bf x}} d{\bf x},$$

where ${\bf x}= (x_1, x_2, ..., x_m)^t$, $\Sigma$ is an $m \times m$ symmetric positive definite covariance matrix and $b_i < \infty$ for all i.

There is reliable and efficient software available for computing $P({\bf b})$for m = 1 and m = 2 (for m = 2 see Donnelly, 1973 and Drezner and Wesolowsky, 1990), so assume m > 2. Perhaps the simplest method uses acceptance-rejection sampling, but this is not expected to be efficient for high accuracy work. Other methods for m > 2 use algorithms developed by Deák(1980, 1986 and 1990), Schervish (1984) and Genz (1992). The Schervish method has been compared with Genz's methods (Genz, 1992), but Deák's methods have not been compared with the other methods. Improvements to Genz's methods have recently been proposed by Beckers and Haegemans(1992), and Gibson, Glasbey and Elston (1992). The purpose of this paper is to present results from tests comparing acceptance-rejection sampling, Deák's methods, Schervish's methods and improved versions of Genz's methods. In Section 2 brief descriptions of the different methods are given, in Section 3 test results are reported and Section 4 provides some concluding remarks.