From 8e14c19b681bbb679fa851a2cac3aaccfa7e1fbc Mon Sep 17 00:00:00 2001 From: Valentin Kaisermayer Date: Fri, 25 Sep 2020 15:55:37 +0200 Subject: [PATCH 1/3] Adds Owen's T function --- src/SpecialFunctions.jl | 9 +++--- src/owens_t.jl | 61 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 66 insertions(+), 4 deletions(-) create mode 100644 src/owens_t.jl diff --git a/src/SpecialFunctions.jl b/src/SpecialFunctions.jl index 814e2a46..d19b1191 100644 --- a/src/SpecialFunctions.jl +++ b/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 diff --git a/src/owens_t.jl b/src/owens_t.jl new file mode 100644 index 00000000..c633f3d5 --- /dev/null +++ b/src/owens_t.jl @@ -0,0 +1,61 @@ +""" + owens_t(x, a) + +Owen's T function + +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 From 0e485c28d41d8e432500bf2866784881d67336a1 Mon Sep 17 00:00:00 2001 From: Valentin Kaisermayer Date: Fri, 25 Sep 2020 16:13:06 +0200 Subject: [PATCH 2/3] added testset for Owen's T function --- test/owens_t.jl | 6 ++++++ test/runtests.jl | 9 +++++---- 2 files changed, 11 insertions(+), 4 deletions(-) create mode 100644 test/owens_t.jl diff --git a/test/owens_t.jl b/test/owens_t.jl new file mode 100644 index 00000000..dd83a0ec --- /dev/null +++ b/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: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 \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 786f2505..e87d0708 100644 --- a/test/runtests.jl +++ b/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 From 225521a95a65f1c5a395e0f39fa4613ca5eb55d7 Mon Sep 17 00:00:00 2001 From: Valentin Kaisermayer Date: Fri, 25 Sep 2020 17:05:50 +0200 Subject: [PATCH 3/3] minor change added comment --- src/owens_t.jl | 2 ++ test/owens_t.jl | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/owens_t.jl b/src/owens_t.jl index c633f3d5..fb99ffc8 100644 --- a/src/owens_t.jl +++ b/src/owens_t.jl @@ -3,6 +3,8 @@ 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, diff --git a/test/owens_t.jl b/test/owens_t.jl index dd83a0ec..996abf0e 100644 --- a/test/owens_t.jl +++ b/test/owens_t.jl @@ -1,6 +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:3), 1 ./ (2 .* pi) .* atan.(-3:0.1:3), rtol=1e-6)) + @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 \ No newline at end of file