Fix inconsistent return types of quantile etc.

src/univariate/discrete/bernoulli.jl

@@ -96,11 +96,13 @@ ccdf(d::Bernoulli, x::Bool) = x ? zero(d.p) : succprob(d)
 ccdf(d::Bernoulli, x::Int) = x < 0 ? one(d.p) :
                              x < 1 ? succprob(d) : zero(d.p)
 
-function quantile(d::Bernoulli{T}, p::Real) where T<:Real
-    0 <= p <= 1 ? (p <= failprob(d) ? zero(T) : one(T)) : T(NaN)
+function quantile(d::Bernoulli, p::Real)
+    _check_quantile_arg(p)
+    p <= failprob(d) ? false : true
 end
-function cquantile(d::Bernoulli{T}, p::Real) where T<:Real
-    0 <= p <= 1 ? (p >= succprob(d) ? zero(T) : one(T)) : T(NaN)
+function cquantile(d::Bernoulli, p::Real)
+    _check_cquantile_arg(p)
+    p >= succprob(d) ? false : true
 end
 
 mgf(d::Bernoulli, t::Real) = failprob(d) + succprob(d) * exp(t)

src/univariate/discrete/bernoullilogit.jl

@@ -87,12 +87,12 @@ logccdf(d::BernoulliLogit, x::Bool) = x ? oftype(float(d.logitp), -Inf) : logsuc
 logccdf(d::BernoulliLogit, x::Int) = x < 0 ? zero(float(d.logitp)) : (x < 1 ? logsuccprob(d) : oftype(float(d.logitp), -Inf))
 
 function quantile(d::BernoulliLogit, p::Real)
-    T = float(partype(d))
-    0 <= p <= 1 ? (p <= failprob(d) ? zero(T) : one(T)) : T(NaN)
+    _check_quantile_arg(p)
+    p <= failprob(d) ? false : true
 end
 function cquantile(d::BernoulliLogit, p::Real)
-    T = float(partype(d))
-    0 <= p <= 1 ? (p >= succprob(d) ? zero(T) : one(T)) : T(NaN)
+    _check_cquantile_arg(q)
+    p >= succprob(d) ? false : true
 end
 
 mgf(d::BernoulliLogit, t::Real) = failprob(d) + exp(t + logsuccprob(d))

src/univariate/discrete/binomial.jl

@@ -126,7 +126,60 @@ end
 
 #### Evaluation & Sampling
 
-@_delegate_statsfuns Binomial binom n p
+# We rely on Rmath for (log)pdf, (log)cdf, and quantile functions
+# We do not use the `@_delegate_statsfuns` macro to work around some issues in the 
+# quantile functions such as `StatsFuns.binominvcdf(10, 0.0, 1.0) = 10`.
+pdf(d::Binomial, x::Real) = binompdf(d.n, d.p, x)
+logpdf(d::Binomial, x::Real) = binomlogpdf(d.n, d.p, x)
+cdf(d::Binomial, x::Real) = binomcdf(d.n, d.p, x)
+ccdf(d::Binomial, x::Real) = binomccdf(d.n, d.p, x)
+logcdf(d::Binomial, x::Real) = binomlogcdf(d.n, d.p, x)
+logccdf(d::Binomial, x::Real) = binomlogccdf(d.n, d.p, x)
+
+function quantile(d::Binomial, q::Real)::Int
+    _check_quantile_arg(q)
+    if iszero(d.p)
+        return 0
+    elseif isone(d.p)
+        return d.n
+    else
+        return binominvcdf(d.n, d.p, q)
+    end
+end
+function cquantile(d::Binomial, q::Real)::Int
+    _check_cquantile_arg(q)
+    if iszero(d.p)
+        return 0
+    elseif isone(d.p)
+        return d.n
+    else
+        return binominvccdf(d.n, d.p, q)
+    end
+end
+function invlogcdf(d::Binomial, lq::Real)::Int
+    _check_invlogcdf_arg(lq)
+    if iszero(d.p)
+        return 0
+    elseif isone(d.p)
+        return d.n
+    elseif isinf(lq)
+        return 0
+    else
+        return binominvlogcdf(d.n, d.p, lq)
+    end
+end
+function invlogccdf(d::Binomial, lq::Real)::Int
+    _check_invlogccdf_arg(lq)
+    if iszero(d.p)
+        return 0
+    elseif isone(d.p)
+        return d.n
+    elseif isinf(lq)
+        return d.n
+    else
+        return binominvlogccdf(d.n, d.p, lq)
+    end
+end
 
 function rand(rng::AbstractRNG, d::Binomial)
     p, n = d.p, d.n

src/univariate/discrete/dirac.jl

