Fix tests (#157)

* Make `generic.jl` tests more similar to `rmath.jl`

* Try to tune tolerances

* Some fixes

* Refactor tests

* Tune tolerances

Author: devmotion

test/generic.jl

@@ -1,38 +1,35 @@
 using StatsFuns
 using StatsFuns: RFunctions
 using ForwardDiff: Dual
+using Test
+
+include("utils.jl")
 
 function check_rmath(fname, statsfun, rmathfun, params, aname, a, isprob, rtol)
     v = @inferred(rmathfun(params..., a))
     rv = @inferred(statsfun(params..., Dual(a))).value
+    @test v isa float(Base.promote_typeof(params..., a))
+    @test rv isa float(Base.promote_typeof(params..., a))
     if isprob
-        rd = abs(v / rv - 1.0)
-        if rd > rtol
-            error("$fname deviates too much from Rmath at " *
-                "params = $params, $aname = $a:\n" *
-                "  v = $v (rv = $rv)\n  |v/rv - 1| = $rd > $rtol.")
-        end
+        @test v ≈ rv rtol=rtol nans=true
     else
-        τ = (1.0 + abs(rv)) * rtol
-        ad = abs(v - rv)
-        if ad > τ
-            error("$fname deviates too much from Rmath at " *
-                "params = $params, $aname = $a:\n" *
-                "  v = $v (rv = $rv)\n  |v - rv| = $ad > $τ.")
-        end
+        @test v ≈ rv atol=rtol rtol=rtol nans=true
     end
 end
 
-function genericcomp(basename, params, X::AbstractArray, rtol=100eps(float(one(eltype(X)))))
+function genericcomp(basename, params, X::AbstractArray, rtol=_default_rtol(params, X))
   pdf = string(basename, "pdf")
   logpdf = string(basename, "logpdf")
   stats_pdf = eval(Symbol(pdf))
   stats_logpdf = eval(Symbol(logpdf))
   rmath_pdf = eval(Meta.parse(string("RFunctions.", pdf)))
   rmath_logpdf = eval(Meta.parse(string("RFunctions.", logpdf)))
-  for i = 1:length(X)
-    x = X[i]
+
+  @testset "pdf with params=$params, x=$x" for x in X
     check_rmath(pdf, stats_pdf, rmath_pdf, params, "x", x, true, rtol)
+  end
+
+  @testset "logpdf with params=$params, x=$x" for x in X
     check_rmath(logpdf, stats_logpdf, rmath_logpdf, params, "x", x, false, rtol)
   end
 end
@@ -121,8 +118,11 @@ end
 
     genericcomp_tests("pois", [
         ((1.0,), 0:30),
-        ((10.0,), 0:42)
+        ((10.0,), 0:37),
+        ((10.0,), 39:42),
     ])
+    # Requires slightly larger tolerance: #157
+    genericcomp("pois", (10.0,), 38:38, 2.5e-14)
 
     genericcomp_tests("tdist", [
         ((1,), -5.0:0.1:5.0),

test/misc.jl

@@ -1,4 +1,5 @@
 using SpecialFunctions, StatsFuns
+using Test
 
 @testset "logmvgamma" begin
     @testset "type behavior" for eltya in (Float32, Float64)

test/rmath.jl

@@ -2,6 +2,8 @@ using StatsFuns
 using StatsFuns: RFunctions
 using Test
 
+include("utils.jl")
+
 function check_rmath(fname, statsfun, rmathfun, params, aname, a, isprob, rtol)
     v = @inferred(statsfun(params..., a))
     rv = @inferred(rmathfun(params..., a))
@@ -17,15 +19,7 @@ end
 get_statsfun(fname) = eval(Symbol(fname))
 get_rmathfun(fname) = eval(Meta.parse(string("RFunctions.", fname)))
 
-function rmathcomp(basename, params, X::AbstractArray)
-    # compute default tolerance:
-    # has to take into account `params` as well since otherwise e.g. `X::Array{<:Rational}`
-    # always uses a tolerance based on `eps(one(Float64))` even when parameters are of type
-    # Float32
-    rtol = 100 * eps(float(one(promote_type(Base.promote_typeof(params...), eltype(X)))))
-    rmathcomp(basename, params, X, rtol)
-end
-function rmathcomp(basename, params, X::AbstractArray, rtol)
+function rmathcomp(basename, params, X::AbstractArray, rtol=_default_rtol(params, X))
     # tackle pdf specially
     has_pdf = true
     if basename == "srdist"

test/utils.jl

@@ -0,0 +1,14 @@
+# default relative tolerance for comparisons with Rmath
+function _default_rtol(params, X::AbstractArray)
+    # has to take into account `params` as well since otherwise e.g. `X::Array{<:Rational}`
+    # always uses a tolerance based on `eps(one(Float64))` even when parameters are of type
+    # Float32
+    return _default_rtol(float(promote_type(Base.promote_typeof(params...), eltype(X))))
+end
+
+# We use less sharp tolerances for Float16 and Float32 since the proportion of significant
+# digits that are equal when evaluating Rmath and StatsFuns is smaller
+# eps^(7//8) means requiring equality of about 7/8 of the significant digits, etc.
+_default_rtol(::Type{Float64}) = eps(Float64)^(7//8) # ~2.0e-14
+_default_rtol(::Type{Float32}) = eps(Float32)^(3//4) # ~6.4e-6
+_default_rtol(::Type{Float16}) = eps(Float16)^(2//3) # ~9.9e-3