Quasi-Monte Carlo Algorithm Tests

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

Tests by Beckers and Haegemans (1992), and by Genz (1993), for MVN problems demonstrated that the performance of MC MVN methods could usually be improved if the sets of (pseudo-)random numbers used by the MC methods were replaced by appropriate sets of quasi-random numbers. If we wish to construct Quasi-Monte Carlo(QMC) MVT methods, we need only replace the $Uniform(0,1)$ random numbers in our MC methods by appropriately chosen sets of $Quasi(0,1)$ random numbers. All of our MC methods are implemented as methods that use only $Uniform(0,1)$ random numbers, However, simple QMC methods do not provide the statistically robust (standard)error estimates that MC methods do provide, so we decided to use randomized QMC algorithms. The QMC methods that we have implemented use approximations to $I(f)$ in the form

$$\hat{\bf T}_{N,P} = \frac{1}{N}\sum_{i=1}^N \frac{1}{2P}\sum... ... + f({\bf 1}- \vert 2\{{\bf p}_j+{\bf w}_i\}-{\bf 1}\vert) ).$$ (15)

In this definition, $\{{\bf x}\}$ denotes the vector obtained by taking the fractional part of each of the components of ${\bf x}$, $w_{k,i} \sim Uniform(0,1)$ and ${\bf p}_j, j = 1, 2, \ldots, P$ is a set of quasi-random points. If the dimension of ${\bf v}$ is $m$ or $m-1$, we use good lattice point sets (see Sloan and Joe, 1994, and Hickernell, 1998) for the required quasi-random point sets. If the dimension of ${\bf v}$ is $\frac{m(m-1)}{2}$, we use ``Richtmeyer sequence'' point sets (see Davis and Rabinowitz, 1984, pp. 482-483) for the dimensions greater than 100. The Richtmeyer points are defined by $p_{k,j} = \{j\sqrt{q_{k-100}}\}$, for $k=101, 102, \ldots,\frac{m(m-1)}{2}$, $j = 1, 2, \ldots, P$, where $q_i$ is the $i^{th}$ prime number (starting with $q_1 =2$). The ``periodizing'' transformation $\vert 2{\bf v}- {\bf 1}\vert$ is included because the quasi-Monte Carlo rules that we use have better convergence properties for periodic integrands. If we let $Q_P({\bf w}) = \frac{1}{2P}\sum_{j=1}^P ( f(\vert 2\{{\bf p}_j+{\bf w}\}-{\bf 1}\vert) + f({\bf 1}- \vert 2\{{\bf p}_j+{\bf w}\}-{\bf 1}\vert) )$, then the standard error $\sigma_N$ for the approximation (15) can be determined using $\sigma_N^2 = \frac{1}{N(N-1)}\sum_{i=1}^N(Q_P({\bf w}_i)-\hat{\bf T}_{N,P})^2$. We have implemented the quasi-Monte Carlo methods using the following algorithm.

Quasi-Monte Carlo Algorithm

