From 29b4076665f92ef3edfc4bcd8c4609a6588e4f72 Mon Sep 17 00:00:00 2001
From: Nicholas Bauer <nicholasbauer@outlook.com>
Date: Thu, 24 Jan 2019 23:46:07 -0500
Subject: [PATCH 01/13] Added studentized range distribution (srdist)

pdf code was written originally by user "Doomphoenix Qxz" on discourse.julialang.org
---
 REQUIRE              |  2 ++
 src/StatsFuns.jl     | 13 +++++++++++++
 src/distrs/fdist.jl  |  4 ++--
 src/distrs/srdist.jl | 32 ++++++++++++++++++++++++++++++++
 src/rmath.jl         | 20 ++++++++++++++------
 5 files changed, 63 insertions(+), 8 deletions(-)
 create mode 100644 src/distrs/srdist.jl

diff --git a/REQUIRE b/REQUIRE
index 02ca4f6..aa33dd5 100644
--- a/REQUIRE
+++ b/REQUIRE
@@ -1,3 +1,5 @@
 julia 0.7-beta
 Rmath 0.4.0
 SpecialFunctions
+QuadGK
+ForwardDiff
diff --git a/src/StatsFuns.jl b/src/StatsFuns.jl
index 41a8a5d..f5a44cb 100644
--- a/src/StatsFuns.jl
+++ b/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
diff --git a/src/distrs/fdist.jl b/src/distrs/fdist.jl
index b0b9fdc..b46c796 100644
--- 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
new file mode 100644
index 0000000..3cf10c9
--- /dev/null
+++ b/src/distrs/srdist.jl
@@ -0,0 +1,32 @@
+# functions related to studentized range distribution
+
+using ForwardDiff: derivative
+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) = (1+erf(x / √2)) / 2
+    ϕ(x) = derivative(Φ, x)
+    function outer(y)
+        function inner(u)
+            return ϕ(u) * ϕ(u - x*y) * (Φ(u) - Φ(u - x*y))^(k-2)
+        end
+        inner_part = quadgk(inner, -Inf, Inf)[1]
+        return inner_part * y^ν * ϕ(y*√ν)
+    end
+    integral = quadgk(outer, 0.0, Inf)[1]
+    return integral * (√τ * k * (k-1) * ν^(ν/2)) / (gamma(ν/2) * 2^(ν/2 - 1))
+end
+
+# logpdf for numbers with generic types
+srdistlogpdf(ν::Real, k::Real, x::Number) = log(pdf(ν, k, x))
diff --git a/src/rmath.jl b/src/rmath.jl
index bb85bdb..9de5982 100644
--- a/src/rmath.jl
+++ b/src/rmath.jl
@@ -34,9 +34,14 @@ function _import_rmath(rname::Symbol, jname::Symbol, pargs)
     invlogcdf = Symbol(jname, "invlogcdf")
     invlogccdf = Symbol(jname, "invlogccdf")
 
+    has_density = true
+    if rname == :tukey
+        has_density = false
+    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
 
@@ -52,11 +57,13 @@ function _import_rmath(rname::Symbol, jname::Symbol, pargs)
 
     # function implementation
     quote
-        $pdf($(pdecls...), x::Union{Float64,Int}) =
-            ccall(($dfun, libRmath), Float64, $dtypes, x, $(pargs...), 0)
+        if $has_density
+            $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 +106,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 λ

From e61cc29d5d692fb1d5865bdb7f2e6effbd594756 Mon Sep 17 00:00:00 2001
From: Nicholas Bauer <nicholasbauer@outlook.com>
Date: Fri, 25 Jan 2019 09:09:56 -0500
Subject: [PATCH 02/13] Fixed interface to libRmath

---
 src/distrs/srdist.jl |  4 +--
 src/rmath.jl         | 80 ++++++++++++++++++++++++++++++--------------
 2 files changed, 56 insertions(+), 28 deletions(-)

diff --git a/src/distrs/srdist.jl b/src/distrs/srdist.jl
index 3cf10c9..1eff1a5 100644
--- a/src/distrs/srdist.jl
+++ b/src/distrs/srdist.jl
@@ -25,8 +25,8 @@ function srdistpdf(ν::Real, k::Real, x::Number)
         return inner_part * y^ν * ϕ(y*√ν)
     end
     integral = quadgk(outer, 0.0, Inf)[1]
