Provides same interface to studentized range distribution functions ([p/q]tukey)

REQUIRE

@@ -1,3 +1,4 @@
 julia 0.7-beta
 Rmath 0.4.0
 SpecialFunctions
+QuadGK

src/StatsFuns.jl

@@ -217,6 +217,18 @@ export
     tdistinvlogcdf,     # inverse-logcdf of student's t distribution
     tdistinvlogccdf,    # inverse-logccdf of student's t distribution
 
+    # distrs/srdist
+    srdistpdf,           # pdf of studentized range distribution
+    srdistlogpdf,        # logpdf of studentized range distribution
+    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 +256,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

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)

src/distrs/srdist.jl

@@ -0,0 +1,29 @@
+# functions related to studentized range distribution
+
+using QuadGK: quadgk
+
+import .RFunctions:
+    srdistcdf,
+    srdistccdf,
+    srdistlogcdf,
+    srdistlogccdf,
+    srdistinvcdf,
+    srdistinvccdf,
+    srdistinvlogcdf,
+    srdistinvlogccdf
+
+# pdf for numbers with generic types
+function srdistpdf(ν::Real, k::Real, x::Number)
+    Φ(x) = normcdf(x)
+    ϕ(x) = normpdf(x)
+    function outer(y)
+        inner(u) = ϕ(u) * ϕ(u - x*y) * (Φ(u) - Φ(u - x*y))^(k-2)
+        inner_part = quadgk(inner, -Inf, Inf) |> first
+        return inner_part * y^ν * ϕ(y*√ν)
+    end
+    integral = quadgk(outer, 0.0, Inf) |> first
+    return integral * (√(2π) * k * (k-1) * ν^(ν/2)) / (gamma(ν/2) * 2^(ν/2 - 1))
+end
+
+# logpdf for numbers with generic types
+srdistlogpdf(ν::Real, k::Real, x::Number) = log(srdistpdf(ν, k, x))

src/rmath.jl

@@ -34,9 +34,14 @@ function _import_rmath(rname::Symbol, jname::Symbol, pargs)
     invlogcdf = Symbol(jname, "invlogcdf")
     invlogccdf = Symbol(jname, "invlogccdf")
 
+    is_tukey = false
+    if rname == :tukey
+        is_tukey = true
+    end
+
     rand = Symbol(jname, "rand")
     has_rand = true
-    if rname == :nbeta || rname == :nf || rname == :nt
+    if rname == :nbeta || rname == :nf || rname == :nt || rname == :tukey
         has_rand = false
     end
 
@@ -46,45 +51,75 @@ function _import_rmath(rname::Symbol, jname::Symbol, pargs)
     dtypes = Expr(:tuple, :Cdouble, _pts..., :Cint)
     ptypes = Expr(:tuple, :Cdouble, _pts..., :Cint, :Cint)
     qtypes = Expr(:tuple, :Cdouble, _pts..., :Cint, :Cint)
+    tukeytypes = Expr(:tuple, :Cdouble, :Cdouble, _pts..., :Cint, :Cint)
     rtypes = Expr(:tuple, _pts...)
 
     pdecls = [Expr(:(::), ps, :(Union{Float64,Int})) for ps in pargs] # [:(p1::Union{Float64, Int}), :(p2::Union{...}), ...]
 
-    # function implementation
-    quote
-        $pdf($(pdecls...), x::Union{Float64,Int}) =
-            ccall(($dfun, libRmath), Float64, $dtypes, x, $(pargs...), 0)
+    # Function implementation
+    if is_tukey
+        pargs = reverse(pargs)
+        quote
+            $cdf($(pdecls...), x::Union{Float64,Int}) =
+                ccall(($pfun, libRmath), Float64, $tukeytypes, x, 1, $(pargs...), 1, 0)
+
+            $ccdf($(pdecls...), x::Union{Float64,Int}) =
+                ccall(($pfun, libRmath), Float64, $tukeytypes, x, 1, $(pargs...), 0, 0)
+
+            $logcdf($(pdecls...), x::Union{Float64,Int}) =
+                ccall(($pfun, libRmath), Float64, $tukeytypes, x, 1, $(pargs...), 1, 1)
+
+            $logccdf($(pdecls...), x::Union{Float64,Int}) =
+                ccall(($pfun, libRmath), Float64, $tukeytypes, x, 1, $(pargs...), 0, 1)
+
+            $invcdf($(pdecls...), q::Union{Float64,Int}) =
+                ccall(($qfun, libRmath), Float64, $tukeytypes, q, 1, $(pargs...), 1, 0)
+
+            $invccdf($(pdecls...), q::Union{Float64,Int}) =
+                ccall(($qfun, libRmath), Float64, $tukeytypes, q, 1, $(pargs...), 0, 0)
+
+            $invlogcdf($(pdecls...), lq::Union{Float64,Int}) =
+                ccall(($qfun, libRmath), Float64, $tukeytypes, lq, 1, $(pargs...), 1, 1)
+
+            $invlogccdf($(pdecls...), lq::Union{Float64,Int}) =
+                ccall(($qfun, libRmath), Float64, $tukeytypes, lq, 1, $(pargs...), 0, 1)
+        end
+    else
+        quote
+            $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)
 
