Merge pull request #66 from BioTurboNick/master

Provides same interface to studentized range distribution functions ([p/q]tukey)
JuliaStats · Feb 4, 2019 · b46dadc · b46dadc
2 parents 450a435 + 259bcda
commit b46dadc
Show file tree

Hide file tree

Showing 6 changed files with 296 additions and 215 deletions.
diff --git a/src/StatsFuns.jl b/src/StatsFuns.jl
@@ -217,6 +217,16 @@ export
     tdistinvlogcdf,     # inverse-logcdf of student's t distribution
     tdistinvlogccdf,    # inverse-logccdf of student's t distribution
 
+    # distrs/srdist
+    srdistcdf,           # cdf of studentized range distribution
+    srdistccdf,          # ccdf of studentized range distribution
+    srdistlogcdf,        # logcdf of studentized range distribution
+    srdistlogccdf,       # logccdf of studentized range distribution
+    srdistinvcdf,        # inverse-cdf of studentized range distribution
+    srdistinvccdf,       # inverse-ccdf of studentized range distribution
+    srdistinvlogcdf,     # inverse-logcdf of studentized range distribution
+    srdistinvlogccdf,    # inverse-logccdf of studentized range distribution
+
     # misc
     logmvgamma,         # logarithm of multivariate gamma function
     lstirling_asym
@@ -244,5 +254,6 @@ include(joinpath("distrs", "norm.jl"))
 include(joinpath("distrs", "ntdist.jl"))
 include(joinpath("distrs", "pois.jl"))
 include(joinpath("distrs", "tdist.jl"))
+include(joinpath("distrs", "srdist.jl"))
 
 end # module
diff --git a/src/distrs/fdist.jl b/src/distrs/fdist.jl
@@ -13,7 +13,7 @@ import .RFunctions:
     fdistinvlogccdf
 
 # pdf for numbers with generic types
-fdistpdf(d1::Real, d2::Real, x::Number) = sqrt((d1 * x)^d1 * d2^d2 / (d1 * x + d2)^(d1 + d2)) / (x * beta(d1 / 2, d2 / 2))
+fdistpdf(ν1::Real, ν2::Real, x::Number) = sqrt((ν1 * x)^ν1 * ν2^ν2 / (ν1 * x + ν2)^(ν1 + ν2)) / (x * beta(ν1 / 2, ν2 / 2))
 
 # logpdf for numbers with generic types
-fdistlogpdf(d1::Real, d2::Real, x::Number) = (d1 * log(d1 * x) + d2 * log(d2) - (d1 + d2) * log(d1 * x + d2)) / 2 - log(x) - lbeta(d1 / 2, d2 / 2)
+fdistlogpdf(ν1::Real, ν2::Real, x::Number) = (ν1 * log(ν1 * x) + ν2 * log(ν2) - (ν1 + ν2) * log(ν1 * x + ν2)) / 2 - log(x) - lbeta(ν1 / 2, ν2 / 2)
diff --git a/src/distrs/srdist.jl b/src/distrs/srdist.jl
@@ -0,0 +1,11 @@
+# functions related to studentized range distribution
+
+import .RFunctions:
+    srdistcdf,
+    srdistccdf,
+    srdistlogcdf,
+    srdistlogccdf,
+    srdistinvcdf,
+    srdistinvccdf,
+    srdistinvlogcdf,
+    srdistinvlogccdf
diff --git a/src/rmath.jl b/src/rmath.jl
@@ -34,29 +34,45 @@ function _import_rmath(rname::Symbol, jname::Symbol, pargs)
     invlogcdf = Symbol(jname, "invlogcdf")
     invlogccdf = Symbol(jname, "invlogccdf")
 
+    is_tukey = rname == :tukey
+
     rand = Symbol(jname, "rand")
     has_rand = true
-    if rname == :nbeta || rname == :nf || rname == :nt
+    if rname == :nbeta || rname == :nf || rname == :nt || is_tukey
         has_rand = false
     end
 
     # arguments & argument types
     _pts = fill(:Cdouble, length(pargs))
 
-    dtypes = Expr(:tuple, :Cdouble, _pts..., :Cint)
-    ptypes = Expr(:tuple, :Cdouble, _pts..., :Cint, :Cint)
-    qtypes = Expr(:tuple, :Cdouble, _pts..., :Cint, :Cint)
-    rtypes = Expr(:tuple, _pts...)
+    if is_tukey
+        dtypes = Expr(:tuple, :Cdouble, :Cdouble, _pts..., :Cint)
+        ptypes = Expr(:tuple, :Cdouble, :Cdouble, _pts..., :Cint, :Cint)
+        qtypes = Expr(:tuple, :Cdouble, :Cdouble, _pts..., :Cint, :Cint)
+        rtypes = Expr(:tuple, :Cdouble, _pts...)
+    else
+        dtypes = Expr(:tuple, :Cdouble, _pts..., :Cint)
+        ptypes = Expr(:tuple, :Cdouble, _pts..., :Cint, :Cint)
+        qtypes = Expr(:tuple, :Cdouble, _pts..., :Cint, :Cint)
+        rtypes = Expr(:tuple, _pts...)
+    end
 
     pdecls = [Expr(:(::), ps, :(Union{Float64,Int})) for ps in pargs] # [:(p1::Union{Float64, Int}), :(p2::Union{...}), ...]
 
