An Algorithm

The problem is now in a form that could be presented as input for a variety of different numerical integration algorithms. Test results which will be reported in the next section show that a simple Monte-Carlo algorithm is surprisingly effective so the details of such an algorithm, which incorporates the transformations discussed in the previous section, are now given.

Multivariate Normal Probabilities Algorithm

1.

Input $\Sigma$, ${\bf a}$, ${\bf b}$, $\epsilon$, $\alpha$ and $N_{max}$.

2.

Compute lower triangular Cholesky factor $C$ for $\Sigma$.

3.

Initialize $Intsum = 0, N = 0, Varsum = 0$,
$d_1 = \Phi(a_1/c_{1,1}), e_1 = \Phi(b_1/c_{1,1})$ and $f_1 = e_1 - d_1$.

4.

Repeat

1. Generate uniform random $w_1, w_2, ..., w_{m-1} \epsilon [0,1]$.

2. For $i = 2, 3, ..., m$ set $y_{i-1} = \Phi^{-1}(d_{i-1}+w_{i-1}(e_{i-1}-d_{i-1}))$,
   $d_i = \Phi((a_i-\sum_{j=1}^{i-1}c_{i,j}y_j)/c_{i,i})$, $e_i = \Phi((b_i-\sum_{j=1}^{i-1}c_{i,j}y_j)/c_{i,i})$
   and $f_i = (e_i-d_i)f_{i-1}$.
3. Set $N = N+1, \delta = ( f_m - Intsum )/N, Intsum = Intsum + \delta,\\ Varsum = (N-2)Varsum/N + \delta^2$ and $Error = \alpha\sqrt{Varsum}$.

Until $Error < \epsilon$ or $N = N_{max}$.

5.

Output $F = Intsum, Error$, and $N$.

The input parameter $\alpha$ is the usual Monte-Carlo confidence factor for the standard error. If for example, $\alpha = 2.5$ is used, we expect the actual error in $F$ to be less than $Error$, 99% of the time. $N_{max}$ is an input parameter that limits the total amount of time allowed for the computation.

For problems where for some i's either $a_i$ is $-\infty$ or $b_i$ is $\infty$, any implementation of the algorithm should explicitly set $d_i = 0$ or $e_i = 1$ to avoid wasteful evaluation of $\Phi$. Preliminary tests with the algorithm showed that there was usually some reduction in the computation time if the variables were reordered (along with appropriate rows and columns of $\Sigma$) so that the variables associated with the largest integration intervals were the innermost variables (this sorting of the variables was suggested by Schervish [SCHRV 84]). Bexause this reordering operation does not take much time compared to the total computation time it was added to the algorithm as an additional step between steps 1 and 2 and used for the tests described in the Section 5.

In some applications it is necessary to compute integrals in the form

$$G(g) = \frac{1}{\sqrt{\vert\Sigma\vert (2\pi)^m}} \int_{a_1... ...ta$}} g(\mbox{\boldmath $\theta$}) d\mbox{\boldmath $\theta$},$$

where $g(\mbox{\boldmath$\theta$})$ might be one or more application specific expectation functions. For these integrals only minor modifications to the algorithm are needed. Because $g(\mbox{\boldmath$\theta$})$ might depend on $\theta_m$, $w_m$ would need to be generated as part of step 4(a) and $y_m = \Phi^{-1}(d_m+w_m(e_m-d_m))$ would need to be added to step 4(b), along with an additional substep to explicitly compute $\mbox{\boldmath$\theta$}= C{\bf y}$. Finally, in step 4(c), $f_m$ would need to be replaced by $f_mg(\mbox{\boldmath$\theta$})$. If the application requires an integral value for more than one $g(\mbox{\boldmath$\theta$})$, then all of the $g$ functions should be done in the same calculation by further expanding steps 3, 4(c) and the error test, in order to avoid duplicating the work required for computing the transformations.