TVN Algorithm Testing

The testing of the algorithms requires a source of highly accurate TVN probabilities was needed for arbitrary ${\bf b}$ and $R$. Initially, Schervish's MULNOR was tried, but a number of inconsistencies were found, including some clear discrepancies between some exactly known results and the output from MULNOR. Kennedy and Wang (1992) also reported erroneous results from MULNOR, although the version of MULNOR used for these tests (emailed from statlib, and supported by the statlib $\Phi(x)$) did not fail for those TVN cases reported as failures by Kennedy and Wang. All tests reported in the rest of this section used an algorithm based on equation (14), implemented in quadruple precision with an adaptive numerical integration algorithm (requested absolute accuracy set at $10^{-18}$) combined with a quadruple precision BVN algorithm based on Owen's (1956) series method, for comparison with single and double precision results.

In order to carry out careful tests of the different algorithms, a large set of general problems was needed with variation in all of the six parameters in ${\bf b}$ and $R$. Let $R = CC^T$, with

$$C = \left[ \begin{array}{ccc} 1 & 0 & 0 \\ \cos(\pi\theta... ...sin(\pi\theta_3) & \sin(\pi\theta_2) \end{array} \right] .$$ (15)

The first tests that were completed used $\theta_1 = 1/258, 17/258, ..., 257/258$, $\theta_2 = 1/258, 17/258, ..., 257/258$, $\theta_3 = 1/258, 17/258, ..., 257/258$, $b_1 = -5, -4, ..., 5$, $b_2 = -5, -4, ..., 5$, $b_3 = b_2, b_2+1, ..., 5$. A complete test with all of the different combinations of parameters checks more than 2 million values for $\Phi({\bf b}, R )$, and includes problems with $\vert R\vert$ as small as $2\cdot10^{-8}$ The methods using equations (10-11), (13-14) were tested using an integration rule with a fixed number, K, of integration rule points for K = 6, 12, 24, and the results are given in the first three rows of Table 1.

The Table 1 (first three row) results are consistent with results reported by Gassmann. The Plackett formula methods (13-14) were the most accurate. None of the TVN tests described so far used $b$ values that were approximately equal. This type of limit combination was what motivated the final modifications described in the previous section for the Drezner-Wesolowsky BVN algorithm. In order to investigate the sensitivity of the TVN algorithms to this type of problems, another test was completed. This test used the same $\theta$'s that were used for the Table 1 first three row results, but used $b_1 = -5, -4, ..., 5$, $b_2 = -5+\epsilon, -4+\epsilon, ..., 5+\epsilon$, $b_3 = b_2, b_2+1, ..., 5$, with $\epsilon = 10^{-2}$. The 24-point rule results from this test are reported in Table 1, where it can be seen that the Drezner-Plackett method has a significant increase in maximum error. Additional tests were completed with other (larger and smaller) $\epsilon$'s but the largest errors occurred near $\epsilon = 10^{-2}$.

The results suggest that a 24-point Gauss-Legendre integration rule could be used with TVN algorithms based on the equations (13) or (14). These algorithms would produce single precision accuracy for most problems. The equation (13) algorithm appears to be less sensitive to the subtractive cancellation loss of accuracy that occurs with nearly equal $b$ values. All of these algorithms were implemented in Fortran double precision. The average time for a TVN probability computation using a Fortran 77 (double precision) implementation with an 800 MHZ Pentium III processor was $O(10^{-5})$ seconds for the Plackett formula and Gauss-Hermite algorithm implementations; the $\Phi^{-1}$ transform algorithm implementation usually took about 10 times longer.

Some tests were also completed to investigate the possibility of fixed rule double precision TVN algorithms. Some 48-point Gauss rule results are also provided in Table 1. There were significant decreases in the maximum errors, but the results suggest that a reliable double precision algorithm might require a Gauss rule with several hundred points. The use of an adaptive integration was also investigated. A simple globally adaptive algorithm similar to the one used in many of the algorithms for QUADPACK (Piessens, deDoncker, Uberhuber and Kahaner, 1983) was selected. This type of algorithm uses a fixed local integration rule that also provides an error estimate. The rule is first applied to the integrand over the whole integration interval. If the error estimate is larger than the requested accuracy, then the interval, with associated error and integral estimates, is used to initialize a list. The algorithm then proceeds in stages. At each stage an interval with largest error is removed from the list and divided in half. The local integration rule is used on each half of the selected interval and results from the two halves are added to the list. The algorithm terminates when the sum of the local error estimates is less than the requested accuracy. The sum of the local integral estimates is returned as the result. After some experimentation, a 23-point Kronrod rule (see Davis and Rabinowitz, 1984, pp. 106-109) was selected for the local integration rule. This degree thirty-five rule includes an imbedded 11-point Gauss-Legendre rule. The difference between results from the two rules is used to provide a local error estimate. This adaptive algorithm was used with requested accuracy level set at $10^{-14}$. The results are given in the last three rows of Table 1. The algorithms based on equations (10), (13) and (14) all had difficulty reliably achieving the requested accuracy level, with the $\Phi^{-1}$ transform method providing the worst performance. These difficulties are probably due to the sensitivity to subtractive cancellations for all three of the methods. The error averages for the three methods were $1\cdot 10^{-12}$, $3\cdot 10^{-15}$ $3\cdot 10^{-17}$, respectively, for the $\epsilon=0$ test. Overall, the tests suggest that the use of equation (14) with an adaptive algorithm can provide double precision results for most TVN problems with times less than $10^{-4}$ seconds using modern computers. Double precision BVN (BVND) and TVN (part of TVTL) Fortran implementations, with supporting functions, are available from the author's website (www.math.wsu.edu/faculty/genz/homepage, in TVPACK)

Table 1: Maximum Errors for TVN Methods for Grid of $\theta$'s

	$\epsilon$	$\Phi^{-1}$ Transform	Gauss-Hermite	Singular Plackett	Drezner-Plackett
6-Point Rule	0	$1\cdot10^{-3}$	$3\cdot10^{-3}$	$6\cdot10^{-5}$	$1\cdot10^{-5}$
12-Point Rule	0	$1\cdot10^{-4}$	$2\cdot10^{-3}$	$2\cdot10^{-5}$	$2\cdot10^{-7}$
24-Point Rule	0	$2\cdot10^{-5}$	$5\cdot10^{-4}$	$2\cdot10^{-6}$	$8\cdot10^{-9}$
24-Point Rule	.01	$2\cdot10^{-5}$	$3\cdot10^{-4}$	$2\cdot10^{-6}$	$8\cdot10^{-7}$
48-Point Rule	0	$7\cdot10^{-6}$	--	$8\cdot10^{-8}$	$6\cdot10^{-11}$
48-Point Rule	.01	$2\cdot10^{-6}$	--	$9\cdot10^{-8}$	$4\cdot10^{-10}$
Adaptive	0	$4\cdot10^{-9}$	--	$1\cdot10^{-11}$	$3\cdot10^{-14}$
Adaptive	.01	$4\cdot10^{-9}$	--	$4\cdot10^{-13}$	$1\cdot10^{-13}$
Adaptive Times	0	$2\cdot10^{-3}$ s	--	$7\cdot10^{-5}$ s	$8\cdot10^{-5}$ s