-    # function implementation
+    if is_tukey
+        # ptukey and qtukey have an extra literal 1 argument
+        pargs = (1, pargs...)
+    end
+
+    # Function implementation
     quote
-        $pdf($(pdecls...), x::Union{Float64,Int}) =
-            ccall(($dfun, libRmath), Float64, $dtypes, x, $(pargs...), 0)
+        if $(!is_tukey)
+            $pdf($(pdecls...), x::Union{Float64,Int}) =
+                ccall(($dfun, libRmath), Float64, $dtypes, x, $(pargs...), 0)
 
-        $logpdf($(pdecls...), x::Union{Float64,Int}) =
-            ccall(($dfun, libRmath), Float64, $dtypes, x, $(pargs...), 1)
+            $logpdf($(pdecls...), x::Union{Float64,Int}) =
+                ccall(($dfun, libRmath), Float64, $dtypes, x, $(pargs...), 1)
+        end
 
         $cdf($(pdecls...), x::Union{Float64,Int}) =
             ccall(($pfun, libRmath), Float64, $ptypes, x, $(pargs...), 1, 0)
@@ -99,13 +115,14 @@ end
 @import_rmath beta beta α β
 @import_rmath binom binom n p
 @import_rmath chisq chisq k
-@import_rmath f fdist d1 d2
+@import_rmath f fdist ν1 ν2
 @import_rmath gamma gamma α β
 @import_rmath hyper hyper ms mf n
 @import_rmath norm norm μ σ
 @import_rmath nbinom nbinom r p
 @import_rmath pois pois λ
 @import_rmath t tdist k
+@import_rmath tukey srdist k ν
 
 @import_rmath nbeta nbeta α β λ
 @import_rmath nchisq nchisq k λ

diff --git a/test/generic.jl b/test/generic.jl
@@ -43,85 +43,95 @@ end
 
 ### Test cases
 
