Monte Carlo Algorithm Tests

All of the methods described in the previous can be implemented as MC methods that use only $Uniform(0,1)$ random numbers. This includes the methods that use $\chi$ random numbers, which can be obtained from $Uniform(0,1)$ numbers via the inverse $\chi$ distribution, and the methods that use sequences of random orthogonal matrices, because each orthogonal matrix can be generated from a sequence of $\frac{m(m-1)}{2}$ $Uniform(0,1)$, (see Fang and Wang, 1994). Therefore all of the MC methods discussed so far, can be viewed as MC methods for integrals of the form

$$I(f) = \int_0^1\int_0^1 \ldots \int_0^1 f({\bf v}) d{\bf v},$$

where $f$ is appropriately chosen, and, for most of the algorithms, ${\bf v}$ has length $m$ or $m-1$, but could have length $m(m-1)/2$ for the methods that use random orthogonal matrices. The standard error can be used to provide an error estimate for all of the algorithms for these methods. To standardize notation, we use $\hat{\bf T}_N$ to denote an approximation to ${\bf T}$ obtained using $N$ simulation or random integrand values $\{f({\bf v}_k)\}_{k=1}^N$,with $\hat{\bf T}_N = \frac{1}{N}\sum_{k=1}^Nf({\bf v}_k)$. The standard error $\sigma_N^2$, for $\hat{\bf T}_N$, is defined using $\sigma_N^2 = \frac{1}{N(N-1)}\sum_{k=1}^N(f({\bf v}_k) - \hat{\bf T}_N)^2$. We will use the error estimate $\hat{\epsilon}= 3\sigma_N$ in order to provide an approximate confidence level of 99% for our algorithms.

Figure 3: Average Monte Carlo Symmetrized MVT Algorithm Times, $\epsilon = 10^{-2}$

When the input for a particular problem is given, there is usually no way of knowing how large $N$ needs to be before we have $3\sigma_N < \epsilon$, so some type of iterative algorithm is required. We have implemented all of our MC algorithms in the following iterative form.

Monte Carlo Algorithm

1. Input ${\bf a}$, ${\bf b}$, $\mbox{\boldmath$\Sigma$}$, $\nu$, $\mbox{\boldmath$\delta$}$, $\epsilon$, and $N_{max}$;.
2. Set $\bar{N}= N_{min}$;
3. Compute $\hat{\bf T}_{\bar{N}}$ and $\sigma_{\bar{N}}$;
4. Set $\bar{{\bf T}}=\hat{\bf T}_{\bar{N}}$, $\bar{\sigma}=\sigma_{\bar{N}}$, $\hat{\epsilon}= 3\bar{\sigma}$;
5. Do While ( $\hat{\epsilon}> \epsilon$ and $\bar{N}< N_{max}$ )
   1. Set $N=\max( \min( \lceil\bar{N}((\hat{\epsilon}/\epsilon)^2-1)\rceil, N_{max}-\bar{N}) , 25 )$;
   2. Compute $\hat{\bf T}_N$ and $\sigma_N$;
   3. Set $\bar{N}= \bar{N}+N$, $\bar{{\bf T}} = \bar{{\bf T}}+\bar{\sigma}^2(\hat{\bf T}_N-\bar{{\bf T}})/(\bar{\sigma}^2+\sigma_N^2)$, $\bar{\sigma}^2 = \bar{\sigma}^2\sigma_N^2/(\bar{\sigma}^2+\sigma_N^2)$, $\hat{\epsilon}= 3\bar{\sigma}$;
   End Do
6. Output $\bar{{\bf T}}$, $\bar{N}$, $\hat{\epsilon}$.

This is usually a two iteration algorithm. We initially use a relatively small $N_{min}$ (we set $N_{min} = 100$ for all MC algorithms) to estimate the variance $\sigma_f^2$ for $f$, assuming $\sigma_f \approx \sigma_{N_{min}}\sqrt{N_{min}}$. Then we try to pick an $N$ so that the next iteration completes the calculation. We choose the smallest $N$ satisfying $N+N_{min}\leq N_{max}$ and $3\bar{\sigma}= 3\sigma_{N+N_{min}} \approx 3\sigma_f/\sqrt{N+N_{min}} \leq \epsilon$. The algebraic rearrangement of this inequality provides the heuristic formula used at step (5a). In practice, we sometimes require more than two iterations, which is why we have implemented the algorithm in the form given. The updating formulas for $\bar{{\bf T}}$ and $\bar{\sigma}$ are algebraically equivalent to $\bar{{\bf T}} = \sum_{i=1}^k\frac{\hat{\bf T}_{N_i}}{\sigma_{N_i}^2}/ \sum_{i=1}^k\frac{1}{\sigma_{N_i}^2}$ and $\bar{\sigma}^2 = 1/\sum_{i=1}^k\frac{1}{\sigma_{N_i}^2},$ if the algorithm requires $k$ iterations, with sample sizes $N_1,N_2,\ldots,N_k$, so that the final $\bar{N}=\sum_{i=1}^kN_i$. This type of weighted combination of values for iterative MC algorithms was used by Lepage (1978). If $\sigma_{N_i}^2 \approx \sigma_f^2/N_i$, then $\bar{{\bf T}} \approx \sum_{i_1}^k\hat{\bf T}_{N_i}/\bar{N}$, and $\bar{\sigma}^2 \approx \sigma_f^2/\bar{N}$, which provides some motivation for the choice of weighting.