-    return integral * (√τ * k * (k-1) * ν^(ν/2)) / (gamma(ν/2) * 2^(ν/2 - 1))
+    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(pdf(ν, k, x))
+srdistlogpdf(ν::Real, k::Real, x::Number) = log(srdistpdf(ν, k, x))
diff --git a/src/rmath.jl b/src/rmath.jl
index 9de5982..ca9f922 100644
--- a/src/rmath.jl
+++ b/src/rmath.jl
@@ -34,9 +34,9 @@ function _import_rmath(rname::Symbol, jname::Symbol, pargs)
     invlogcdf = Symbol(jname, "invlogcdf")
     invlogccdf = Symbol(jname, "invlogccdf")
 
-    has_density = true
+    is_tukey = false
     if rname == :tukey
-        has_density = false
+        is_tukey = true
     end
 
     rand = Symbol(jname, "rand")
@@ -51,47 +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
-        if $has_density
+    # 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)
-        end
 
-        $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
@@ -113,7 +141,7 @@ end
 @import_rmath nbinom nbinom r p
 @import_rmath pois pois λ
 @import_rmath t tdist k
-@import_rmath tukey srdist k ν
+@import_rmath tukey srdist ν k
 
 @import_rmath nbeta nbeta α β λ
 @import_rmath nchisq nchisq k λ

From 7ff830456d7df176d05f808fab6096ef00487118 Mon Sep 17 00:00:00 2001
From: Nicholas Bauer <nicholasbauer@outlook.com>
Date: Fri, 25 Jan 2019 12:53:21 -0500
Subject: [PATCH 03/13] Partial tests, but need to better understand what
 they're doing

---
 test/generic.jl |  8 ++++++++
 test/rmath.jl   | 38 +++++++++++++++++++++++++++++---------
 2 files changed, 37 insertions(+), 9 deletions(-)

diff --git a/test/generic.jl b/test/generic.jl
index cd5f717..ee5c914 100644
--- a/test/generic.jl
+++ b/test/generic.jl
@@ -125,3 +125,11 @@ genericcomp_tests("tdist", [
     ((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))
+])
diff --git a/test/rmath.jl b/test/rmath.jl
index 4926213..bd7cdf2 100644
--- a/test/rmath.jl
+++ b/test/rmath.jl
@@ -15,6 +15,12 @@ get_statsfun(fname) = eval(Symbol(fname))
 get_rmathfun(fname) = eval(Meta.parse(string("RFunctions.", fname)))
 
 function rmathcomp(basename, params, X::AbstractArray, rtol=100eps(float(one(eltype(X)))))
+    # tackle pdf specially
+    has_rmathpdf = true
+    if basename == "srdist"
+        has_rmathpdf = false
+    end
+
     pdf     = string(basename, "pdf")
     logpdf  = string(basename, "logpdf")
     cdf     = string(basename, "cdf")