-        $cdf($(pdecls...), x::Union{Float64,Int}) =
-            ccall(($pfun, libRmath), Float64, $ptypes, x, $(pargs...), 1, 0)
+            $cdf($(pdecls...), x::Union{Float64,Int}) =
+                ccall(($pfun, libRmath), Float64, $ptypes, x, $(pargs...), 1, 0)
 
-        $ccdf($(pdecls...), x::Union{Float64,Int}) =
-            ccall(($pfun, libRmath), Float64, $ptypes, x, $(pargs...), 0, 0)
+            $ccdf($(pdecls...), x::Union{Float64,Int}) =
+                ccall(($pfun, libRmath), Float64, $ptypes, x, $(pargs...), 0, 0)
 
-        $logcdf($(pdecls...), x::Union{Float64,Int}) =
-            ccall(($pfun, libRmath), Float64, $ptypes, x, $(pargs...), 1, 1)
+            $logcdf($(pdecls...), x::Union{Float64,Int}) =
+                ccall(($pfun, libRmath), Float64, $ptypes, x, $(pargs...), 1, 1)
 
-        $logccdf($(pdecls...), x::Union{Float64,Int}) =
-            ccall(($pfun, libRmath), Float64, $ptypes, x, $(pargs...), 0, 1)
+            $logccdf($(pdecls...), x::Union{Float64,Int}) =
+                ccall(($pfun, libRmath), Float64, $ptypes, x, $(pargs...), 0, 1)
 
-        $invcdf($(pdecls...), q::Union{Float64,Int}) =
-            ccall(($qfun, libRmath), Float64, $qtypes, q, $(pargs...), 1, 0)
+            $invcdf($(pdecls...), q::Union{Float64,Int}) =
+                ccall(($qfun, libRmath), Float64, $qtypes, q, $(pargs...), 1, 0)
 
-        $invccdf($(pdecls...), q::Union{Float64,Int}) =
-            ccall(($qfun, libRmath), Float64, $qtypes, q, $(pargs...), 0, 0)
+            $invccdf($(pdecls...), q::Union{Float64,Int}) =
+                ccall(($qfun, libRmath), Float64, $qtypes, q, $(pargs...), 0, 0)
 
-        $invlogcdf($(pdecls...), lq::Union{Float64,Int}) =
-            ccall(($qfun, libRmath), Float64, $qtypes, lq, $(pargs...), 1, 1)
+            $invlogcdf($(pdecls...), lq::Union{Float64,Int}) =
+                ccall(($qfun, libRmath), Float64, $qtypes, lq, $(pargs...), 1, 1)
 
-        $invlogccdf($(pdecls...), lq::Union{Float64,Int}) =
-            ccall(($qfun, libRmath), Float64, $qtypes, lq, $(pargs...), 0, 1)
+            $invlogccdf($(pdecls...), lq::Union{Float64,Int}) =
+                ccall(($qfun, libRmath), Float64, $qtypes, lq, $(pargs...), 0, 1)
 
-        if $has_rand
-            $rand($(pdecls...)) =
-                ccall(($rfun, libRmath), Float64, $rtypes, $(pargs...))
+            if $has_rand
+                $rand($(pdecls...)) =
+                    ccall(($rfun, libRmath), Float64, $rtypes, $(pargs...))
+            end
         end
     end
 end
@@ -99,13 +134,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 λ

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