1. Introduction
Let
be the standard normal cumulative distribution function (cdf) and
the standard normal probability density function (pdf). The bivariate standard normal cdf
(1)
with a correlation coefficient
has the bivariate standard normal pdf
(2)
If
, then
, while for
we have the degenerated cases
(3)
Additionally, for
define
and
. The functions
and
are symmetric regarding x and y, and because they are complementary, it suffices to analyze only one of them. In this respect, L is often preferred, but not in this paper where some equations for L are also rewritten using
and some other auxiliary equations.
Noting that (3) implies
and
, the equations from ( [1] , p. 937, 26.3.19 and 26.3.20), which do not cover these cases if
, as we get an indefinite expression
for the auxiliary correlation coefficients
and
, can be transformed into
(4)
where
if
or
and
, and
otherwise,
(5)
where
, if
and
, if
. Note that although
is not close to 1,
or
can be, and that
and the lower Equations (5) imply
and
.
Since ( [1] , p. 936, 26.3.9 and 26.3.8) imply
(6)
and
, the problem how to compute
is reduced to the computation of
for
and
.
Extensive older bibliography exists about how
can be computed (see Gupta [2] ). Many old methods, e.g. Donnelly’s method [3] based on Owen’s results [4] , Divgi’s method [5] and the method of Drezner and Wesolowsky [6] , are sufficiently effective and accurate for the majority of present needs, but new methods still appear, e.g. [7] [8] [9] [10] . More details and comparisons of some methods can be found in [11] [12] [13] .
The Owen’s T function is defined by ( [4] , pp. 1078-1079)
(7)
for
. It appears frequently in the extensive Owen’s table of integrals involving Gaussian functions ( [14] , pp. 393-413) and can be used for computing
, as will be explained in the next section. The series (7) is used in Donnelly’s algorithm [3] and, based on it, is also realized in pbinorm() function from the VGAM package [15] . Since it alternates, there are problems with numerical stability, e.g. using double precision arithmetic, it may provide only single precision results. More about Owen’s T function and older ways of calculating it are given by Patefield and Tandy [7] , who developed a hybrid method to compute it. To achieve double precision accuracy, they use the series (7) or one of five other methods. Which method is chosen and which its parameters are used depends on which of the 18 regions
belongs to—see ( [7] , p. 7, Fig. 2). Their method is also implemented in OwenT() function from the OwenQ package [16] , hereafter OwenQ::OwenT().
The method presented in this paper extends the theoretical knowledge and can be easily implemented. Its great advantage over the already mentioned methods, which enable
computation via Owen’s T function, is its simplicity, faster convergence and ensured numerical stability. Regarding other methods that have proven to be good or even excellent in terms of
computation speed and accuracy, e.g. the one implemented in the pmvnorm() function from the mvtnorm package [17] , the main advantages are simplicity and independence from various constants. Both enable easy portability between different environments and using high-precision arithmetic. There is no need to adjust the number and precision of various constants, such as weights and knots in quadrature formulas, which are used in many
computational methods.
In Section 2, the background results are presented, which are the starting point for the development of a new method, as well as for the acceleration of the tetrachoric series, which was used for the benchmark values computation. In Section 3, auxiliary results are derived. In core Section 4, new series for
are derived in Theorems 3 and 4, and in Corollary 2 they are combined in a unified series for
. However, in all cases, the series must be supplemented by the standard normal cdf and some of them also by the arctangent function computation. Four series for computing the Owen’s T function are presented in Theorem 5. Among them (40b) can be viewed as a modified Euler’s arctangent series. In Sections 5 and 6, algorithmic and numeric issues are discussed, and testing results are presented.
Using double precision arithmetic and testing on a random sample of one million triplets
with uniformly distributed components, the maximum absolute error was 3.45 × 10−16, which is approximately 3.11-times the machine epsilon 2−53. On the same random sample, but with
transformed so that more than half of
were greater than 0.9999, and by equivalent computing of two function values instead of
in rare cases with
, the maximum absolute error was even slightly smaller.
All new series are suitable for vector and high-precision computation. The
fastest one asymptotically converges as
, where
. Using the Rmpfr package [18] and 1024-bit precision computation, 199 iterations were needed to compute
, achieving an absolute error smaller than
.
2. Background Results
The bivariate standard normal cdf can be expanded into tetrachoric series ( [19] , p. 196), ( [1] , p. 940, 26.3.29)
(8)
where
are the “probabilist’s” Hermite polynomials defined by
The tetrachoric series (8) converges a bit faster than a geometric series with a quotient
and slowly converges for
close to or equal to 1.
Let
. Extending the usual double factorial function definition for
to
with
, we have
Setting
in (8), transforming the summation index
, and using
for
, we get a slightly faster converging series
(9)
Vasicek ( [20] , p. 7) found an alternative series for the computation of
, which is a good alternative to the series (9) if
.
Integrating the well known equation ( [21] , p. 352)
(10)
with respect to
, using equations
and (3), we get
(11)
Due to historical reasons, the variates h and k are often used instead of x and y. We will use h instead of x but will not need k in the role of y.
In the sequel, let
and
be permanently connected by the
one-to-one correspondence
and
. Notice that (4) and
imply
. Setting
,
, and writing
instead of
, the Equations (11) can be rewritten to
(12)
where, using a substitution
,
Analogously,
and
are defined for
by integrals with limits
and
, and
and
respectively, which do not exist for
. The computation of
is reduced to the computation
of one of the integrals
,
and
.
The special cases (12) of the Equations (11) in combination with (4) are a starting point for some computation methods for
which are based on
the integration. Owen’s is one of them—note that
. Using r instead of a in the sequel and resolving the problem with
as before, Owen’s equation ( [14] , p. 416, 3.1) can be rewritten as
(13)
where
if
or
and
, and
otherwise, and
(14)
Another Owen’s important equation is ( [14] , p. 414, 2.7)
(15)
where
if
and
if
. It implies
. If
, then
. Using (12) and the right Equation (6), we get
(16)
Let
,
,
,
and
. By adding
on both sides of (15) and using (7) and (12), it follows
(17)
where
if
and
if
.
If
, then
,
and
.
can be computed by computing
and using (17). Consequently, the right Equation (6) and the transformation of the parameters
and
, if needed, reduce the condition
to
with a corresponding
. The problem how to compute
is reduced to the computation of
for
and
.
The theorems in the sequel do not depend on the restriction to
with a corresponding
found in [7] . However, if
, the transformation of
and
speeds up convergence significantly. This is, for example,
also the case for the Vasicek’s series and the tetrachoric series (9), which transforms into
(18)
where
,
if
and
if
. Since the
results calculated with the Vasicek’s series in Section 6 only supplement benchmark values calculated with the tetrachoric series, details about it are not given here.
3. Auxiliary Results
This paper pays attention to the estimation of the method error due to truncation of infinite series, while the analysis of rounding errors is out of its scope. The following lemma is intended for this. In it and all the theorems below, the phrase “truncated after the nth term” also means omitting the entire series. In these cases, set
and assume the empty sum is 0.
Lemma 1. If the Euler’s series
(19)
is truncated after the nth term, then the remainder
has the same sign as r and is bounded in absolute value by
(20)
Proof. Let
. Noting that
and
for
( [22] , p. 397, 3.621, 4.), and using Tonelli’s theorem, it follows
(21)
The denominator in the last integrand is an increasing function of x. Inserting 1 and 0 instead of x implies
(22)
and, together with
, (20). □
If
,
, (19) with
instead of r, and
imply
and (20) with
. This series for
converges faster than (19), so we can already sense some kind of similarity with (15).
In the sequel, the lower and upper regularized gamma functions
and
respectively, where
and
are the lower and upper incomplete gamma functions respectively, both defined for
and
, will be very important, as well as their lower and upper bounds.
In the Alzer’s inequality
(23)
where
( [23] , p. 772, Theorem 1), ( [24] , p. 221, (5.4)), allowing
resulted in strict inequalities becoming non-strict. It implies
(24)
For
and
, we have
( [22] , p. 908, 8.352, 1.), where
is the remainder of the Taylor series of the exponential function. The Lagrangian form of the remainder
, where
, implies
for
. Since
for
( [25] , pp. 323-324, 3.8.25), it follows
(25)
where the equality holds only for
.
Theorem 2. Let
and
be absolutely convergent series for
. Then
(26)
and
(27)
for
and
.
Proof. Let
and
for
. Using
and a substitution
, we get
Setting
implies
for
. Then
, where the first sum is finite. Using the fact that
is a decreasing function of a ( [26] , p. 154), (25) and the absolute convergence assumption of
, it follows
and
Noting that
, we proved that
. By Fubini’s theorem it follows that
, hence (26).
Let
and
for
. Analogously to above, and using
, we get
,
and (27). □
Note that in the same way as the Equations (26) and (27) were derived, analogous equations can be derived if the integration limits 0 and x are replaced by x and
respectively, and lower regularized gamma function is replaced by the upper one. However, for each case it should be checked whether the order of integration and summation can be reversed.
The usefulness of Theorem 2 is that if the sequences
and
can be computed recursively, then the terms of the series (26) and (27) can be too. Let
. Crucial is the calculation of
and
,
, the rest is trivial. Since
, we will calculate
instead of
.
Using
( [1] , p. 934, 26.2.30), we get
, and
. Since
( [1] , p. 262, 6.5.21) implies
, using auxiliary variables
and
, it follows
for
and
, and
(28)
for
4. Main Results
The core of this article are Theorems 3, 4 and 5. Since the first two do not specify variants for the transformed parameters, the last one, together with (13) and (14), is more practically useful. It also justifies the title of the article.
Theorem 3. Let
,
,
and
. Then
(29)
where
. The series
converges for every
. If it is truncated after the nth term, then
(30)
where
and
is bounded by
(31)
where
.
Proof. From (7) and
for
( [1] , p. 932, 26.2.11), it follows
(32)
, where
and
. Since
the series on the left side converges for every
and is an odd function of
. If
, then all its terms are positive, implying that it converges absolutely.
All requirements of Theorem 2 are fulfilled. Assuming
and using
, (32) and (27), it follows
which, together with (12), (7),
,
and
, imply (29) for
. Since (7) implies
, (29) is also valid for
. Let
. Using (23), (21) and (22), it follows
and (31). □
Note that if
, then
and
. Hence the sequence
is decreasing.
Corollary 1. For every
we have
(33)
Proof. Using (29) with
, hence
,
and the right Equation (16), it follows
and a quadratic equation
. If
, then
and (33) is valid. If
, since
if
and
if
, only
is a proper root of the two. Note that sending
implies
. □
Theorem 4. Let
,
,
and
. Then
(34)
where
. The series
converges for every
. If it is truncated after the nth term, then
(35)
where
and
is bounded by
(36)
Proof. Using (29), (19),
,
and
, it follows (34). Noting that
and (22) imply (36). □
Note that the upper bound (24) is not suitable for a simple estimation of
here because
, hence
was implicitly used instead.
is an increasing function of a. The ratio
can be greater than 1, e.g. if
and
, we get
. In Theorem 4, the sequence
is not always decreasing, but, in any case, the well-behaved
series (19) divided by 2π is a majorant series for
from Theorems 3 and 4.
Corollary 2. Let
,
,
,
,
and
. Let
and
for
. Then
(37)
where
is S with
instead of
, and
if
or
and
, and
otherwise.
Proof. Since
and
,
and
, (4), (5), (29) and (34) imply (37). □
Note that
and
, therefore they are invariant on the transformation of
and
, but
,
and
from Theorems 3 and 4 are not because
and
if
. They transform to
,
and
. By transforming the parameters and using (17), the Equations (30) and (35) can be rewritten to
(38)
and
(39)
respectively, where
if
and
if
.
The transformation makes sense if
, implying
and faster con
vergence. For example, the inequality
from Theorem 4 transforms into
, where the upper bound
is
-times smaller than a corresponding upper bound
. The noteworthy cost for faster convergence is the additional computation of
.
Theorem 5. If
, the Owen’s T function (7) can be written as
(40a)
(40b)
(40c)
(40d)
where
,
, and
if
and
if
.
Proof. Using (12) and (7), the Equations (29) and (34) imply (40a) and (40b) respectively, which, together with (17), imply (40c) and (40d) respectively. □
Note that when truncating the series (40a) and (40b), the truncation error in absolute value can be estimated by (31) and (36) respectively. The same upper
bounds with
instead of r also apply to (40c) and (40d).
Because of (40b) and
, the Owen’s T function can be seen as a modified Euler’s arctangent series (19) divided by 2π, which is its majorant series.
5. Recursion and Asymptotic Speed of Convergence
Let
and
, where S and
relate to Theorem 4. Using (34) and (28), the series
can be calculated by Recursion 1 (tildes are omitted). Analogously, let
and
, where S and
relate to Theorem
3. The series
can be calculated by Recursion 2 (hats are omitted). In practice, both recursions can be easily combined, since only four steps, marked with
, are different.
In both cases, given
, the recursion is terminated when the condition
is met and
is taken as its result. However, if
is too small, or even intentionally negative, the calculation proceeds up to the actual accuracy threshold—the recursion is terminated when the condition
is met. This criterion was used in the calculations presented in the next section.
In both recursions, all variables are positive and are bounded upward. If
, the sequence
strictly increases and, due to (24), converges to 1. The sequence
is strictly increasing if
. If
, the recursion terminates
with
. Note that (7) implies
. If step 7 in Recursion 2 is replaced by
, the sequence
changes its sign but keeps the direction. However, the modified Recursion 2 computes
.
Let
and
be the results of Recursion 1 and modified Recursion 2 respectively, calculated with the same h and r if
, or with
and
if
. Let
be calculated by (35) or (39), where unknown remainders are ignored. Analogously, let
be calculated by (30) or (38), where also
is ignored. If
, then
if
and
if
. Even if
, one of
and
underestimates and the other overestimates
. Consequently, the upper bounds (31) and (36) are not needed if
and modified
are calculated in the same recursion which terminates when
become small enough. Upper bounds are also not needed if in Recursions 1 and 2 step 17 is replaced so that
is the only termination condition. In this case, steps 3, 8 and 15 are unnecessary.
Given q, which one of the series (29) and (34) converges faster? Since
is a cdf of the gamma distribution with the median
, it follows
if
. The an
swer depends on the number of considered terms, which depends on the prescribed accuracy. On this basis, it seems that series (29) is more suitable than (34) for “small” q and vice versa for “large” q. Table 1 and Table 4 in the next section show how this is reflected in practice.
Gautschi’s inequality
( [27] , p. 138, 5.6.4) with
implies
. Combining it with (25), sending
and neglecting non-essential constants imply that the series (29) asymptotically converges as
. Assuming calculation with
,
and (38) instead of h, r and (30) respectively if
, hence
, and recalling that q is invariant on the transformation of
and
, we get an asymptotic convergence as
. The series (29) or its corresponding transformed version, unwritten and based on (40c), asymptotically converges as
, where
. On the other hand, since (24) implies
, the series (34) or its corresponding transformed version, unwritten and based on (40d), asymptotically converges as
, hence in the worst case as
.
Based on both asymptotic convergences, comparing (34) to (29) would be unfair because the calculation of
should also be taken into account for (29). Since it is implicitly embedded in Recursion 1, it is calculated on the fly at the cost of worse speed of convergence. Each of the recursions is executable with the basic four arithmetic operations, two calls of elementary functions (if the calculation of r from
is also considered), and
in the modified Recursion 2. However, the computational cost of
(and
if needed) should also be taken in account, but it is the same in (29) and (34).
6. Numeric Issues and Testing Results
A common cause of numerical instability, where two absolutely large values give an absolutely small result when added or subtracted, cannot occur in Recursions 1 and 2. The calculation by them is stable, but something may happen that needs
to be pointed out. The side effect of transforming
and
is that
can become unusually large. When computing
, according to Marsaglia ( [28] , p. 2), “[...] the region of interest may not be the statistician’s
, say, but rather
[...]’’. Assume the same needs when computing
,
and using transformed parameters if
. If
and
, then
,
and
, implying that the
region of interest could be
. Even the usual
transformed to
is well beyond the established limits.
The “problem” of too large
is mentioned because we have to be aware of it if we are going to perform high-precision computation. In any case, the recursion must be complemented by a suitable computation of
and
if needed. As an example that not all are like that, the computation by the well
known series
( [1] , p. 932, 26.2.10) and the
standard double precision arithmetic is numerically unstable even for moderately large
, e.g. using the programming language R [29] , for
we get a completely wrong result, approximately −1683. A similar problem can lead to a significant decrease in accuracy if
is calculated with (7), as will be seen in Subsection 6.1.
In double precision arithmetic, the minimal positive non-zero number is
. Assuming
, it follows
and
. If
, then
is between
and
. If
, then all elements of the sequence
are zero. Consequently, all elements of
are one and all elements of
in Recursion 2 are zero.
From the integral in (7), it follows
if
, where the upper bound is a decreasing function of
. If
and
, then
and
. Analogously to the case when
, the recursion terminates with
, which is practically a correct result for
, as well as
if
, and
if
.
The Recursions 1 and 2 with adjustment of their results according to (40a)—(40d) are implemented in the package Phi2rho [30] . It contains the functions OwenT(), hereafter Phi2rho::OwenT(), and Phi2xy(), which compute
and
respectively. They also enable computation with the tetrachoric and Vasicek’s series, original or accelerated. The latter means that h and r are transformed if
and
for the tetrachoric and Vasicek’s series respectively, and the Equation (15) is used, so (18) instead of (9) in the case of tetrachoric series. For the sake of comparison, all series were considered during testing in Subsection 6.1. In Subsection 6.2, the original tetrachoric and Vasicek’s series were excluded because they converge too slowly for unfavorable
.
During testing, the values of
and
were computed for a selected and randomly generated set of parameters. Reference values were computed using the accelerated tetrachoric series and quadruple (128-bit) precision. Absolute differences between tested and reference values are declared as absolute errors. The maximum absolute error is considered a measure of accuracy. Since the mean values and some quantiles are also important, they are shown in the tables in Subsection 6.2. An integral part of the testing was also a comparison of the results obtained with the Phi2rho package and competing packages available on CRAN.
6.1. Numeric Issues and Results When Computing
Numerical calculations of
were performed on 39,999 grid points
,
where h ranges from −10 to 10 in increments of 0.1,
and
ranges
from −0.99 to 0.99 in increments of 0.01. Since for
we get
, the actual range of h for which
must be computed is approximately from −70 to 70, so well beyond the
.
Note atanExt = no in tables means that Recursion 1 was executed, hence (40b) was used if
and (40d) if
, and atanExt = yes means that Recursion 2 was executed, hence (40a) was used if
and (40c) if
. Both recursions are well suited for vectors, so each of the parameters h and r can be a scalar or a vector. If they are both vectors, they must be of the same length, if one of h and r is a vector and the other is a scalar, the latter is replicated to the length of the former. Since the number of iterations is determined by the “worst’’ component, we would expect execution to be faster if function calls are done in a loop for each component of the vector separately, but usually the run time is significantly increased due to the overhead of the loop, more initialization, etc.
For the consistency test,
was calculated on the test grid with h and r both scalars, one a scalar and the other a vector. For all three calculations with Recursion 1, the results matched perfectly, as well as for all three with Recursion 2. The maximum absolute difference between calculations with Recursions 1 and 2 was approximately 9.49 × 10−17.
When testing accuracy, the tetrachoric and Vasicek’s series, original and accelerated, were also tested, using
. By transforming the
parameters, if appropriate, both comparison methods gained a lot, as can be concluded from Table 1. The absolute errors, the maximum values of which are
Table 1. Double precision computation of
on the test grid.
collected in this table, were calculated against the reference values computed with the accelerated tetrachoric series and 128-bit precision.
Since the vector computations were performed, there are two average numbers of iterations. The first is empirical because the iterations were terminated according to the “worst’’ component of the vector, and the second is hypothetical, if the computations were performed for each component separately. Since all iterations ran to the accuracy limit and with an unbounded number of iterations, the original Vasicek’s series needed the enormous number of iterations, not only for the worse case, but also in average, when calculating with scalar h and vector
r. However, it is not intended for cases with
, and its error can be well estimated, which was ignored here.
Comparative computations with competing packages were also performed on
the test grid, using
. The pbinorm() function from
the VGAM package [15] , i.e. using (7), had an average and maximum absolute error of 1.12 × 10−8 and 3.62 × 10−7 respectively. Similar maximum absolute errors were obtained by the function pbvnorm() from the pbv package [31] , which is based on [6] , and pmvnorm() from the mvtnorm package, but only with the parameter algorithm = Miwa.
6.2. Numeric Issues and Results When Computing
In all
computational methods based on (4) and (5), or (13) and (14), exceptional cases with
,
and
are not problematic because (3) can be used, however, nearby cases are their weak point. Recalling
that
and using (10) and (2), it follows
and
, except if
. If
is close to 1 and
, then
and q may be small enough that
is moderate
or even large. As a result, even a small rounding error in the calculation of key variables that depend on
can cause a much larger error in the result. The calculation of
and
with (14) is particularly problematic. In a similar situation, Meyer ( [9] , p. 5) recommends an equivalent but less sensitive calculation
by
if
and
if
, and
analogously for
. In our case, his proposal improves accuracy, but does not solve the problem sufficiently.
The problem could be solved by a hybrid computation along the lines of some other methods, e.g. the series from ( [32] , pp. 242-245) could be used if
and
. Using the initial idea to develop these series, i.e. axis rotation in a way that removes the cross product term
in (1), by setting
and
the Equation (1) transforms to
and, using
and
, to
By transforming
in the second integral and using ( [14] , p. 402, 10,010.1),
it follows
, where
. This equation is suitable for
. From it and
we obtain the equation
where
, that is suitable for
. Both are combined in
(41)
where
,
and
. Noting that
hence
, we stop here because one potentially critical case is replaced with two non-critical ones.
If we continued to simplify the right side of (41) by (13), we would get (13). This actually happens when calculating with the functions used in tests, but it turns out that the results are more accurate than if we start with (13) if
and
. Of course, the Equation (41) could be used as a starting point instead of (13) in all cases, but this is not rational because it requires 4 calculations of the T function instead of 2, two of which cancel each other out, at least in theory, and as a result, in non-critical cases the error may be larger than when starting with (13). When testing, (41) was used as a starting point instead
of (13) if
. However, since
, using (41) may
also make sense in some other circumstances, e.g. when computing directly with (8).
In this case, testing was performed on the randomly generated set of parameters. After setting the random generator seed to 123 to allow replication and independent verification, one million x values, one million y values, and one million ρ values were drawn from a uniform distribution on the intervals
,
, and
respectively.
For each triplet
, auxiliary values
and
were calculated. For the new series, the tested
values were computed using (41) as a branch point and (13) for both branches if
, and using (13) directly otherwise. The same applies for the tetrachoric and Vasicek’s series, only (4) must be used instead of (13). However, the described procedure is built in the Phi2xy() function. The tested
values needed for the preparation of Table 2 were computed with:
row 1: Phi2xy(x, y, ρ, fun = "tOwenT");
row 2: Phi2xy(x, y, ρ, fun = "vOwenT");
row 3: Phi2xy(x, y, ρ, fun = "mOwenT", opt = FALSE);
row 4: Phi2xy(x, y, ρ, fun = "mOwenT", opt = TRUE).
The absolute errors, which statistics are collected in Table 2, were calculated against reference values computed with the accelerated tetrachoric series as described, but with 128-bit precision. From a user’s point of view, only the parameters x, y and ρ needed to be converted into 128-bit “mpfr-numbers’’ before calling Phi2xy(x, y, ρ, fun = "tOwenT"), hence x <- mpfr(x, precBits = 128) and analogously for y and ρ. In this case, the computation took significantly more time, since all intermediate and final results were 128-bit “mpfr-numbers’’.
In this test,
,
, there are 5 cases with
and
.
In order to better check the accuracy if
is close to 1, the computation was repeated with the same x and y, and
. In this case, more than half of
are greater than 0.9999 and none is equal to 1. The results are in Table 3.
In this test, there are 178 cases with
and
. For both novel series in Table 3, the maximum absolute error is 1.54 × 10−14 if (41) is not used, improving to 1.21 × 10−15 if Meyer’s advice is followed. Similarly applies to the other two series.
Table 2.
computation—absolute error comparison.
Table 3.
computation—absolute error comparison.
Due to the time-consuming 128-bit precision computation, the presented tests are the only performed tests based on the 128-bit precision reference values. However, additional tests were made by comparing the results of Phi2xy() with the results obtained by other functions from packages available on CRAN. For three functions, the accuracy was found to be worse than 2 × 10−7 as measured by the maximum absolute error on the test grid from Subsection 6.1. Regarding accuracy, only the function pmvnorm() with the parameter algorithm = TVPACK from the mvtnorm package proved to be equivalent to Phi2xy() and even slightly more accurate with the maximum absolute errors 2.58 × 10−16 and 2.19 × 10−16 on the test sets of triplets
and
respectively. The same function with algorithm = GenzBretz achieved 2.58 × 10−16 and 2.20 × 10−11.
The OwenQ::OwenT() function also proved to be equivalent to Phi2xy() as measured by the maximum absolute error on the test sets of triplets, but only when upgraded to compute
in the manner used in Phi2xy(). This function and Phi2rho::OwenT() are essentially wrappers for the internal functions tOwenT(), vOwenT() and mOwenT(), which as workhorses compute series.
In terms of reliability, stability and accuracy of the new methods, no problems or significant differences according to the parameters from the different data sets were detected, except those already described and related to
and
. To eliminate them, the Equation (44) was derived, which also proved to be successful when using the OwenQ::OwenT() function, upgraded to compute
. Regarding parameters from different data sets, the only detected big difference is in the number of iterations, as can be seen from Figure 1.
Two million
computations on Windows 10 and 64-bit R 4.3.1 on the old Intel i7-6500U CPU @ 2.50 GHz and 8 GB RAM, using only one thread
Figure 1.
computation with the new series (atanExt = yes) with double (left) and 128-bit (right) precision—number of iterations.
of four, lasted 26, 31 and 286 seconds for Phi2xy(), upgraded OwenQ::OwenT() and pmvnorm() with algorithm = TVPACK respectively. The 128-bit precision reference values computation took over 17 hours, compared to 43 seconds for the double precision one. Based on Table 1 and Table 4, we can conclude that the increased number of iterations has only an insignificant effect on the huge time extension factor.
6.3. High-Precision Computation
All functions in the Phi2rho package are ready to use the Rmpfr package, which enables using arbitrary precision numbers instead of double precision ones and provides all the high-precision functions needed. It interfaces R to the widely used C Library MPFR [33] . Assuming that the Rmpfr package is loaded, the functions should only be called with the parameters, which are “mpfr-numbers’’ of the same precision. All 128-bit precision benchmark values used in Subsections 6.1 and 6.2 were calculated using Rmpfr.
To get a sense of how the number of iterations depends on the required precision, the test from Subsection 6.1 was partially repeated using the 128-bit precision computation. The results are collected in Table 4. Note that
and that in this case the values in the table are not the maximum absolute errors because the comparison values are computed with the same precision.
From Table 1 and Table 4 can be concluded that there is no significant difference in accuracy between the compared series if the accelerated series are considered, but one of the new series converges significantly faster than the others. It can also be concluded that the Vasicek’s series with
instead of
is a
reasonable alternative to the tetrachoric series (8) also for
and remains a reasonable alternative for
if compared to the transformed tetrachoric series (9).
Table 4. 128-bit precision computation of
on the test grid.
Table 5. A special case of high-precision computing—absolute error comparison.
The number of iterations for the computation with the new series with double and 128-bit precision can be compared in Figure 1. Only the first quadrant is presented because others are its mirror images.
We also provide examples with
and
for which reference values can be computed alternatively.
and
were computed by (13) and (14), and by the right side of (16) by the pnorm() function from the package Rmpfr as a benchmark. In all cases, the latter was computed by double the precision expressed in bits and used for the former. Both reference values were used to compare how quickly the compared series converge when computing with more than 128-bit precision. In Table 5 and Table 6, for each of them the absolute error and the number of iterations, which is the same for both computations, is presented. Since
, there is no difference between the original and accelerated tetrachoric series, and the same applies to Vasicek’s series.
Time-consuming high precision testing was less intensive than double precision one. Since there is no competing package on CRAN ready for high-precision computation of
, it could not be supplemented by comparison computations with competing packages.
Table 6. A special case of high-precision computing—absolute error comparison.
7. Conclusions
If the generally applicable Theorem 2 is not considered, Theorems 3 and 4, and Corollary 2 are interesting mainly theoretically, because Theorem 5, together with (13) and (14), replaces and supplements them in practice. The series (40a)—(40d) enable fast and numerically stable computation, which is often more important than the speed of convergence. From Table 1 and Table 4 can be concluded that the computation with (40a) or (40c) needs significantly fewer iterations than the computation with the accelerated tetrachoric and Vasicek’s series, and from Table 5 and Table 6 can be assumed that the advantage increases with increasing precision. For 53-bit precision (double precision), it required a half, and for 1024-bit precision only a fifth, of the iterations demanded by the competing series. However, the cost of calculating
is not taken into account. In connection with this, let us just mention two more questions.
Recalling that
, the Taylor series of the arctangent function and (7) imply
which Fayed and Atiya ( [32] , p. 241) have already noticed. If the above series, (29) and (34) are viewed as functions of the variable h, they have a similar structure, but the coefficients, depending on r, belong to different arctangent series. Whether the similarity could be explained by a similar Euler transform as the one which transforms the Taylor arctangent series to the Euler’s arctangent series seems interesting question, but was not deeply explored.
In Recursion 1, the external calculation of the arctangent function is avoided, assuming that the Euler’s series (19) is used for it, and the expectation that there would be no difference in the number of iterations for Recursions 1 and 2 if those for calculating the arctangent were also taken into account. However, the use of the external arctangent calculation allows the use of already known and possible future faster converging series. Such is the case with
(42)
which is based on the Taylor series of
and
.
Indeed, the number of iterations is not significantly smaller if
and the accuracy obtained by double precision computation is sufficient, and even the square root must be calculated, but it could be significantly smaller in high-precision computation. The Euler’s series (19) is included in [1] and [27] , but (42) is not, even though it converges faster. It can be found in ( [22] , p. 61, 1.644 1.) as
with a reference to the 1922 source. Due to equations
and
, the coefficients are the same as those in (42). Whether we could also find a corresponding series for the Owen’s T function, which would converge faster than (40b), remains an open challenge.