Numerical evaluation of multinormal probabilities

The multivariate normal probability integral with a product correlation structure can be transformed to a one dimensional integral and easily evaluated when the correlation matrix is non-singular and well conditioned. However, the nearly singular case is much more difficult and previous methods fail to compute it with high numerical precision. This paper demonstrates that the (nearly) singular case can be computed to almost double precision using a three step adaptive Simpson method with the epsilon-algorithm by Wynn (1956). Tests using randomly chosen problems show that the method gives more reliable results than the adaptive Simpson method of Dunnett (1989) as well as the globally adaptive integration routine DQAGPE from QUADPACK
(Piessens et al., 1983).