@@ -27,8 +33,10 @@ function rmathcomp(basename, params, X::AbstractArray, rtol=100eps(float(one(elt
     invlogccdf = string(basename, "invlogccdf")
     rand       = string(basename, "rand")
 
-    stats_pdf     = get_statsfun(pdf)
-    stats_logpdf  = get_statsfun(logpdf)
+    if has_rmathpdf
+        stats_pdf     = get_statsfun(pdf)
+        stats_logpdf  = get_statsfun(logpdf)
+    end
     stats_cdf     = get_statsfun(cdf)
     stats_ccdf    = get_statsfun(ccdf)
     stats_logcdf  = get_statsfun(logcdf)
@@ -38,8 +46,10 @@ function rmathcomp(basename, params, X::AbstractArray, rtol=100eps(float(one(elt
     stats_invlogcdf  = get_statsfun(invlogcdf)
     stats_invlogccdf = get_statsfun(invlogccdf)
 
-    rmath_pdf     = get_rmathfun(pdf)
-    rmath_logpdf  = get_rmathfun(logpdf)
+    if has_rmathpdf
+        rmath_pdf     = get_rmathfun(pdf)
+        rmath_logpdf  = get_rmathfun(logpdf)
+    end
     rmath_cdf     = get_rmathfun(cdf)
     rmath_ccdf    = get_rmathfun(ccdf)
     rmath_logcdf  = get_rmathfun(logcdf)
@@ -51,16 +61,18 @@ function rmathcomp(basename, params, X::AbstractArray, rtol=100eps(float(one(elt
 
     # tackle rand specially
     has_rand = true
-    if basename == "nbeta" || basename == "nfdist" || basename == "ntdist"
+    if basename == "nbeta" || basename == "nfdist" || basename == "ntdist" || basename == "srdist"
         has_rand = false
     end
     rmath_rand = has_rand ? get_rmathfun(rand) : nothing
 
     for x in X
-        check_rmath(pdf, stats_pdf, rmath_pdf,
-            params, "x", x, true, rtol)
-        check_rmath(logpdf, stats_logpdf, rmath_logpdf,
-            params, "x", x, false, rtol)
+        if has_rmathpdf
+            check_rmath(pdf, stats_pdf, rmath_pdf,
+                params, "x", x, true, rtol)
+            check_rmath(logpdf, stats_logpdf, rmath_logpdf,
+                params, "x", x, false, rtol)
+        end
         check_rmath(cdf, stats_cdf, rmath_cdf,
             params, "x", x, true, rtol)
         check_rmath(ccdf, stats_ccdf, rmath_ccdf,
@@ -209,3 +221,11 @@ rmathcomp_tests("tdist", [
     ((5,), -5.0:0.1:5.0),
     ((5,), -5:5),
 ])
+
+rmathcomp_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)),
+])

From 85e17185ef408a923857041d1f0ba0afbea03b64 Mon Sep 17 00:00:00 2001
From: Nicholas Bauer <nicholasbauer@outlook.com>
Date: Fri, 25 Jan 2019 18:31:28 -0500
Subject: [PATCH 04/13] Wrapped in TestSet, documented case where functions
 diverge

---
 test/rmath.jl | 244 ++++++++++++++++++++++++++------------------------
 1 file changed, 127 insertions(+), 117 deletions(-)

diff --git a/test/rmath.jl b/test/rmath.jl
index bd7cdf2..c552e92 100644
--- a/test/rmath.jl
+++ b/test/rmath.jl
@@ -112,120 +112,130 @@ end
 
 ### Test cases
 
-rmathcomp_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),
-    ((10, 2), 0.01:0.01:0.99),
-])
-
-rmathcomp_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),
-])
-
-rmathcomp_tests("chisq", [
-    ((1,), 0.0:0.1:8.0),
-    ((4,), 0.0:0.1:8.0),
-    ((9,), 0.0:0.1:8.0),
-    ((9,), 0:8),
-])
-
-rmathcomp_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)),
-    ((10, 3), (0:5)),
-])
-
-rmathcomp_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)),
-    ((2, 3), (1:12)),
-])
-
-rmathcomp_tests("hyper", [
-    ((2, 3, 4), 0.0:4.0),
-    ((2, 3, 4), 0:4)
-])
-
-rmathcomp_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),
-    ((3, 4, 2), 0.01:0.01:0.99),
-])
-
-rmathcomp_tests("nbinom", [
-    ((1, 0.5), 0.0:20.0),
-    ((3, 0.5), 0.0:20.0),
-    ((3, 0.2), 0.0:20.0),
-    ((3, 0.8), 0.0:20.0),
-    ((3, 0.8), 0:20),
-])
-
-rmathcomp_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),
-    ((4, 3), 0:8),
-])
-
-rmathcomp_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),
-    ((2, 3, 1), 1:10),
-])
-
-rmathcomp_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),
-    ((0, 1), -3:3),
-    ((0.0, 0.0), -6.0:0.1:6.0),
-])
-
-rmathcomp_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),
-    ((2, 4), -4:10),
-])
-
-rmathcomp_tests("pois", [
-    ((0.5,), 0.0:20.0),
-    ((1.0,), 0.0:20.0),
-    ((2.0,), 0.0:20.0),
-    ((10.0,), 0.0:20.0),
-    ((10,), 0:20),
-])
-
-rmathcomp_tests("tdist", [
-    ((1,), -5.0:0.1:5.0),
-    ((2,), -5.0:0.1:5.0),
-    ((5,), -5.0:0.1:5.0),
-    ((5,), -5:5),
-])
-
-rmathcomp_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)),
-])
+@testset "Distributions" begin
+    rmathcomp_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),
+        ((10, 2), 0.01:0.01:0.99),
+    ])
+
+    rmathcomp_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),
+    ])
+
+    rmathcomp_tests("chisq", [
+        ((1,), 0.0:0.1:8.0),
+        ((4,), 0.0:0.1:8.0),
+        ((9,), 0.0:0.1:8.0),
+        ((9,), 0:8),
+    ])
+
+    rmathcomp_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)),
+        ((10, 3), (0:5)),
+    ])
+
+    rmathcomp_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)),
+        ((2, 3), (1:12)),
+    ])
+
+    rmathcomp_tests("hyper", [
+        ((2, 3, 4), 0.0:4.0),
+        ((2, 3, 4), 0:4)
+    ])
+
+    rmathcomp_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),
+        ((3, 4, 2), 0.01:0.01:0.99),
+    ])
+
+    rmathcomp_tests("nbinom", [
+        ((1, 0.5), 0.0:20.0),
+        ((3, 0.5), 0.0:20.0),
+        ((3, 0.2), 0.0:20.0),
+        ((3, 0.8), 0.0:20.0),
+        ((3, 0.8), 0:20),
+    ])
+
+    rmathcomp_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),
+        ((4, 3), 0:8),
+    ])
+
+    rmathcomp_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),
+        ((2, 3, 1), 1:10),
+    ])
+
+    rmathcomp_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),
+        ((0, 1), -3:3),
+        ((0.0, 0.0), -6.0:0.1:6.0),
+    ])
+
+    rmathcomp_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),
+        ((2, 4), -4:10),
+    ])
+
+    rmathcomp_tests("pois", [
+        ((0.5,), 0.0:20.0),
+        ((1.0,), 0.0:20.0),
+        ((2.0,), 0.0:20.0),
+        ((10.0,), 0.0:20.0),
+        ((10,), 0:20),
+    ])
+
+    rmathcomp_tests("tdist", [
+        ((1,), -5.0:0.1:5.0),
+        ((2,), -5.0:0.1:5.0),
+        ((5,), -5.0:0.1:5.0),
+        ((5,), -5:5),
+    ])
+
+    rmathcomp_tests("srdist", [
+        ((1,2), (0.0:0.2:5.0)),
+        ((2,2), (0.0:0.2:5.0)),
+        ((5,3), (0.0:0.2:5.0)),
+        ((10,2), (0.0:0.2:5.0)),
+        ((10,5), (0.0:0.2:5.0))
+    ])
+    # Note: Convergence fails in srdist with cdf values below 0.16 with df = 10, k = 5.
+    # Reduced df or k allows convergence. This test documents this behavior.
+    x = 0.15
+    q = srdistcdf(10, 5, 0.15)
+    rx = srdistinvcdf(10, 5, q)
+    rtol = 100eps(1.0)
+    @test_broken x ≈ rx atol=rtol rtol=rtol nans=true
+    ])
+end

