Abstract:The Owen's T function is presented in four new ways, one of them as a series similar to the Euler's arctangent series divided by $2\pi$, which is its majorant series. All possibilities enable numerically stable and fast convergent computation of the bivariate normal integral with simple recursion. When tested $\Phi_\varrho^2(x,y)$ computation on a random sample of one million parameter triplets with uniformly distributed components and using double precision arithmetic, the maximum absolute error was $3.45\cdot 10^{-16}$. In additional testing, focusing on cases with correlation coefficients close to one in absolute value, when the computation may be very sensitive to small rounding errors, the accuracy was retained. In rare potentially critical cases, a simple adjustment to the computation procedure was performed - one potentially critical computation was replaced with two equivalent non-critical ones. All new series are suitable for vector and high-precision computation, assuming they are supplemented with appropriate efficient and accurate computation of the arctangent and standard normal cumulative distribution functions. They are implemented by the R package Phi2rho, available on CRAN. Its functions allow vector arguments and are ready to work with the Rmpfr package, which enables the use of arbitrary precision instead of double precision numbers. A special test with up to 1024-bit precision computation is also presented.
What problem does this paper attempt to address?
The problem that this paper attempts to solve is to provide a new method for calculating the bivariate normal integral, especially by using Owen's T - function to represent it in a form similar to the Euler arctangent series. Specifically, the author proposes four new ways to express Owen's T - function, one of which is to represent it as the form of the Euler arctangent series divided by \(2\pi\), which helps to achieve numerically stable and rapidly convergent calculations.
### Core of the problem
1. **Calculation of the bivariate normal integral**: The bivariate normal integral is very important in statistics and probability theory, especially when dealing with correlated variables. Although traditional methods are effective, in some cases (such as when the correlation coefficient is close to 1), problems such as numerical instability or slow convergence may occur.
2. **Numerical stability**: When the correlation coefficient \(\rho\) is close to 1, the calculation may be very sensitive to small rounding errors, resulting in inaccurate results. Therefore, it is crucial to find a calculation method that can maintain high precision and is numerically stable.
3. **Rapid convergence**: In order to improve the computational efficiency, a series expansion form that can converge rapidly is required. The new method proposed by the author based on the Euler arctangent series aims to achieve this.
### Main contributions
- **New series representation**: The author proposes four new series representation methods, one of which is to represent Owen's T - function in the form of the Euler arctangent series, which not only simplifies the calculation but also improves the numerical stability.
- **Numerical stability**: By using these new series representations, the author ensures that even when the correlation coefficient is close to 1, the calculation can still maintain high precision.
- **Rapid convergence**: The new method shows a faster convergence speed in most cases, thus improving the computational efficiency.
- **High - precision calculation**: The new method is suitable for vector and high - precision calculations, and can be implemented through the R package `Phi2rho`, supporting arbitrary - precision calculations.
### Test results
The author conducted tests on 1 million randomly generated parameter triples, using double - precision arithmetic, and the maximum absolute error was only \(3.45\times 10^{-16}\), which is approximately 3.11 times the machine precision \(2^{-53}\). In addition, in the case where the correlation coefficient is close to 1, high precision can still be maintained by adjusting the calculation process.
### Formula summary
Owen's T - function is defined as:
\[ T(h, a)=\frac{1}{2\pi}\int_{0}^{a}\frac{e^{-\frac{1}{2}h^{2}(1 + x^{2})}}{1 + x^{2}}\,dx \]
One of the new series representations is:
\[ T(h, r)=\frac{\arctan(r)}{2\pi}-\frac{r}{2\pi(1 + r^{2})}\sum_{k = 0}^{\infty}\frac{(2k)!!}{(2k+1)!!}P(k + 1, q)\left(\frac{r^{2}}{1 + r^{2}}\right)^{k} \]
where \(q=\frac{1}{2}(1 + r^{2})h^{2}\), and \(P(k + 1, q)\) is the normalized form of the lower incomplete gamma function.
These improvements make the calculation of the bivariate normal integral more efficient, stable and accurate.