1. Input ${\bf a}$, ${\bf b}$, $\mbox{\boldmath$\Sigma$}$, $\nu$, $\mbox{\boldmath$\delta$}$, $\epsilon$, and $N_{max}$;.
2. Set $N = N_{min}$, $i=0$;
3. Compute $\hat{\bf T}_{\bar{N}}$ and $\sigma_{\bar{N}}$;
4. Set $\bar{N}=2NP_0$, $\bar{{\bf T}}=\hat{\bf T}_{N,P_0}$, $\bar{\sigma}=\sigma_N$, $\hat{\epsilon}= 3\bar{\sigma}$;
5. Do While ( $\hat{\epsilon}> \epsilon$ and $\bar{N}< N_{max}$ )
   1. Set $i = i+1$;
   2. Compute $\hat{\bf T}_{N,P_i}$ and $\sigma_N$.
   3. Set $\bar{N}= \bar{N}+2NP_i$, $\bar{{\bf T}} = \bar{{\bf T}}+\bar{\sigma}^2(\hat{\bf T}_{N,P_i}-\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.5\bar{\sigma}$;
   End Do
6. Output $\bar{\hat{\bf T}}$, $\bar{N}$, $\hat{\epsilon}$.

The sequence $P_0, P_1, \ldots$ is a sequence of primes starting with $P_0 = 31$, with $P_{i+1} \approx 3P_i/2$. All of our QMC algorithms use $N_{min} = 8$, with $\hat{\epsilon}= 3.5\bar{\sigma}$ (using 3.5 instead of 3 because of the smaller $N_{min}$).

Figure 8: Average MVT Algorithm Times, $\epsilon = 10^{-4}$

We first tested the quasi-Monte Carlo algorithms (with prioritization included) at the $\epsilon = 10^{-3}$ accuracy level. Some of the results (based on 100 samples) are given in Figure 6, with QRSVN results omitted because they were similar to QRSVT results. The times are significantly lower than the MC times for the QRSR, QRSVT and QRSVN algorithms, with clear winners QRSVT and QRSVN. We believe that there was no significant improvement in the QRSR1 and QRSR2 algorithm times, compared to MCSR1 and MCSR2 times because the MCSR1 and MCSR2 algorithms can themselves be considered randomized quasi-random algorithms (based on the $S_n(Z)$ rules which use evenly distributed spherical surface points), so the additional ``quasi-randomization'' of the $Z$ matrices does not produce a significant improvement in algorithm performance.

We conducted an additional test of the quasi-Monte Carlo algorithms QRSVT and QRSVN at the $\epsilon = 10^{-4}$ accuracy level. The results (based on 100 samples) are given in Figure 7. The times are significantly lower for the QRSVN algorithm. We were surprised by this result, because the SVN method is based on replacing an $m$-dimensional problem with an $m$+1-dimensional problem. An explanation for this difference comes from the fact that each QRSVN $f$ value requires $m$ $\Phi$ values, $m-1$ $\Phi^{-1}$ values, and one $\chi^{-1}$ value, but each QRSVT $f$ value requires $m$ $t$ values and $m-1$ $t^{-1}$ values. The $t$ values, using $\nu, \nu+1, \ldots, \nu+m-1$ degrees of freedom, are computed using integration by parts so that $t_{k}$ uses a sum of $k/2$ terms. Some $t$ values are also used in the $t^{-1}$ evaluations, so the the work for the 1-dimensional distribution evaluation for one $f$ value for QRSVT is $O(m^2)$. For the QRSVN method, the $\chi^{-1}$ time is $O(\nu)$, but the individual $\Phi$ and $\Phi^{-1}$ times are independent of $m$, so the time one $f$ value for the QRSVN method is only $O(m+\nu)$. We think these $f$ value time complexity differences explain the increasing difference between the SVN and SVT times as $m$ increases.

We also conducted some tests with some subregion adaptive algorithms, SASVN and SASVT, at the $\epsilon = 10^{-4}$ accuracy level. These algorithms use a subregion adaptive integration method, similar to the one that was effective for the lower dimensional MVN problems (see Genz, 1992, 1993, and Berntsen, Espelid and Genz, 1991), applied to the respective SV-Chi-Normal and SV-t formulations of the MVT problem. The results (based on 100 samples, for $m =2, \ldots, 11$) are given in Figure 8. The times for the SASVN are significantly lower than the QRSVN for dimensions 2-8, but after that the SASVN times increase rapidly, usually exceeding the QRSVN times for $m > 10$. The SASVT times (not shown in Figure 8) exhibited similar behavior, but they were consistently larger than the SASVN times, usually exceeding the QRSVN times for $m> 8$.

These results provide strong evidence that multivariate t-probabilities can be robustly and reliably computed at low to moderate accuracy levels in less than a second of workstation time for problems with up to twenty dimensions. The symmetrization, variable prioritization, and quasi-randomization techniques all produce clearly observable improvements in algorithm performance. When moderate accuracy is required, the QRSVN method can be significantly faster than the other methods, except for $m < 10$, where the SASVN method can be faster. Software for all of the methods discussed here is available from the authors.