From 1e8e3af7c38d75af2efeaf9a9faaccefa29ba489 Mon Sep 17 00:00:00 2001
From: Nicholas Bauer <nicholasbauer@outlook.com>
Date: Sat, 26 Jan 2019 00:01:37 -0500
Subject: [PATCH 05/13] Removed stray bracket

---
 test/rmath.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/rmath.jl b/test/rmath.jl
index c552e92..2353a43 100644
--- a/test/rmath.jl
+++ b/test/rmath.jl
@@ -230,6 +230,7 @@ end
         ((10,2), (0.0:0.2:5.0)),
         ((10,5), (0.0:0.2:5.0))
     ])
+
     # Note: Convergence fails in srdist with cdf values below 0.16 with df = 10, k = 5.
     # Reduced df or k allows convergence. This test documents this behavior.
     x = 0.15
@@ -237,5 +238,4 @@ end
     rx = srdistinvcdf(10, 5, q)
     rtol = 100eps(1.0)
     @test_broken x ≈ rx atol=rtol rtol=rtol nans=true
-    ])
 end

From 694dbdedc8dcfcb55d0652026437b5d1af5860cd Mon Sep 17 00:00:00 2001
From: Nicholas Bauer <nicholasbauer@outlook.com>
Date: Sun, 27 Jan 2019 11:27:55 -0500
Subject: [PATCH 06/13] Corrected pdf