@@ -47,7 +47,10 @@ logcdf(d::Dirac, x::Real) = x < d.value ? -Inf : isnan(x) ? NaN : 0.0
 ccdf(d::Dirac, x::Real) = x < d.value ? 1.0 : isnan(x) ? NaN : 0.0
 logccdf(d::Dirac, x::Real) = x < d.value ? 0.0 : isnan(x) ? NaN : -Inf
 
-quantile(d::Dirac{T}, p::Real) where {T} = 0 <= p <= 1 ? d.value : T(NaN)
+function quantile(d::Dirac, p::Real)
+    _check_quantile_arg(p)
+    return d.value
+end
 
 mgf(d::Dirac, t) = exp(t * d.value)
 cgf(d::Dirac, t) = t*d.value

src/univariate/discrete/discreteuniform.jl

@@ -85,7 +85,10 @@ function cdf(d::DiscreteUniform, x::Int)
     end
 end
 
-quantile(d::DiscreteUniform, p::Real) = iszero(p) ? d.a : d.a - 1 + ceil(Int, p * span(d))
+function quantile(d::DiscreteUniform, p::Real)
+    _check_quantile_arg(p)
+    iszero(p) ? d.a : d.a - 1 + ceil(Int, p * span(d))
+end
 
 function mgf(d::DiscreteUniform, t::Real)
     a, b = d.a, d.b

src/univariate/discrete/geometric.jl

@@ -105,21 +105,23 @@ logcdf(d::Geometric, x::Int) = log1mexp(log1p(-d.p) * max(x + 1, 0))
 logccdf(d::Geometric, x::Real) = logccdf_int(d, x)
 logccdf(d::Geometric, x::Int) =  log1p(-d.p) * max(x + 1, 0)
 
-quantile(d::Geometric, p::Real) = invlogccdf(d, log1p(-p))
-
-cquantile(d::Geometric, p::Real) = invlogccdf(d, log(p))
-
-invlogcdf(d::Geometric, lp::Real) = invlogccdf(d, log1mexp(lp))
+function quantile(d::Geometric, p::Real)
+    _check_quantile_arg(p)
+    return invlogccdf(d, log1p(-p))
+end
+function cquantile(d::Geometric, p::Real)
+    _check_cquantile_arg(p)
+    return invlogccdf(d, log(p))
+end
 
-function invlogccdf(d::Geometric{T}, lp::Real) where T<:Real
-    if (lp > zero(d.p)) || isnan(lp)
-        return T(NaN)
-    elseif isinf(lp)
-        return T(Inf)
-    elseif lp == zero(d.p)
-        return zero(T)
-    end
-    max(ceil(lp/log1p(-d.p)) - 1, zero(T))
+function invlogcdf(d::Geometric, lp::Real)
+    _check_invlogcdf_arg(lp)
+    return invlogccdf(d, log1mexp(lp))
+end
+function invlogccdf(d::Geometric, lp::Real)
+    _check_invlogccdf_arg(lp)
+    z = lp/log1p(-d.p) - 1
+    return ceil(Int, max(z, zero(z)))
 end
 
 function laplace_transform(d::Geometric, t)

src/univariate/discrete/negativebinomial.jl

@@ -112,10 +112,30 @@ cdf(d::NegativeBinomial, x::Real) = nbinomcdf(d.r, d.p, x)
 ccdf(d::NegativeBinomial, x::Real) = nbinomccdf(d.r, d.p, x)
 logcdf(d::NegativeBinomial, x::Real) = nbinomlogcdf(d.r, d.p, x)
 logccdf(d::NegativeBinomial, x::Real) = nbinomlogccdf(d.r, d.p, x)
-quantile(d::NegativeBinomial, q::Real) = convert(Int, nbinominvcdf(d.r, d.p, q))
-cquantile(d::NegativeBinomial, q::Real) = convert(Int, nbinominvccdf(d.r, d.p, q))
-invlogcdf(d::NegativeBinomial, lq::Real) = convert(Int, nbinominvlogcdf(d.r, d.p, lq))
-invlogccdf(d::NegativeBinomial, lq::Real) = convert(Int, nbinominvlogccdf(d.r, d.p, lq))
+function quantile(d::NegativeBinomial, q::Real)::Int
+    _check_quantile_arg(q)
+    return nbinominvcdf(d.r, d.p, q)
+end
+function cquantile(d::NegativeBinomial, q::Real)::Int
+    _check_cquantile_arg(q)
+    return nbinominvccdf(d.r, d.p, q)
+end
+function invlogcdf(d::NegativeBinomial, lq::Real)::Int
+    _check_invlogcdf_arg(lq)
+    if isinf(lq)
+        return nbinominvcdf(d.r, d.p, zero(lq))
+    else
+        return nbinominvlogcdf(d.r, d.p, lq)
+    end
+end
+function invlogccdf(d::NegativeBinomial, lq::Real)::Int
+    _check_invlogccdf_arg(lq)
+    if isinf(lq)
+        return nbinominvccdf(d.r, d.p, zero(lq))
+    else
+        return nbinominvlogccdf(d.r, d.p, lq)
+    end
+end
 
 ## sampling
 function rand(rng::AbstractRNG, d::NegativeBinomial)

