HISTORY
In the past 20 or so years, most implementations of Owen's T function have followed the algorithms given in "Fast and accurate Calculation of Owen's
T-Function", by M. Patefield and D. Tandy, Journal of Statistical Software, 5 (5), 1 - 25 (2000)
Six algorithms were given, and which is was used depends on the values of (h,a)
T1: first m terms of series expansion of Owen (1956)
T2: approximates 1/(1+x^2) by power series expansion up to order 2m
T3: approximates 1/(1+x^2) by chebyshev polynomials of degree 2m in x
T4: new expression for zi from T2
T5: Gauss 2m-point quadrature; 30 figures accuracy with m=48 (p. 18)
T6: For when a is very close to 1, use formula derived from T(h,1) = 1/2 Φ(h)[1-Φ(h)]
They developed code for these algorithms on a DEC VAX 750. The VAX 750 came out in 1980 and had a processor clock speed of 3.125 MHz.
The reason for 6 algorithms was to speed up the function when possible, with T1 being faster than T2, T2 faster than T3, etc.
THIS FUNCTION
A native Julia implementation, based on the equations in the paper. The FORTRAN source code was not analyzed, translated, or used. This is a new
implementation that takes advantages of Julia's unique capabilities (and those of modern computers).
T1 through T4 are not implemented. Instead, if a < 0.999999, T5 is used to calculate Owen's T (using 48 point Gauss-Legendre quadrature)
For 0.999999 < a < 1.0, T6 is implemented.
REFERENCES
[1] "Fast and accurate Calculation of Owen's T-Function", by M. Patefield and D. Tandy, Journal of Statistical Software, 5 (5), 1 - 25 (2000)
[2] "Tables for Computing Bivariate Normal Probabilities", by Donald P. Owen, The Annals of Mathematical Statistics, Vol. 27, No. 4 (Dec 1956), pp. 1075-1090
blackeneth
commented
3 years ago
owent(h,a) returns Owen's T function for (h,a)
HISTORY In the past 20 or so years, most implementations of Owen's T function have followed the algorithms given in "Fast and accurate Calculation of Owen's T-Function", by M. Patefield and D. Tandy, Journal of Statistical Software, 5 (5), 1 - 25 (2000)
Six algorithms were given, and which is was used depends on the values of (h,a)
T1: first m terms of series expansion of Owen (1956) T2: approximates 1/(1+x^2) by power series expansion up to order 2m T3: approximates 1/(1+x^2) by chebyshev polynomials of degree 2m in x T4: new expression for zi from T2 T5: Gauss 2m-point quadrature; 30 figures accuracy with m=48 (p. 18) T6: For when a is very close to 1, use formula derived from T(h,1) = 1/2 Φ(h)[1-Φ(h)]
They developed code for these algorithms on a DEC VAX 750. The VAX 750 came out in 1980 and had a processor clock speed of 3.125 MHz.
The reason for 6 algorithms was to speed up the function when possible, with T1 being faster than T2, T2 faster than T3, etc.
THIS FUNCTION A native Julia implementation, based on the equations in the paper. The FORTRAN source code was not analyzed, translated, or used. This is a new implementation that takes advantages of Julia's unique capabilities (and those of modern computers).
T1 through T4 are not implemented. Instead, if a < 0.999999, T5 is used to calculate Owen's T (using 48 point Gauss-Legendre quadrature) For 0.999999 < a < 1.0, T6 is implemented.
REFERENCES [1] "Fast and accurate Calculation of Owen's T-Function", by M. Patefield and D. Tandy, Journal of Statistical Software, 5 (5), 1 - 25 (2000) [2] "Tables for Computing Bivariate Normal Probabilities", by Donald P. Owen, The Annals of Mathematical Statistics, Vol. 27, No. 4 (Dec 1956), pp. 1075-1090
Partial Derivatives (FYI) D[owent[x,a],x] = -exp(-0.5x^2)erf(ax/sqrt2)/(2sqrt2π) D[owent[x,a],a] = exp(-0.5(1+a^2)(x^2))/((1+a^2)*2π)