---
 src/distrs/srdist.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/distrs/srdist.jl b/src/distrs/srdist.jl
index 1eff1a5..b127be0 100644
--- a/src/distrs/srdist.jl
+++ b/src/distrs/srdist.jl
@@ -25,7 +25,7 @@ function srdistpdf(ν::Real, k::Real, x::Number)
         return inner_part * y^ν * ϕ(y*√ν)
     end
     integral = quadgk(outer, 0.0, Inf)[1]
-    return integral * (√2π * k * (k-1) * ν^(ν/2)) / (gamma(ν/2) * 2^(ν/2 - 1))
+    return integral * (√(2π) * k * (k-1) * ν^(ν/2)) / (gamma(ν/2) * 2^(ν/2 - 1))
 end
 
 # logpdf for numbers with generic types

From 2487075478b5631bf2492ca9afbb4214339315da Mon Sep 17 00:00:00 2001
From: Nicholas Bauer <nicholasbauer@outlook.com>
Date: Sun, 27 Jan 2019 11:45:44 -0500
Subject: [PATCH 07/13] Simplified pdf code

---
 src/distrs/srdist.jl | 8 +++-----
 1 file changed, 3 insertions(+), 5 deletions(-)

diff --git a/src/distrs/srdist.jl b/src/distrs/srdist.jl
index b127be0..be7ef2f 100644
--- a/src/distrs/srdist.jl
+++ b/src/distrs/srdist.jl
@@ -18,13 +18,11 @@ function srdistpdf(ν::Real, k::Real, x::Number)
     Φ(x) = (1+erf(x / √2)) / 2
     ϕ(x) = derivative(Φ, x)
     function outer(y)
-        function inner(u)
-            return ϕ(u) * ϕ(u - x*y) * (Φ(u) - Φ(u - x*y))^(k-2)
-        end
-        inner_part = quadgk(inner, -Inf, Inf)[1]
+        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)[1]
+    integral = quadgk(outer, 0.0, Inf) |> first
     return integral * (√(2π) * k * (k-1) * ν^(ν/2)) / (gamma(ν/2) * 2^(ν/2 - 1))
 end
 

From 5e0d949cd2fa7c32666760e5c4bb960008d4284d Mon Sep 17 00:00:00 2001
From: Nicholas Bauer <nicholasbauer@outlook.com>
Date: Mon, 28 Jan 2019 21:31:48 -0500
Subject: [PATCH 08/13] Switched srdistpdf internal functions to use
 normpdf/cdf

Vast improvement in performance.
---
 REQUIRE              | 1 -
 src/distrs/srdist.jl | 4 ++--
 2 files changed, 2 insertions(+), 3 deletions(-)

diff --git a/REQUIRE b/REQUIRE
index aa33dd5..0d480d0 100644
--- a/REQUIRE
+++ b/REQUIRE
@@ -2,4 +2,3 @@ julia 0.7-beta
 Rmath 0.4.0
 SpecialFunctions
 QuadGK
-ForwardDiff
diff --git a/src/distrs/srdist.jl b/src/distrs/srdist.jl
index be7ef2f..84a7d0b 100644
--- a/src/distrs/srdist.jl
+++ b/src/distrs/srdist.jl
@@ -15,8 +15,8 @@ import .RFunctions:
 
 # pdf for numbers with generic types
 function srdistpdf(ν::Real, k::Real, x::Number)