Figure 4: Average Monte Carlo Prioritized MVT Algorithm Times, $\epsilon = 10^{-2}$

The type of testing that we use consists of applying several algorithms to a set of randomly selected MVT problems, $\{ {\bf T}( {\bf a}_j, {\bf b}_j, \mbox{\boldmath$\Sigma$}_j, \nu_j ) \}_{j=1}^L$, for selected values of $\epsilon$, $m$ and $L$, and with $N_{max}$ chosen large enough to allow the algorithms to terminate. We generate random covariance matrices $\mbox{\boldmath$\Sigma$}_j$ which are random correlation matrices using the algorithm described by Marsaglia and Olkin (1984). We generate random limit vectors ${\bf a}_j$ with components $a_{i,j} = -3v_{i,j}\sqrt{m}$, ${\bf b}_j$ with components $b_{i,j} = 3w_{i,j}\sqrt{m}$, and random $\nu_j=\max(1,\lfloor 10u_j\sqrt{m}\rfloor)$, where all $v_{i,j}$, $w_{i,j}$ and $u_j$ are $Uniform(0,1)$ random numbers. The distributions that we have chosen for ${\bf b}_j$ and $\nu_j$ were picked after some experimentation to produce average MVT values in the range .2-.8 for $m$ in the range 2-20. We compute the average time taken by each method for each $m$, and the average, denoted by $\bar{c}$, for $\log(\vert{\bf T}_j - \hat{\bf T}_j\vert)/\log(\epsilon)$ (the ratio of the number of correct digits produced to the number of digits requested). The ``correct'' ${\bf T}_j$ values used for our tests were determined using our most efficient algorithm (QSMVN, described later in this section) with input error tolerance set at $\epsilon/10$. We also count the number of times $\hat{\epsilon}_j > \epsilon$ and the average of $\log_{10}(\vert{\bf T}_j - \hat{\bf T}_j\vert)-\log_{10}(\epsilon)$ (the number of wrong digits) for these events. A 433 MHZ DEC Alpha workstation was used for the computations, with all algorithms implemented in double precision in FORTRAN.

Figure 5: Average Monte Carlo MVT Algorithm Times, $\epsilon = 10^{-3}$

We first present some test results, for $m = 2, 3, \ldots, 20$, for three methods which will be denoted MCAR, MCSR, MCSVT. These methods use simple Monte Carlo algorithms based on equations (8), (9) and (12), respectively. The results, obtained using $L = 100$ samples, are given in Figure 2. The number given in parentheses next to the name of each method is the $\bar{c}$ statistic for that method, averaged over all of the $m$ values. All of the methods were reliable, terminating typically with a results that had approximately $2.8$ correct decimal digits. We also tested the algorithms based on equation (10, with $n=1$), denoted by MCSR1, and equation (14), denoted by MCSVN, but the results were similar to (with MCSR1 results slightly better than) the results for MCSVN at this accuracy level.

A standard method for improving for overall performance of MC algorithms is to use simple antithetic variates. This method replaces $f({\bf v}_k)$ by $(f({\bf v}_k)+f({\bf 1}-{\bf v}_k))/2$, with ${\bf 1}= (1,1,\ldots,1)'$, in the MC sum. In this case we redefine $\hat{\bf T}_N$ and $\sigma_N$ to be given by

$$\hat{\bf T}_N = \frac{2}{N}\sum_{k=1}^{N/2}\frac{f({\bf v}_k... ...rac{f({\bf v}_k)+f({\bf 1}-{\bf v}_k)}{2} - \hat{\bf T}_N)^2.$$

The methods that use equation (10) are already symmetrized. We tested the other algorithms, but with this ``symmetrized'' modification included, at the $\epsilon = 0.01$ accuracy level and the results (again with $L = 100$ samples) are given in Figure 3. The times are clearly lower for the symmetrized algorithms, with more significant improvements for the MCSR and MCSVN methods. The times for then MCAR method also showed some improvement, but they were still significantly higher than the times for the other methods. The clear loser at this accuracy level is the MCAR method (acceptance-rejection). We will not present any further results for this method, although we have carried out additional tests of MCAR and obtained similar results, when compared with other methods. We do not recommend this method for the efficient computation of MVT probabilities.

Figure 6: Average Quasi-Monte Carlo MVT Algorithm Times, $\epsilon = 10^{-3}$

For a final test at this accuracy level we added the variable prioritization preconditioning step, described near the end of Section 2.3, to the symmetrized algorithms. Some of the results (based on 100 samples) are given in Figure 4. The times are even lower for all of the prioritized algorithms, with the most significant improvement for the MCSVN method. All of the tests described in the rest of this paper used symmetrized algorithms with prioritization.

In order to more carefully compare the MC methods we carried out another test at a higher accuracy level, using the methods MCSVT, MCSVN, MCSR, MVSR1 and MCSR2 (based on equation (10) with $n = 2$). Some of the results are given in Figure 5. All of the SR test results appear to exhibit the overall increase in time taken as function of dimension, but with times for the odd $m$ larger than expected. This occurs because the $F$ values needed for these methods require the evaluation of a series (determined from a succession of integrations by parts) which ends with an extra $t_\nu$ value when $m$ is odd, and the computation of this $t_\nu$ value is what requires increases the time when $m$ is an odd integer.