Adds Owen's T function

src/SpecialFunctions.jl

@@ -2,8 +2,7 @@ module SpecialFunctions
 
 using OpenSpecFun_jll
 
-export
-    airyai,
+export airyai,
     airyaiprime,
     airybi,
     airybiprime,
@@ -57,7 +56,8 @@ export
     zeta,
     sinint,
     cosint,
-    lbinomial
+    lbinomial,
+    owens_t
 
 include("bessel.jl")
 include("erf.jl")
@@ -68,12 +68,13 @@ include("gamma_inc.jl")
 include("betanc.jl")
 include("beta_inc.jl")
 include("deprecated.jl")
+include("owens_t.jl")
 
 for f in (:digamma, :erf, :erfc, :erfcinv, :erfcx, :erfi, :erfinv, :logerfc, :logerfcx,
           :eta, :gamma, :invdigamma, :logfactorial, :lgamma, :trigamma, :ellipk, :ellipe)
     @eval $(f)(::Missing) = missing
 end
-for f in (:beta, :lbeta)
+for f in (:beta, :lbeta, :owens_t)
     @eval $(f)(::Number, ::Missing) = missing
     @eval $(f)(::Missing, ::Number) = missing
     @eval $(f)(::Missing, ::Missing) = missing

src/owens_t.jl

@@ -0,0 +1,63 @@
+"""
+    owens_t(x, a)
+
+Owen's T function
+
+Ported from Matlab implementation of John Burkardt, which is distributed under the GNU LGPL license.
+
+Reference:
+    MA Porter, DJ Winstanley, Remark AS R30: A Remark on Algorithm AS76: An Integral Useful in 
+    Calculating Noncentral T and Bivariate Normal Probabilities, Applied Statistics, Volume 28, 
+    Number 1, 1979, page 113.
+
+    JC Young, Christoph Minder, Algorithm AS 76: An Algorithm Useful in Calculating Non-Central T and
+    Bivariate Normal Distributions, Applied Statistics, Volume 23, Number 3, 1974, pages 455-457.
+"""
+function owens_t(x::Real, a::Real)
+    ng = 5
+    r = [0.1477621, 0.1346334, 0.1095432, 0.0747257, 0.0333357]
+    tp = 0.159155
+    tv1 = 1.0E-35
+    tv2 = 15.0
+    tv3 = 15.0
+    tv4 = 1.0E-05
+    u = [0.0744372, 0.2166977, 0.3397048, 0.4325317, 0.4869533]
+    # Test for X near zero.
+    if abs(x) < tv1
+        return tp * atan(a)
+    end
+    # Test for large values of abs(X).
+    if tv2 < abs(x)
+        return 0.
+    end
+    # Test for `a` near zero.
+    if abs(a) < tv1
+        return 0.
+    end
+    # Test whether abs(`a`) is so large that it must be truncated.
+    xs = -0.5 * x * x
+    x2 = a
+    fxs = a * a
+    #  Computation of truncation point by Newton iteration.  
+    if tv3 <= log(1.0 + fxs) - xs * fxs
+        x1 = 0.5 * a
+        fxs = 0.25 * fxs
+        while true
+            rt = fxs + 1.
+            x2 = x1 + (xs * fxs + tv3 - log(rt)) / (2. * x1 * ( 1. / rt - xs))
+            fxs = x2 * x2
+            if abs(x2 - x1) < tv4
+                break
+            end
+            x1 = x2
+        end
+    end
+    # Gaussian quadrature.
+    rt = 0.0
+    for i in 1:ng
+        r1 = 1.0 + fxs * (0.5 + u[i])^2
+        r2 = 1.0 + fxs * (0.5 - u[i])^2
+        rt = rt + r[i] * (exp(xs * r1) / r1 + exp(xs * r2) / r2)
+    end
+    return rt * x2 * tp
+end

test/owens_t.jl

@@ -0,0 +1,6 @@
+@testset "Owen's T function" begin
+    @test all(owens_t.(-3:0.1:3, 0) .== 0.)
+    @test all(isapprox.(owens_t.(0, -3:0.1:3), 1 ./ (2 .* pi) .* atan.(-3:0.1:3), rtol=1e-6))
+    @test all(isapprox.([owens_t(-h, a) for h in -3:0.1:3, a in -3:0.1:3], [owens_t(h, a) for h in -3:0.1:3, a in -3:0.1:3], rtol=1e-6))
+    @test all(isapprox.([owens_t(h, -a) for h in -3:0.1:3, a in -3:0.1:3], [-owens_t(h, a) for h in -3:0.1:3, a in -3:0.1:3], rtol=1e-6))
+end

test/runtests.jl

@@ -9,8 +9,8 @@ using SpecialFunctions: AmosException, f64
 # useful test functions for relative error, which differ from isapprox
 # in that relerrc separately looks at the real and imaginary parts
 relerr(z, x) = z == x ? 0.0 : abs(z - x) / abs(x)
-relerrc(z, x) = max(relerr(real(z),real(x)), relerr(imag(z),imag(x)))
-≅(a,b) = relerrc(a,b) ≤ 1e-13
+relerrc(z, x) = max(relerr(real(z), real(x)), relerr(imag(z), imag(x)))
+≅(a,b) = relerrc(a, b) ≤ 1e-13
 
 tests = [
     "bessel",
@@ -21,7 +21,8 @@ tests = [
     "gamma_inc",
     "gamma",
     "sincosint",
-    "other_tests"
+    "other_tests",
+    "owens_t"
 ]
 
 const testdir = dirname(@__FILE__)
@@ -30,6 +31,6 @@ const testdir = dirname(@__FILE__)
 for t in tests
     tp = joinpath(testdir, "$(t).jl")
     @testset "$(t)" begin
-      include(tp)
+        include(tp)
     end
 end