-    Φ(x) = (1+erf(x / √2)) / 2
-    ϕ(x) = derivative(Φ, x)
+    Φ(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

From dbc1d6332311c3bb79668a143c5f2f9507e84d00 Mon Sep 17 00:00:00 2001
From: Nicholas Bauer <nicholasbauer@outlook.com>
Date: Mon, 28 Jan 2019 21:33:31 -0500
Subject: [PATCH 09/13] Removed ForwardDiff

---
 src/distrs/srdist.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/distrs/srdist.jl b/src/distrs/srdist.jl
index 84a7d0b..9e1e77a 100644
--- a/src/distrs/srdist.jl
+++ b/src/distrs/srdist.jl
@@ -1,6 +1,5 @@
 # functions related to studentized range distribution
 
-using ForwardDiff: derivative
 using QuadGK: quadgk
 
 import .RFunctions:

From 34a4b92efb11acbf20de3e0202c5b07379c22c63 Mon Sep 17 00:00:00 2001
From: Nicholas Bauer <nicholasbauer@outlook.com>
Date: Mon, 28 Jan 2019 21:38:46 -0500
Subject: [PATCH 10/13] Wrapped generic tests in TestSet too

---
 test/generic.jl | 182 ++++++++++++++++++++++++------------------------
 test/rmath.jl   |   2 +-
 2 files changed, 93 insertions(+), 91 deletions(-)

diff --git a/test/generic.jl b/test/generic.jl
index ee5c914..ccb07df 100644
--- a/test/generic.jl
+++ b/test/generic.jl
@@ -43,93 +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),
-])
-
-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))
-])
+@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
diff --git a/test/rmath.jl b/test/rmath.jl
index 2353a43..1a5e5d1 100644
--- a/test/rmath.jl
+++ b/test/rmath.jl
@@ -112,7 +112,7 @@ end
 
 ### Test cases
 
-@testset "Distributions" begin
+@testset "RMath" begin
     rmathcomp_tests("beta", [
         ((1.0, 1.0), 0.01:0.01:0.99),
         ((2.0, 3.0), 0.01:0.01:0.99),

From 087c0aa4a905af0ff2e371f15fd73b56a9c7af02 Mon Sep 17 00:00:00 2001
From: Nicholas Bauer <nicholasbauer@outlook.com>
Date: Tue, 29 Jan 2019 14:08:59 -0500
Subject: [PATCH 11/13] Removed pdf and associated tests

---
 src/StatsFuns.jl     |  2 --
 src/distrs/srdist.jl | 18 ------------------
 test/generic.jl      | 14 +++++++-------
 test/rmath.jl        | 16 +++++++++-------
 4 files changed, 16 insertions(+), 34 deletions(-)

diff --git a/src/StatsFuns.jl b/src/StatsFuns.jl
index f5a44cb..3e51787 100644
--- a/src/StatsFuns.jl
+++ b/src/StatsFuns.jl
@@ -218,8 +218,6 @@ export
     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
diff --git a/src/distrs/srdist.jl b/src/distrs/srdist.jl
index 9e1e77a..8d1e180 100644
--- a/src/distrs/srdist.jl
+++ b/src/distrs/srdist.jl
@@ -1,7 +1,5 @@
 # functions related to studentized range distribution
 
-using QuadGK: quadgk
-
 import .RFunctions:
     srdistcdf,
     srdistccdf,
@@ -11,19 +9,3 @@ import .RFunctions:
     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))
diff --git a/test/generic.jl b/test/generic.jl
index ccb07df..1b4bf0e 100644
--- a/test/generic.jl
+++ b/test/generic.jl
@@ -127,11 +127,11 @@ end
         ((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))
-    ])
+    #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
diff --git a/test/rmath.jl b/test/rmath.jl
index 1a5e5d1..138cce8 100644
--- a/test/rmath.jl
+++ b/test/rmath.jl
@@ -16,13 +16,15 @@ get_rmathfun(fname) = eval(Meta.parse(string("RFunctions.", fname)))
 
 function rmathcomp(basename, params, X::AbstractArray, rtol=100eps(float(one(eltype(X)))))
     # tackle pdf specially
-    has_rmathpdf = true
+    has_pdf = true
     if basename == "srdist"
-        has_rmathpdf = false
+        has_pdf = false
     end
 