src/univariate/discrete/noncentralhypergeometric.jl

@@ -9,23 +9,23 @@ abstract type NoncentralHypergeometric{T<:Real} <: DiscreteUnivariateDistributio
 
 # Functions
 
-function quantile(d::NoncentralHypergeometric{T}, q::Real) where T<:Real
-    if !(zero(q) <= q <= one(q))
-        T(NaN)
-    else
-        range = support(d)
-        if q > 1/2
-            q = 1 - q
-            range = reverse(range)
-        end
+function quantile(d::NoncentralHypergeometric, q::Real)
+    _check_quantile_arg(q)
 
-        qsum, i = zero(T), 0
-        while qsum < q
-            i += 1
-            qsum += pdf(d, range[i])
-        end
-        range[i]
+    range = support(d)
+    if q > 1/2
+        q = 1 - q
+        range = reverse(range)
+    end
+
+    # Support must not be empty
+    r, state = iterate(range)
+    qsum = pdf(d, r)
+    while qsum < q && (r_state = iterate(range, state)) !== nothing
+        r, state = r_state
+        qsum += pdf(d, r)
     end
+    return r
 end
 
 params(d::NoncentralHypergeometric) = (d.ns, d.nf, d.n, d.ω)

src/univariates.jl

@@ -392,21 +392,30 @@ quantile(d::UnivariateDistribution, p::Real)
 
 The complementary quantile value, i.e. `quantile(d, 1-q)`.
 """
-cquantile(d::UnivariateDistribution, p::Real) = quantile(d, 1.0 - p)
+function cquantile(d::UnivariateDistribution, p::Real)
+    _check_cquantile_arg(p)
+    quantile(d, 1 - p)
+end
 
 """
     invlogcdf(d::UnivariateDistribution, lp::Real)
 
 The inverse function of logcdf.
 """
-invlogcdf(d::UnivariateDistribution, lp::Real) = quantile(d, exp(lp))
+function invlogcdf(d::UnivariateDistribution, lp::Real)
+    _check_invlogcdf_arg(lp)
+    quantile(d, exp(lp))
+end
 
 """
     invlogccdf(d::UnivariateDistribution, lp::Real)
 
 The inverse function of logccdf.
 """
-invlogccdf(d::UnivariateDistribution, lp::Real) = quantile(d, -expm1(lp))
+function invlogccdf(d::UnivariateDistribution, lp::Real)
+    _check_invlogccdf_arg(lp)
+    quantile(d, -expm1(lp))
+end
 
 # gradlogpdf
 
@@ -609,6 +618,38 @@ function integerunitrange_logccdf(d::DiscreteUnivariateDistribution, x::Integer)
     return result
 end
 
+### Error messages
+
+const QUANTILE_ARG_ERROR = "`Distributions.quantile(d, p)` only accepts arguments `p` in `[0, 1]`." 
+const CQUANTILE_ARG_ERROR = "`Distributions.cquantile(d, p)` only accepts arguments `p` in `[0, 1]`." 
+const INVLOGCDF_ARG_ERROR = "`Distributions.invlogcdf(d, lp)` only accepts non-negative arguments `lp`." 
+const INVLOGCCDF_ARG_ERROR = "`Distributions.invlogccdf(d, lp)` only accepts non-negative arguments `lp`." 
+
+function _check_quantile_arg(p::Real)
+    if !(zero(p) <= p <= oneunit(p))
+        throw(DomainError(p, QUANTILE_ARG_ERROR))
+    end
+    nothing
+end
+function _check_cquantile_arg(p::Real)
+    if !(zero(p) <= p <= oneunit(p))
+        throw(DomainError(p, CQUANTILE_ARG_ERROR))
+    end
+    nothing
+end
+function _check_invlogcdf_arg(lp::Real)
+    if !(lp <= zero(lp))
+        throw(DomainError(lp, INVLOGCDF_ARG_ERROR))
+    end
+    nothing
+end
+function _check_invlogccdf_arg(lp::Real)
+    if !(lp <= zero(lp))
+        throw(DomainError(lp, INVLOGCDF_ARG_ERROR))
+    end
+    nothing
+end
+
 ### macros to use StatsFuns for method implementation
 
 macro _delegate_statsfuns(D, fpre, psyms...)
@@ -639,10 +680,34 @@ macro _delegate_statsfuns(D, fpre, psyms...)
         $Distributions.ccdf(d::$D, x::Real) = $(fccdf)($(pargs...), x)
         $Distributions.logccdf(d::$D, x::Real) = $(flogccdf)($(pargs...), x)
 
-        $Distributions.quantile(d::$D, q::Real) = convert($T, $(finvcdf)($(pargs...), q))
-        $Distributions.cquantile(d::$D, q::Real) = convert($T, $(finvccdf)($(pargs...), q))
-        $Distributions.invlogcdf(d::$D, lq::Real) = convert($T, $(finvlogcdf)($(pargs...), lq))
-        $Distributions.invlogccdf(d::$D, lq::Real) = convert($T, $(finvlogccdf)($(pargs...), lq))
+        function $Distributions.quantile(d::$D, q::Real)::$T
+            _check_quantile_arg(q)
+            return $(finvcdf)($(pargs...), q)
+        end
+        function $Distributions.cquantile(d::$D, q::Real)::$T
+            _check_cquantile_arg(q)
+            return $(finvccdf)($(pargs...), q)
+        end
+
+        # There is a bug in some functions in StatsFuns.RFunctions/Rmath:
+        # For `lq = -Inf` they return `NaN`
+        # We work around this issue by falling back to quantile/cquantile
+        function $Distributions.invlogcdf(d::$D, lq::Real)::$T
+            _check_invlogcdf_arg(lq)
+            if isinf(lq)
+                return $(finvcdf)($(pargs...), zero(lq))
+            else
+                return $(finvlogcdf)($(pargs...), lq)
+            end
+        end
+        function $Distributions.invlogccdf(d::$D, lq::Real)::$T
+            _check_invlogccdf_arg(lq)
+            if isinf(lq)
+                return $(finvccdf)($(pargs...), zero(lq))
+            else
+                return $(finvlogccdf)($(pargs...), lq)
+            end
+        end
     end
 end
 

test/testutils.jl

@@ -403,6 +403,7 @@ end
 
 
 function test_evaluation(d::DiscreteUnivariateDistribution, vs::AbstractVector, testquan::Bool=true)
+    T = eltype(vs)
     nv  = length(vs)
     p   = Vector{Float64}(undef, nv)
     c   = Vector{Float64}(undef, nv)
@@ -433,11 +434,10 @@ function test_evaluation(d::DiscreteUnivariateDistribution, vs::AbstractVector,
         if testquan
             ep = 1.0e-8
             if p[i] > 2 * ep   # ensure p[i] is large enough to guarantee a reliable result
-                @test quantile(d, c[i] - ep) == v
-                @test cquantile(d, cc[i] + ep) == v
-                @test invlogcdf(d, lc[i] - ep) == v
-                if 0.0 < c[i] < 1.0
-                    @test invlogccdf(d, lcc[i] + ep) == v
+                for (f, z) in ((quantile, c[i] - ep), (cquantile, cc[i] + ep), (invlogcdf, lc[i] - ep), (invlogccdf, lcc[i] + ep))
+                    fz = @inferred(f(d, z))
+                    @test fz isa T
+                    @test fz == v
                 end
             end
         end
@@ -455,6 +455,7 @@ end
 
 
 function test_evaluation(d::ContinuousUnivariateDistribution, vs::AbstractVector, testquan::Bool=true)
+    T = eltype(vs)
     nv  = length(vs)
     p   = Vector{Float64}(undef, nv)
     c   = Vector{Float64}(undef, nv)
@@ -490,10 +491,11 @@ function test_evaluation(d::ContinuousUnivariateDistribution, vs::AbstractVector
             qtol = isa(d, InverseGaussian) ? 1.0e-4 : 1.0e-10
             qtol = isa(d, StudentizedRange) ? 1.0e-5 : qtol
             if p[i] > 1.0e-6
-                @test isapprox(quantile(d, c[i])    , v, atol=qtol * (abs(v) + 1.0))
-                @test isapprox(cquantile(d, cc[i])  , v, atol=qtol * (abs(v) + 1.0))
-                @test isapprox(invlogcdf(d, lc[i])  , v, atol=qtol * (abs(v) + 1.0))
-                @test isapprox(invlogccdf(d, lcc[i]), v, atol=qtol * (abs(v) + 1.0))
+                for (f, z) in ((quantile, c[i]), (cquantile, cc[i]), (invlogcdf, lc[i]), (invlogccdf, lcc[i]))
+                    fz = @inferred(f(d, z))
+                    @test fz isa T
+                    @test fz ≈ v atol = qtol * (abs(v) + 1.0)
+                end
             end
         end
     end