-genericcomp_tests("beta", [
-    ((1.0, 1.0), 0.01:0.01:0.99),
-    ((2.0, 3.0), 0.01:0.01:0.99),
-    ((10.0, 2.0), 0.01:0.01:0.99),
-])
-
-genericcomp_tests("binom", [
-    ((1, 0.5), 0.0:1.0),
-    ((1, 0.7), 0.0:1.0),
-    ((8, 0.6), 0.0:8.0),
-    ((20, 0.1), 0.0:20.0),
-    ((20, 0.9), 0.0:20.0),
-    ((20, 0.9), 0:20),
-])
-
-genericcomp_tests("chisq", [
-    ((1,), 0.0:0.1:8.0),
-    ((4,), 0.0:0.1:8.0),
-    ((9,), 0.0:0.1:8.0),
-])
-
-genericcomp_tests("fdist", [
-    ((1, 1), (0.0:0.1:5.0)),
-    ((2, 1), (0.0:0.1:5.0)),
-    ((5, 2), (0.0:0.1:5.0)),
-    ((10, 1), (0.0:0.1:5.0)),
-    ((10, 3), (0.0:0.1:5.0)),
-])
-
-genericcomp_tests("gamma", [
-    ((1.0, 1.0), (0.05:0.05:12.0)),
-    ((0.5, 1.0), (0.05:0.05:12.0)),
-    ((3.0, 1.0), (0.05:0.05:12.0)),
-    ((9.0, 1.0), (0.05:0.05:12.0)),
-    ((2.0, 3.0), (0.05:0.05:12.0)),
-])
-
-# genericcomp_tests("nbeta", [
-#     ((1.0, 1.0, 0.0), 0.01:0.01:0.99),
-#     ((2.0, 3.0, 0.0), 0.01:0.01:0.99),
-#     ((1.0, 1.0, 2.0), 0.01:0.01:0.99),
-#     ((3.0, 4.0, 2.0), 0.01:0.01:0.99),
-# ])
-
-# genericcomp_tests("nchisq", [
-#     ((2, 1), 0.0:0.2:8.0),
-#     ((2, 3), 0.0:0.2:8.0),
-#     ((4, 1), 0.0:0.2:8.0),
-#     ((4, 3), 0.0:0.2:8.0),
-# ])
-
-# genericcomp_tests("nfdist", [
-#     ((1.0, 1.0, 0.0), 0.1:0.1:10.0),
-#     ((1.0, 1.0, 2.0), 0.1:0.1:10.0),
-#     ((2.0, 3.0, 1.0), 0.1:0.1:10.0),
-# ])
-
-genericcomp_tests("norm", [
-    ((0.0, 1.0), -6.0:0.01:6.0),
-    ((2.0, 1.0), -3.0:0.01:7.0),
-    ((0.0, 0.5), -3.0:0.01:3.0),
-    ((0, 1), -3.0:3.0),
-    ((0, 1), -3.0:0.01:3.0)
-])
-
-# genericcomp_tests("ntdist", [
-#     ((0, 1), -4.0:0.1:10.0),
-#     ((0, 4), -4.0:0.1:10.0),
-#     ((2, 1), -4.0:0.1:10.0),
-#     ((2, 4), -4.0:0.1:10.0),
-# ])
-
-genericcomp_tests("pois", [
-    ((1.0,), 0:30),
-    ((10.0,), 0:42)
-])
-
-genericcomp_tests("tdist", [
-    ((1,), -5.0:0.1:5.0),
-    ((2,), -5.0:0.1:5.0),
-    ((5,), -5.0:0.1:5.0),
-])
+@testset "Generic" begin
+    genericcomp_tests("beta", [
+        ((1.0, 1.0), 0.01:0.01:0.99),
+        ((2.0, 3.0), 0.01:0.01:0.99),
+        ((10.0, 2.0), 0.01:0.01:0.99),
+    ])
+
+    genericcomp_tests("binom", [
+        ((1, 0.5), 0.0:1.0),
+        ((1, 0.7), 0.0:1.0),
+        ((8, 0.6), 0.0:8.0),
+        ((20, 0.1), 0.0:20.0),
+        ((20, 0.9), 0.0:20.0),
+        ((20, 0.9), 0:20),
+    ])
+
+    genericcomp_tests("chisq", [
+        ((1,), 0.0:0.1:8.0),
+        ((4,), 0.0:0.1:8.0),
+        ((9,), 0.0:0.1:8.0),
+    ])
+
+    genericcomp_tests("fdist", [
+        ((1, 1), (0.0:0.1:5.0)),
+        ((2, 1), (0.0:0.1:5.0)),
+        ((5, 2), (0.0:0.1:5.0)),
+        ((10, 1), (0.0:0.1:5.0)),
+        ((10, 3), (0.0:0.1:5.0)),
+    ])
+
+    genericcomp_tests("gamma", [
+        ((1.0, 1.0), (0.05:0.05:12.0)),
+        ((0.5, 1.0), (0.05:0.05:12.0)),
+        ((3.0, 1.0), (0.05:0.05:12.0)),
+        ((9.0, 1.0), (0.05:0.05:12.0)),
+        ((2.0, 3.0), (0.05:0.05:12.0)),
+    ])
+
+    # genericcomp_tests("nbeta", [
+    #     ((1.0, 1.0, 0.0), 0.01:0.01:0.99),
+    #     ((2.0, 3.0, 0.0), 0.01:0.01:0.99),
+    #     ((1.0, 1.0, 2.0), 0.01:0.01:0.99),
+    #     ((3.0, 4.0, 2.0), 0.01:0.01:0.99),
+    # ])
+
+    # genericcomp_tests("nchisq", [
+    #     ((2, 1), 0.0:0.2:8.0),
+    #     ((2, 3), 0.0:0.2:8.0),
+    #     ((4, 1), 0.0:0.2:8.0),
+    #     ((4, 3), 0.0:0.2:8.0),
+    # ])
+
+    # genericcomp_tests("nfdist", [
+    #     ((1.0, 1.0, 0.0), 0.1:0.1:10.0),
+    #     ((1.0, 1.0, 2.0), 0.1:0.1:10.0),
+    #     ((2.0, 3.0, 1.0), 0.1:0.1:10.0),
+    # ])
+
+    genericcomp_tests("norm", [
+        ((0.0, 1.0), -6.0:0.01:6.0),
+        ((2.0, 1.0), -3.0:0.01:7.0),
+        ((0.0, 0.5), -3.0:0.01:3.0),
+        ((0, 1), -3.0:3.0),
+        ((0, 1), -3.0:0.01:3.0)
+    ])
+
+    # genericcomp_tests("ntdist", [
+    #     ((0, 1), -4.0:0.1:10.0),
+    #     ((0, 4), -4.0:0.1:10.0),
+    #     ((2, 1), -4.0:0.1:10.0),
+    #     ((2, 4), -4.0:0.1:10.0),
+    # ])
+
+    genericcomp_tests("pois", [
+        ((1.0,), 0:30),
+        ((10.0,), 0:42)
+    ])
+
+    genericcomp_tests("tdist", [
+        ((1,), -5.0:0.1:5.0),
+        ((2,), -5.0:0.1:5.0),
+        ((5,), -5.0:0.1:5.0),
+    ])
+
+    #genericcomp_tests("srdist", [
+    #    ((1,2), (0.0:0.1:5.0)),
+    #    ((2,2), (0.0:0.1:5.0)),
+    #    ((5,3), (0.0:0.1:5.0)),
+    #    ((10,2), (0.0:0.1:5.0)),
+    #    ((10,5), (0.0:0.1:5.0))
+    #])
+end