Owen's T function

src/SpecialFunctions.jl

@@ -7,6 +7,7 @@ using IrrationalConstants:
     sqrt2π,
     invπ,
     inv2π,
+    inv4π,
     invsqrt2,
     invsqrt2π,
     logtwo,
@@ -58,6 +59,7 @@ export
     logerfcx,
     faddeeva,
     eta,
+    owent,
 
     # Gamma functions
     gamma,
@@ -103,6 +105,7 @@ include("gamma.jl")
 include("gamma_inc.jl")
 include("betanc.jl")
 include("beta_inc.jl")
+include("owent.jl")
 if !isdefined(Base, :get_extension)
     include("../ext/SpecialFunctionsChainRulesCoreExt.jl")
 end

src/owent.jl

@@ -0,0 +1,133 @@
+# Owen's T Function
+#    Written by Andy Gough; August 2021 (see https://github.com/JuliaStats/StatsFuns.jl/issues/99#issuecomment-1124581689)
+#    Edited by Johanni Brea; January 2025
+#    Rev 1.09
+#    MIT License
+#
+# dependencies
+#   IrrationalConstants
+#   SpecialFunctions
+#
+# 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.
+#
+# Partial Derivatives (FYI)
+# D[owent[x,a],x] = -exp(-0.5*x^2)*erf(a*x/sqrt2)/(2*sqrt2π)
+# D[owent[x,a],a] = exp(-0.5*(1+a^2)*(x^2))/((1+a^2)*2π)
+#
+@doc raw"""
+    owent(h, a)
+
+Returns the value of Owen's T function
+```math
+T(h,a) = \frac{1}{2\pi } \int_{0}^{a} \frac{e^{-\frac{1}{2}h^2(1+x^2)}}{1+x^2}dx\quad(-\infty < h,a < +\infty)
+```
+
+For *h* and *a* > 0, *T(h,a)* gives the volume of the uncorrelated bivariate normal distribution with zero mean and unit variance over the area from *y = ax* and *y = 0* and to the right of *x = h*.
+
+## Example
+```
+julia> owent(0.0625, 0.025)
+0.003970281304296922
+```
+
+## References
+"Fast and accurate Calculation of Owen's T-Function", by M. Patefield and D. Tandy, Journal of Statistical Software, 5 (5), 1 - 25 (2000)
+
+"Tables for Computing Bivariate Normal Probabilities", by Donald P. Owen, The Annals of Mathematical Statistics, Vol. 27, No. 4 (Dec 1956), pp. 1075-1090
+#
+"""
+function owent(h::T, a::T) where {T <: Real}
+
+    invsqrt2_T = T(invsqrt2)
+    inv2π_T = T(inv2π)
+
+    #*********************
+    # shortcut evaluations
+    #*********************
+
+    if h < 0
+        return owent(abs(h),a)
+    end
+
+    if h == 0
+        return atan(a)*inv2π_T
+    end
+
+    if a < 0
+        return -owent(h,abs(a))
+    end
+
+    if a == 0
+        return zero(a)
+    end
+
+    if a == 1
+        return T(0.125)*erfc(-h*invsqrt2_T)*erfc(h*invsqrt2_T)
+    end
+
+    if a == Inf
+        return T(0.25)*erfc(sqrt(h^2)*invsqrt2_T)
+    end
+
+    # below reduces the range from -inf < h,a < +inf to h ≥ 0, 0 ≤ a ≤ 1
+    if a > 1
+        return T(0.25)*(erfc(-h*invsqrt2_T) + erfc(-a*h*invsqrt2_T)) - T(0.25)*erfc(-h*invsqrt2_T)*erfc(-a*h*invsqrt2_T) - owent(a*h,one(a)/a)
+    end
+
+    # calculate Owen's T
+
+    if a ≤ T(0.999999)
+        x, w = gauss_legendre(T)
+        return sum(w .* t2.(h, a, x))
+    else
+        # a > 0.999999, T6 from paper (quadrature using QuadGK would also work, but be slower)
+
+        j = T(0.5)*erfc(-h*invsqrt2_T)
+        k = atan((one(a)-a)/(one(a)+a))
+        towen = T(0.5)*j*(one(h)-j)-inv2π_T*k*exp((-T(0.5)*(one(a)-a)*h^2)/k)
+
+        return towen
+    end
+end
+
+t2(h::T, a, x) where T = T(inv4π)*a*exp(-T(0.5)*(h^2)*(one(h)+(a*x)^2))/(one(h)+(a*x)^2)
+
+owent(h::Real, a::Real) = owent(promote(h,a)...)
+
+# 48-point Gauss-Legendre quadrature (Arblib.hypgeom_legendre_p_ui_root!)
+gauss_legendre(::Type{Float64}) =
+    (0.9987710072524261, 0.9935301722663508, 0.9841245837228269, 0.9705915925462473, 0.9529877031604309, 0.9313866907065543, 0.9058791367155696, 0.8765720202742479, 0.8435882616243935, 0.8070662040294426, 0.7671590325157404, 0.7240341309238146, 0.6778723796326639, 0.6288673967765136, 0.5772247260839727, 0.523160974722233, 0.4669029047509584, 0.4086864819907167, 0.34875588629216075, 0.28736248735545555, 0.22476379039468905, 0.1612223560688917, 0.0970046992094627, 0.03238017096286936, -0.03238017096286936, -0.0970046992094627, -0.1612223560688917, -0.22476379039468905, -0.28736248735545555, -0.34875588629216075, -0.4086864819907167, -0.4669029047509584, -0.523160974722233, -0.5772247260839727, -0.6288673967765136, -0.6778723796326639, -0.7240341309238146, -0.7671590325157404, -0.8070662040294426, -0.8435882616243935, -0.8765720202742479, -0.9058791367155696, -0.9313866907065543, -0.9529877031604309, -0.9705915925462473, -0.9841245837228269, -0.9935301722663508, -0.9987710072524261),
+    (0.0031533460523058385, 0.0073275539012762625, 0.01147723457923454, 0.015579315722943849, 0.01961616045735553, 0.02357076083932438, 0.027426509708356948, 0.03116722783279809, 0.03477722256477044, 0.03824135106583071, 0.04154508294346475, 0.04467456085669428, 0.04761665849249048, 0.05035903555385447, 0.05289018948519367, 0.055199503699984165, 0.057277292100403214, 0.059114839698395635, 0.06070443916589388, 0.062039423159892665, 0.06311419228625402, 0.06392423858464819, 0.06446616443595009, 0.06473769681268392, 0.06473769681268392, 0.06446616443595009, 0.06392423858464819, 0.06311419228625402, 0.062039423159892665, 0.06070443916589388, 0.059114839698395635, 0.057277292100403214, 0.055199503699984165, 0.05289018948519367, 0.05035903555385447, 0.04761665849249048, 0.04467456085669428, 0.04154508294346475, 0.03824135106583071, 0.03477722256477044, 0.03116722783279809, 0.027426509708356948, 0.02357076083932438, 0.01961616045735553, 0.015579315722943849, 0.01147723457923454, 0.0073275539012762625, 0.0031533460523058385)
+# 24-point Gauss-Legendre quadrature (Arblib.hypgeom_legendre_p_ui_root!)
+gauss_legendre(::Type{Float32}) =
+    (0.9951872f0, 0.9747286f0, 0.93827456f0, 0.88641554f0, 0.82000196f0, 0.74012417f0, 0.64809364f0, 0.5454215f0, 0.43379351f0, 0.31504267f0, 0.19111887f0, 0.064056896f0, -0.064056896f0, -0.19111887f0, -0.31504267f0, -0.43379351f0, -0.5454215f0, -0.64809364f0, -0.74012417f0, -0.82000196f0, -0.88641554f0, -0.93827456f0, -0.9747286f0, -0.9951872f0),
+    (0.01234123f0, 0.02853139f0, 0.044277437f0, 0.059298586f0, 0.07334648f0, 0.086190164f0, 0.097618654f0, 0.10744427f0, 0.115505666f0, 0.12167047f0, 0.12583746f0, 0.1279382f0, 0.1279382f0, 0.12583746f0, 0.12167047f0, 0.115505666f0, 0.10744427f0, 0.097618654f0, 0.086190164f0, 0.07334648f0, 0.059298586f0, 0.044277437f0, 0.02853139f0, 0.01234123f0)
+gauss_legendre(::Type{Float16}) =
+    (Float16(0.9814), Float16(0.9043), Float16(0.77), Float16(0.5874), Float16(0.368), Float16(0.1252), Float16(-0.1252), Float16(-0.368), Float16(-0.5874), Float16(-0.77), Float16(-0.9043), Float16(-0.9814)),
+    (Float16(0.04718), Float16(0.10693), Float16(0.16), Float16(0.2031), Float16(0.2335), Float16(0.2491), Float16(0.2491), Float16(0.2335), Float16(0.2031), Float16(0.16), Float16(0.10693), Float16(0.04718))
+# 48-point Gauss-Legendre quadrature (Arblib.hypgeom_legendre_p_ui_root!)
+gauss_legendre(::Type{BigFloat}) =
+    (big"0.9987710072524261186005414915631136400889376502767210386129404813754588436074878", big"0.9935301722663507575479287508490741183566147495946719296171518380987546182067713", big"0.9841245837228268577445836000265988305892392234173847299576501679855297780009794", big"0.9705915925462472504614119838006600573024339116308837060283723521653233091284874", big"0.9529877031604308607229606660257183432085413318239187368639476034939458705853333", big"0.9313866907065543331141743801016012677199970856189504298706048642530730422171056", big"0.9058791367155696728220748356710117883122621998274108453524854254710168231209838", big"0.8765720202742478859056935548050967545616485337299619927478757518746727101403824", big"0.8435882616243935307110898445196560498708870117375524015149131998988410546898503", big"0.8070662040294426270825530430245384459730130294604153865758629418121821540044232", big"0.7671590325157403392538554375229690536226423308482073722351285886640508368078524", big"0.7240341309238146546744822334936652465850928122807223627293663025733514606200864", big"0.6778723796326639052118512806759090588499546790260486130710406429754946468798164", big"0.6288673967765136239951649330699946520249089997901617709817329945195319139770715", big"0.5772247260839727038178092385404787728539972861401955280523973994277369963343583", big"0.5231609747222330336782258691375085262891876218118841075802295472194144547473473", big"0.4669029047509584045449288616507985092368121042585169441818691951347943934426029", big"0.4086864819907167299162254958146332864599228429948880647711509833256205384841253", big"0.3487558862921607381598179372704079161343096499683925760321229677812815940686729", big"0.2873624873554555767358864613167976878515583058010397789085000321689998442687597", big"0.224763790394689061224865440174692277438561804041654806164742641045181941897513", big"0.161222356068891718056437390783497694774374379741895117703242637556516342099581", big"0.09700469920946269893005395585362452015273622930093698643058076594480403626262214", big"0.03238017096286936203332224315213444204596280236151809242500322001737781920338223", big"-0.03238017096286936203332224315213444204596280236151809242500322001737781920338223", big"-0.09700469920946269893005395585362452015273622930093698643058076594480403626262214", big"-0.161222356068891718056437390783497694774374379741895117703242637556516342099581", big"-0.224763790394689061224865440174692277438561804041654806164742641045181941897513", big"-0.2873624873554555767358864613167976878515583058010397789085000321689998442687597", big"-0.3487558862921607381598179372704079161343096499683925760321229677812815940686729", big"-0.4086864819907167299162254958146332864599228429948880647711509833256205384841253", big"-0.4669029047509584045449288616507985092368121042585169441818691951347943934426029", big"-0.5231609747222330336782258691375085262891876218118841075802295472194144547473473", big"-0.5772247260839727038178092385404787728539972861401955280523973994277369963343583", big"-0.6288673967765136239951649330699946520249089997901617709817329945195319139770715", big"-0.6778723796326639052118512806759090588499546790260486130710406429754946468798164", big"-0.7240341309238146546744822334936652465850928122807223627293663025733514606200864", big"-0.7671590325157403392538554375229690536226423308482073722351285886640508368078524", big"-0.8070662040294426270825530430245384459730130294604153865758629418121821540044232", big"-0.8435882616243935307110898445196560498708870117375524015149131998988410546898503", big"-0.8765720202742478859056935548050967545616485337299619927478757518746727101403824", big"-0.9058791367155696728220748356710117883122621998274108453524854254710168231209838", big"-0.9313866907065543331141743801016012677199970856189504298706048642530730422171056", big"-0.9529877031604308607229606660257183432085413318239187368639476034939458705853333", big"-0.9705915925462472504614119838006600573024339116308837060283723521653233091284874", big"-0.9841245837228268577445836000265988305892392234173847299576501679855297780009794", big"-0.9935301722663507575479287508490741183566147495946719296171518380987546182067713", big"-0.9987710072524261186005414915631136400889376502767210386129404813754588436074878"),
+    (big"0.003153346052305838632677311543891487578283938831693622295209493250319586438316842", big"0.007327553901276262102383979621786550058707902559201353274881829548806980072502799", big"0.01147723457923453948959266760909162808642050630874764065376681674103503658508731", big"0.01557931572294384872817695583446031397637626899155246951309343105269243335619984", big"0.01961616045735552781446071965221270969581303773413223918112083050740924629812146", big"0.02357076083932437914051930137844923022172973852218859873423906486456506379639118", big"0.0274265097083569482000738362625058204511841551616509759972809374993765019410236", big"0.03116722783279808890206575684635441945428534148356953550954371886143141262424302", big"0.03477722256477043889254858596380241059728139690706809871800663617967672335903626", big"0.03824135106583070631721725652371561786382396835498228892925819103405053922410909", big"0.04154508294346474921405882236106479775347282603403806308273482122272582562965843", big"0.04467456085669428041944858712585039498846278686250200843292144633919149051230188", big"0.04761665849249047482590662347892983015799806674344968539676989627880988507905503", big"0.05035903555385447495780761908786560603299409302590633069379205724693441466024811", big"0.05289018948519366709550505626469891466172648563310918638649123384829276249063296", big"0.05519950369998416286820349519163543900445092560756100054805625793058523675145725", big"0.05727729210040321570515023468470057624152712300411207753884993747681745421856388", big"0.05911483969839563574647481743351991065965560255705499855629113348583514270048057", big"0.06070443916589388005296923202782047788526086425647775511151144466063789427975123", big"0.06203942315989266390419778413759851830638339966509146156903781450273903590161649", big"0.06311419228625402565712602275023331812741364337110079121114724790803811921086588", big"0.06392423858464818662390620182551540891897408498264299989087420749955378258611148", big"0.06446616443595008220650419365770506572569192445553030876055845653739235337295456", big"0.06473769681268392250302493873659155355208191894663651001456309552308307891126462", big"0.06473769681268392250302493873659155355208191894663651001456309552308307891126462", big"0.06446616443595008220650419365770506572569192445553030876055845653739235337295456", big"0.06392423858464818662390620182551540891897408498264299989087420749955378258611148", big"0.06311419228625402565712602275023331812741364337110079121114724790803811921086588", big"0.06203942315989266390419778413759851830638339966509146156903781450273903590161649", big"0.06070443916589388005296923202782047788526086425647775511151144466063789427975123", big"0.05911483969839563574647481743351991065965560255705499855629113348583514270048057", big"0.05727729210040321570515023468470057624152712300411207753884993747681745421856388", big"0.05519950369998416286820349519163543900445092560756100054805625793058523675145725", big"0.05289018948519366709550505626469891466172648563310918638649123384829276249063296", big"0.05035903555385447495780761908786560603299409302590633069379205724693441466024811", big"0.04761665849249047482590662347892983015799806674344968539676989627880988507905503", big"0.04467456085669428041944858712585039498846278686250200843292144633919149051230188", big"0.04154508294346474921405882236106479775347282603403806308273482122272582562965843", big"0.03824135106583070631721725652371561786382396835498228892925819103405053922410909", big"0.03477722256477043889254858596380241059728139690706809871800663617967672335903626", big"0.03116722783279808890206575684635441945428534148356953550954371886143141262424302", big"0.0274265097083569482000738362625058204511841551616509759972809374993765019410236", big"0.02357076083932437914051930137844923022172973852218859873423906486456506379639118", big"0.01961616045735552781446071965221270969581303773413223918112083050740924629812146", big"0.01557931572294384872817695583446031397637626899155246951309343105269243335619984", big"0.01147723457923453948959266760909162808642050630874764065376681674103503658508731", big"0.007327553901276262102383979621786550058707902559201353274881829548806980072502799", big"0.003153346052305838632677311543891487578283938831693622295209493250319586438316842")