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,142 @@
+# 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 to make type stable; January 2025
+#    Rev 1.09
+#    MIT License
+#
+# dependencies
+#   IrrationalConstants
+#   SpecialFunctions
+#   LinearAlgebra
+#
+# 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.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 for (h,a)
+
+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
+```
+
+Worst case accuracy is about 2e-16.
+"""
+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 = (-T(0.9987710072524261), -T(0.9935301722663508), -T(0.9841245837228269), -T(0.9705915925462473), -T(0.9529877031604309), -T(0.9313866907065543), -T(0.9058791367155696)
+        , -T(0.8765720202742479), -T(0.8435882616243935), -T(0.8070662040294426), -T(0.7671590325157404), -T(0.7240341309238146), -T(0.6778723796326639), -T(0.6288673967765136)
+        , -T(0.5772247260839727), -T(0.5231609747222331), -T(0.4669029047509584), -T(0.4086864819907167), -T(0.34875588629216075), -T(0.28736248735545555), -T(0.22476379039468905)
+        , -T(0.1612223560688917), -T(0.0970046992094627), -T(0.03238017096286937), T(0.03238017096286937), T(0.0970046992094627), T(0.1612223560688917), T(0.22476379039468905)
+        , T(0.28736248735545555), T(0.34875588629216075), T(0.4086864819907167), T(0.4669029047509584), T(0.5231609747222331), T(0.5772247260839727), T(0.6288673967765136)
+        , T(0.6778723796326639), T(0.7240341309238146), T(0.7671590325157404), T(0.8070662040294426), T(0.8435882616243935), T(0.8765720202742479), T(0.9058791367155696)
+        , T(0.9313866907065543), T(0.9529877031604309), T(0.9705915925462473), T(0.9841245837228269), T(0.9935301722663508), T(0.9987710072524261))
+
+        w = (T(0.0031533460523059122), T(0.0073275539012762885), T(0.011477234579234613), T(0.015579315722943824), T(0.01961616045735561), T(0.023570760839324363)
+        , T(0.027426509708356944), T(0.031167227832798003), T(0.03477722256477052), T(0.038241351065830737), T(0.04154508294346467), T(0.0446745608566943), T(0.04761665849249045)
+        , T(0.05035903555385445), T(0.05289018948519363), T(0.055199503699984116), T(0.05727729210040322), T(0.05911483969839564), T(0.06070443916589387), T(0.06203942315989268)
+        , T(0.06311419228625402), T(0.06392423858464813), T(0.06446616443594998), T(0.06473769681268386), T(0.06473769681268386), T(0.06446616443594998), T(0.06392423858464813)
+        , T(0.06311419228625402), T(0.06203942315989268), T(0.06070443916589387), T(0.05911483969839564), T(0.05727729210040322), T(0.055199503699984116)
+        , T(0.05289018948519363), T(0.05035903555385445), T(0.04761665849249045), T(0.0446745608566943), T(0.04154508294346467), T(0.038241351065830737), T(0.03477722256477052)
+        , T(0.031167227832798003), T(0.027426509708356944), T(0.023570760839324363), T(0.01961616045735561), T(0.015579315722943824), T(0.011477234579234613), T(0.0073275539012762885)
+        , T(0.0031533460523059122))
+
+        towen = zero(a)
+        @inbounds for i in eachindex(w)
+            towen += w[i] * t2(h,a,x[i])
+        end
+
+        return towen
+    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)...)

test/owent.jl

@@ -0,0 +1,101 @@
+# Owen's T function tests
+
+using IrrationalConstants: inv2π, invsqrt2
+
+# test values for accurate and precise calculation
+hvec = [0.0625, 6.5, 7.0, 4.78125, 2.0, 1.0, 0.0625, 1, 1, 1, 1, 0.5, 0.5, 0.5, 0.5, 0.25, 0.25, 0.25, 0.25, 0.125, 0.125, 0.125, 0.125, 0.0078125
+, 0.0078125, 0.0078125, 0.0078125, 0.0078125, 0.0078125, 0.0625, 0.5, 0.9, 2.5, 7.33, 0.6, 1.6, 2.33, 2.33]
+avec = [0.25, 0.4375, 0.96875, 0.0625, 0.5, 0.9999975, 0.999999125, 0.5, 1, 2, 3, 0.5, 1, 2, 3, 0.5, 1, 2, 3, 0.5, 1, 2, 3, 0.5, 1, 2, 3, 10, 100
+, 0.999999999999999, 0.999999999999999, 0.999999999999999, 0.999999999999999, 0.999999999999999, 0.999999999999999, 0.999999999999999, 0.999999999999999
+, 0.99999]
+cvec = [big"0.0389119302347013668966224771378", big"2.00057730485083154100907167685e-11", big"6.399062719389853083219914429e-13"
+, big"1.06329748046874638058307112826e-7", big"0.00862507798552150713113488319155", big"0.0667418089782285927715589822405"
+, big"0.1246894855262192"
+, big"0.04306469112078537", big"0.06674188216570097", big"0.0784681869930841", big"0.0792995047488726", big"0.06448860284750375", big"0.1066710629614485"
+, big"0.1415806036539784", big"0.1510840430760184", big"0.07134663382271778", big"0.1201285306350883", big"0.1666128410939293", big"0.1847501847929859"
+, big"0.07317273327500386", big"0.1237630544953746", big"0.1737438887583106", big"0.1951190307092811", big"0.07378938035365545"
+, big"0.1249951430754052", big"0.1761984774738108", big"0.1987772386442824", big"0.2340886964802671", big"0.2479460829231492"
+, big"0.1246895548850743676554299881345328280176736760739893903915691894"
+, big"0.1066710629614484543187382775527753945264849005582264731161129477"
+, big"0.0750909978020473015080760056431386286348318447478899039422181015"
+, big"0.0030855526911589942124216949767707430484322201889568086810922629"
+, big"5.7538182971139187466647478665179637676531179007295252996453e-14", big"0.0995191725772188724714794470740785702586033387786949658229016920"
+, big"0.0258981646643923680014142514442989928165349517076730515952020227"
+, big"0.0049025023268168675126146823752680242063832053551244071400100690"
+, big"0.0049024988349089450612896251009169062698683918433614542387524648"]
+
+# test values for type stability checking
+ht = [0.0625, 0.0, 0.5, 0.5, 0.5]
+at = [0.025, 0.5, 0.0, 1.0, +Inf]
+
+#=
+# Interactive tests, displayin error and error magnitude
+mvbck = "\033[1A\033[58G"
+mva = "\033[72G"
+mve = "\033[90G"
+
+warmup = owent(0.0625,0.025)
+warmup = owent(0.0625,0.999999125)
+
+for i in 1:size(hvec,1)
+    if i == 1
+        println("\t\tExecution Time","\033[58G","h",mva,"a",mve,"error\t\t","log10 error")
+    end
+    h = hvec[i]
+    a = avec[i]
+    c = cvec[i]
+    print(i,"\t")
+    @time tx = owent(h,a)
+    err = Float64(tx-c)
+    logerr = Float64(round(log10(abs(tx-c)),sigdigits=3))
+    println(mvbck,h,mva,a,mve,err,"\t",logerr)
+end
+=#
+
+@testset "Owen's T" begin
+
+    # check that error for calc is within specification
+    @testset "Owen's T value checks" begin
+
+        for i in 1:size(hvec,1)
+            h = hvec[i] # h test value
+            a = avec[i] # a test value
+            c = cvec[i] # the "correct" answer
+            t = owent(h,a)
+
+            err = round(log10(abs(t-c)),digits=0)
+
+            @test err ≤ -16.0
+        end
+
+    end
+
+    @testset "Owen's T Shortcut Evaluations" begin
+        @test owent(-0.0625,0.025) == owent(0.0625,0.025) # if h<0 equals owent(abs(h),a)
+        @test owent(0.0,0.025) == atan(0.025)*inv2π  # if h=0, t = atan(a)*inv2π
+        @test owent(0.0625,0.0) == 0.0  # when a=0, owent=0
+        @test owent(0.0625,-0.025) == -owent(0.0625,0.025) # if a<0, t = -owent(h,abs(a))
+        @test owent(0.0625,1.0) == 0.125*erfc(-0.0625*invsqrt2)*erfc(0.0625*invsqrt2) # if a=1, t = (formula)
+        @test owent(0.0625,+Inf) == 0.25*erfc(sqrt(0.0625^2)*invsqrt2) # if a=∞, t = (formula)
+        @test owent(0.0625,-Inf) == -0.25*erfc(sqrt(0.0625^2)*invsqrt2) # should also work, due to a<0 condition
+    end
+
+    @testset "Owen's T type stability" begin
+
+        for T1 in [Float16, Float32, Float64, BigFloat]
+            for T2 in [BigFloat, Float64, Float32, Float16]
+                for i in 1:size(ht,1)
+                    h=T1(ht[i])
+                    a=T2(at[i])
+                    (p1, p2) = promote(h,a)
+                    t=owent(h,a)
+                    @test typeof(t) == typeof(p1)
+                    #println("T1: ",T1," ",typeof(p1),"\tT2: ",T2," ",typeof(p2),"\t",i,"\th ",h,"\ta ",a,"\tt ",t,"\t",typeof(t)," ",typeof(t)==typeof(p1))
+                end
+            end
+        end
+
+    end # test type stability
+
+end
+# Owen's T Function Tests

test/runtests.jl

@@ -33,7 +33,8 @@ tests = [
     "logabsgamma",
     "sincosint",
     "other_tests",
-    "chainrules"
+    "chainrules",
+    "owent"
 ]
 
 const testdir = dirname(@__FILE__)