twiddle with rounding of signrank quantile

src/distrs/signrank.jl

@@ -30,44 +30,39 @@ function signrankpdf(n::Int, W::Union{Float64,Int})
     max_W = (n * (n + 1)) >> 1
     if W < 0
         return 0.0
-    elseif W >= max_W + 1
-        return 0.0
     elseif W > max_W >> 1
         return signrankpdf(n, max_W - W)
     end
     DP = signrankDP(n, W)
     return ldexp(float(DP[1]), -n)
 end
 
-function signranklogpdf(n::Int, W::Union{Float64,Int}, )
+function signranklogpdf(n::Int, W::Union{Float64,Int})
     return log(signrankpdf(n, W))
 end
 
-function signrankcdf(n::Int, W::Union{Float64,Int}, )
+function signrankcdf(n::Int, W::Union{Float64,Int})
     W = round(Int, W)
     max_W = (n * (n + 1)) >> 1
     if W < 0
         return 0.0
-    elseif W >= max_W
-        return 1.0
     elseif W > max_W >> 1
-        return 1.0 - signrankcdf(n, max_W - W - 1, )
+        return 1.0 - signrankcdf(n, max_W - W - 1)
     end
     DP = signrankDP(n, W)
     return sum(Base.Fix2(ldexp, -n) ∘ float, DP)
 end
 
-function signranklogcdf(n::Int, W::Union{Float64,Int}, )
+function signranklogcdf(n::Int, W::Union{Float64,Int})
     return log(signrankcdf(n, W))
 end
 
-function signrankccdf(n::Int, W::Union{Float64,Int}, )
-    W = round(Int, W)
+function signrankccdf(n::Int, W::Union{Float64,Int})
     max_W = (n * (n + 1)) >> 1
     return signrankcdf(n, max_W - W - 1, )
 end
 
-function signranklogccdf(n::Int, W::Union{Float64,Int},)
+function signranklogccdf(n::Int, W::Union{Float64,Int})
     return log(signrankccdf(n, W))
 end
 
@@ -83,13 +78,20 @@ function signrankinvcdf(n::Int, p::Float64)
 end
 
 function signrankinvlogcdf(n::Int, logp::Float64)
-    signrankinvcdf(n, exp(logp))
+    if logp > 0.0
+        return NaN
+    end
+    W = 0
+    while signranklogcdf(n, W) < logp # TODO binary search and symmetry
+        W += 1
+    end
+    return float(W)
 end
 
 function signrankinvccdf(n::Int, p::Float64)
     signrankinvcdf(n, 1 - p)
 end
 
 function signrankinvlogccdf(n::Int, logp::Float64)
-    signrankinvccdf(n, exp(logp))
+    signrankinvlogcdf(n, log1mexp(logp))
 end

test/rmath.jl

@@ -218,7 +218,8 @@ end
 
     rmathcomp_tests("signrank", [
         ((4),-2:12),
-        #((10),-2:57),
+        ((10),-2:57),
+#=         ((50),-2:1277), =#
     ])
 
     rmathcomp_tests("binom", [