diff --git a/src/SpecialFunctions.jl b/src/SpecialFunctions.jl index 0d296cc0..10927697 100644 --- a/src/SpecialFunctions.jl +++ b/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 diff --git a/src/owent.jl b/src/owent.jl new file mode 100644 index 00000000..f048e45d --- /dev/null +++ b/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") diff --git a/test/owent.jl b/test/owent.jl new file mode 100644 index 00000000..4a9347f0 --- /dev/null +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index 69dcafe7..197b9825 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -33,7 +33,8 @@ tests = [ "logabsgamma", "sincosint", "other_tests", - "chainrules" + "chainrules", + "owent" ] const testdir = dirname(@__FILE__)