-    pdf     = string(basename, "pdf")
-    logpdf  = string(basename, "logpdf")
+    if has_pdf
+        pdf     = string(basename, "pdf")
+        logpdf  = string(basename, "logpdf")
+    end
     cdf     = string(basename, "cdf")
     ccdf    = string(basename, "ccdf")
     logcdf  = string(basename, "logcdf")
@@ -33,7 +35,7 @@ function rmathcomp(basename, params, X::AbstractArray, rtol=100eps(float(one(elt
     invlogccdf = string(basename, "invlogccdf")
     rand       = string(basename, "rand")
 
-    if has_rmathpdf
+    if has_pdf
         stats_pdf     = get_statsfun(pdf)
         stats_logpdf  = get_statsfun(logpdf)
     end
@@ -46,7 +48,7 @@ function rmathcomp(basename, params, X::AbstractArray, rtol=100eps(float(one(elt
     stats_invlogcdf  = get_statsfun(invlogcdf)
     stats_invlogccdf = get_statsfun(invlogccdf)
 
-    if has_rmathpdf
+    if has_pdf
         rmath_pdf     = get_rmathfun(pdf)
         rmath_logpdf  = get_rmathfun(logpdf)
     end
@@ -67,7 +69,7 @@ function rmathcomp(basename, params, X::AbstractArray, rtol=100eps(float(one(elt
     rmath_rand = has_rand ? get_rmathfun(rand) : nothing
 
     for x in X
-        if has_rmathpdf
+        if has_pdf
             check_rmath(pdf, stats_pdf, rmath_pdf,
                 params, "x", x, true, rtol)
             check_rmath(logpdf, stats_logpdf, rmath_logpdf,

From bc3758164e1d524e331eb5bb0ab3f53171c7162d Mon Sep 17 00:00:00 2001
From: Nicholas Bauer <nicholasbauer@outlook.com>
Date: Tue, 29 Jan 2019 14:10:01 -0500
Subject: [PATCH 12/13] Removed QuadGK requirement

---
 REQUIRE | 1 -
 1 file changed, 1 deletion(-)

diff --git a/REQUIRE b/REQUIRE
index 0d480d0..02ca4f6 100644
--- a/REQUIRE
+++ b/REQUIRE
@@ -1,4 +1,3 @@
 julia 0.7-beta
 Rmath 0.4.0
 SpecialFunctions
-QuadGK

From 259bcda5a1fb9c580a84818224614d1dca16f66c Mon Sep 17 00:00:00 2001
From: Nicholas Bauer <nicholasbauer@outlook.com>
Date: Tue, 29 Jan 2019 15:56:32 -0500
Subject: [PATCH 13/13] Consolidated tukey cases with others

---
 src/rmath.jl | 99 +++++++++++++++++++++-------------------------------
 1 file changed, 40 insertions(+), 59 deletions(-)

diff --git a/src/rmath.jl b/src/rmath.jl
index ca9f922..e2b4c77 100644
--- a/src/rmath.jl
+++ b/src/rmath.jl
@@ -34,92 +34,73 @@ 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
+    is_tukey = rname == :tukey
 
     rand = Symbol(jname, "rand")
     has_rand = true
-    if rname == :nbeta || rname == :nf || rname == :nt || rname == :tukey
+    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)
-    tukeytypes = Expr(:tuple, :Cdouble, :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
-        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)
+        # ptukey and qtukey have an extra literal 1 argument
+        pargs = (1, pargs...)
+    end
 
-            $invlogccdf($(pdecls...), lq::Union{Float64,Int}) =
-                ccall(($qfun, libRmath), Float64, $tukeytypes, lq, 1, $(pargs...), 0, 1)
-        end
-    else
-        quote
+    # Function implementation
+    quote
+        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)
+        end
 
-            $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...))
-            end
+        if $has_rand
+            $rand($(pdecls...)) =
+                ccall(($rfun, libRmath), Float64, $rtypes, $(pargs...))
         end
     end
 end
@@ -141,7 +122,7 @@ end
 @import_rmath nbinom nbinom r p
 @import_rmath pois pois λ
 @import_rmath t tdist k
-@import_rmath tukey srdist ν k
+@import_rmath tukey srdist k ν
 
 @import_rmath nbeta nbeta α β λ
 @import_rmath nchisq nchisq k λ