Fix (log)(c)cdf with Inf, -Inf and NaN

src/mixtures/mixturemodel.jl

@@ -277,15 +277,12 @@ function insupport(d::AbstractMixtureModel, x::AbstractVector)
     return any(insupport(component(d, i), x) for (i, pi) in enumerate(p) if !iszero(pi))
 end
 
-function _cdf(d::UnivariateMixture, x::Real)
+function cdf(d::UnivariateMixture, x::Real)
     p = probs(d)
     r = sum(pi * cdf(component(d, i), x) for (i, pi) in enumerate(p) if !iszero(pi))
     return r
 end
 
-cdf(d::UnivariateMixture{Continuous}, x::Real) = _cdf(d, x)
-cdf(d::UnivariateMixture{Discrete}, x::Integer) = _cdf(d, x)
-
 function _mixpdf1(d::AbstractMixtureModel, x)
     p = probs(d)
     return sum(pi * pdf(component(d, i), x) for (i, pi) in enumerate(p) if !iszero(pi))
@@ -362,10 +359,8 @@ function _mixlogpdf!(r::AbstractArray, d::AbstractMixtureModel, x)
     return r
 end
 
-pdf(d::UnivariateMixture{Continuous}, x::Real) = _mixpdf1(d, x)
-pdf(d::UnivariateMixture{Discrete}, x::Int) = _mixpdf1(d, x)
-logpdf(d::UnivariateMixture{Continuous}, x::Real) = _mixlogpdf1(d, x)
-logpdf(d::UnivariateMixture{Discrete}, x::Int) = _mixlogpdf1(d, x)
+pdf(d::UnivariateMixture, x::Real) = _mixpdf1(d, x)
+logpdf(d::UnivariateMixture, x::Real) = _mixlogpdf1(d, x)
 
 _pdf!(r::AbstractArray, d::UnivariateMixture{Discrete}, x::UnitRange) = _mixpdf!(r, d, x)
 _pdf!(r::AbstractArray, d::UnivariateMixture, x::AbstractArray) = _mixpdf!(r, d, x)
@@ -374,7 +369,7 @@ _logpdf!(r::AbstractArray, d::UnivariateMixture, x::AbstractArray) = _mixlogpdf!
 _pdf(d::MultivariateMixture, x::AbstractVector) = _mixpdf1(d, x)
 _logpdf(d::MultivariateMixture, x::AbstractVector) = _mixlogpdf1(d, x)
 _pdf!(r::AbstractArray, d::MultivariateMixture, x::AbstractMatrix) = _mixpdf!(r, d, x)
-_lodpdf!(r::AbstractArray, d::MultivariateMixture, x::AbstractMatrix) = _mixlogpdf!(r, d, x)
+_logpdf!(r::AbstractArray, d::MultivariateMixture, x::AbstractMatrix) = _mixlogpdf!(r, d, x)
 
 
 ## component-wise pdf and logpdf

src/truncate.jl

@@ -91,94 +91,59 @@ end
 
 quantile(d::Truncated, p::Real) = quantile(d.untruncated, d.lcdf + p * d.tp)
 
-function _pdf(d::Truncated, x::T) where {T<:Real}
-    if d.lower <= x <= d.upper
-        pdf(d.untruncated, x) / d.tp
-    else
-        zero(T)
-    end
-end
-
-function pdf(d::Truncated{<:ContinuousUnivariateDistribution}, x::T) where {T<:Real}
-    _pdf(d, float(x))
-end
-
-function pdf(d::Truncated{D}, x::T) where {D<:DiscreteUnivariateDistribution, T<:Real}
-    isinteger(x) || return zero(float(T))
-    _pdf(d, x)
+function pdf(d::Truncated, x::Real)
+    result = pdf(d.untruncated, x) / d.tp
+    return d.lower <= x <= d.upper ? result : zero(result)
 end
 
-function pdf(d::Truncated{D}, x::T) where {D<:DiscreteUnivariateDistribution, T<:Integer}
-    _pdf(d, float(x))
+function logpdf(d::Truncated, x::Real)
+    result = logpdf(d.untruncated, x) - d.logtp
+    return d.lower <= x <= d.upper ? result : oftype(result, -Inf)
 end
 
-function _logpdf(d::Truncated, x::T) where {T<:Real}
-    if d.lower <= x <= d.upper
-        logpdf(d.untruncated, x) - d.logtp
+function cdf(d::Truncated, x::Real)
+    result = (cdf(d.untruncated, x) - d.lcdf) / d.tp
+    return if x < d.lower
+        zero(result)
+    elseif x >= d.upper
+        one(result)
     else
-        TF = float(T)
-        -TF(Inf)
+        result
     end
 end
 
-function logpdf(d::Truncated{D}, x::T) where {D<:DiscreteUnivariateDistribution, T<:Real}
-    TF = float(T)
-    isinteger(x) || return -TF(Inf)
-    return _logpdf(d, x)
-end
-
-function logpdf(d::Truncated{D}, x::Integer) where {D<:DiscreteUnivariateDistribution}
-    _logpdf(d, x)
-end
-
-function logpdf(d::Truncated{D, Continuous}, x::T) where {D<:ContinuousUnivariateDistribution, T<:Real}
-    _logpdf(d, x)
-end
-
-# fallback to avoid method ambiguities
-_cdf(d::Truncated, x::T) where {T<:Real} =
-    x <= d.lower ? zero(T) :
-    x >= d.upper ? one(T) :
-    (cdf(d.untruncated, x) - d.lcdf) / d.tp
-
-cdf(d::Truncated, x::Real) = _cdf(d, x)
-cdf(d::Truncated, x::Integer) = _cdf(d, float(x)) # float conversion for stability
-
-function _logcdf(d::Truncated, x::T) where {T<:Real}
-    TF = float(T)
-    if x <= d.lower
-        -TF(Inf)
+function logcdf(d::Truncated, x::Real)
+    result = logsubexp(logcdf(d.untruncated, x), log(d.lcdf)) - d.logtp
+    return if x < d.lower
+        oftype(result, -Inf)
     elseif x >= d.upper
-        zero(TF)
+        zero(result)
     else
-        log(cdf(d.untruncated, x) - d.lcdf) - d.logtp
+        result
     end
 end
 
-logcdf(d::Truncated, x::Real) = _logcdf(d, x)
-logcdf(d::Truncated, x::Integer) = _logcdf(d, x)
-
-_ccdf(d::Truncated, x::T) where {T<:Real} =
-    x <= d.lower ? one(T) :
-    x >= d.upper ? zero(T) :
-    (d.ucdf - cdf(d.untruncated, x)) / d.tp
-
-ccdf(d::Truncated, x::Real) = _ccdf(d, x)
-ccdf(d::Truncated, x::Integer) = _ccdf(d, float(x))
-
-function _logccdf(d::Truncated, x::T) where {T<:Real}
-    TF = float(T)
-    if x <= d.lower
-        zero(TF)
-    elseif x >= d.upper
-        -TF(Inf)
+function ccdf(d::Truncated, x::Real)
+    result = (d.ucdf - cdf(d.untruncated, x)) / d.tp
+    return if x <= d.lower
+        one(result)
+    elseif x > d.upper
+        zero(result)
     else
-        log(d.ucdf - cdf(d.untruncated, x)) - d.logtp
+        result
     end
 end
 
-logccdf(d::Truncated, x::Real) = _logccdf(d, x)
-logccdf(d::Truncated, x::Integer) = _logccdf(d, x)
+function logccdf(d::Truncated, x::Real)
+    result = logsubexp(logccdf(d.untruncated, x), log1p(-d.ucdf)) - d.logtp
+    return if x <= d.lower
+        zero(result)
+    elseif x > d.upper
+        oftype(result, -Inf)
+    else
+        result
+    end
+end
 
 ## random number generation
 

src/univariate/continuous/betaprime.jl

@@ -85,25 +85,28 @@ end
 
 #### Evaluation
 
-function logpdf(d::BetaPrime{T}, x::Real) where T<:Real
-    (α, β) = params(d)
-    if x < 0
-        T(-Inf)
-    else
-        (α - 1) * log(x) - (α + β) * log1p(x) - logbeta(α, β)
-    end
+function logpdf(d::BetaPrime, x::Real)
+    α, β = params(d)
+    _x = max(0, x)
+    z = xlogy(α - 1, _x) - (α + β) * log1p(_x) - logbeta(α, β)
+    return x < 0 ? oftype(z, -Inf) : z
 end
 
-cdf(d::BetaPrime{T}, x::Real) where {T<:Real} = x <= 0 ? zero(T) : betacdf(d.α, d.β, x / (1 + x))
-ccdf(d::BetaPrime{T}, x::Real) where {T<:Real} = x <= 0 ? one(T) : betaccdf(d.α, d.β, x / (1 + x))
-logcdf(d::BetaPrime{T}, x::Real) where {T<:Real} =  x <= 0 ? T(-Inf) : betalogcdf(d.α, d.β, x / (1 + x))
-logccdf(d::BetaPrime{T}, x::Real) where {T<:Real} =  x <= 0 ? zero(T) : betalogccdf(d.α, d.β, x / (1 + x))
+function zval(::BetaPrime, x::Real)
+    y = max(x, 0)
+    z = y / (1 + y)
+    # map `Inf` to `Inf` (otherwise it returns `NaN`)
+    return isinf(x) && x > 0 ? oftype(z, Inf) : z
+end
+xval(::BetaPrime, z::Real) = z / (1 - z)
 
-quantile(d::BetaPrime, p::Real) = (x = betainvcdf(d.α, d.β, p); x / (1 - x))
-cquantile(d::BetaPrime, p::Real) = (x = betainvccdf(d.α, d.β, p); x / (1 - x))
-invlogcdf(d::BetaPrime, p::Real) = (x = betainvlogcdf(d.α, d.β, p); x / (1 - x))
-invlogccdf(d::BetaPrime, p::Real) = (x = betainvlogccdf(d.α, d.β, p); x / (1 - x))
+for f in (:cdf, :ccdf, :logcdf, :logccdf)
+    @eval $f(d::BetaPrime, x::Real) = $f(Beta(d.α, d.β; check_args=false), zval(d, x))
+end
 
+for f in (:quantile, :cquantile, :invlogcdf, :invlogccdf)
+    @eval $f(d::BetaPrime, p::Real) = xval(d, $f(Beta(d.α, d.β; check_args=false), p))
+end
 
 #### Sampling
 

src/univariate/continuous/chi.jl

@@ -83,16 +83,13 @@ logpdf(d::Chi, x::Real) = (ν = d.ν;
 
 gradlogpdf(d::Chi{T}, x::Real) where {T<:Real} = x >= 0 ? (d.ν - 1) / x - x : zero(T)
 
-cdf(d::Chi, x::Real) = chisqcdf(d.ν, x^2)
-ccdf(d::Chi, x::Real) = chisqccdf(d.ν, x^2)
-logcdf(d::Chi, x::Real) = chisqlogcdf(d.ν, x^2)
-logccdf(d::Chi, x::Real) = chisqlogccdf(d.ν, x^2)
-
-quantile(d::Chi, p::Real) = sqrt(chisqinvcdf(d.ν, p))
-cquantile(d::Chi, p::Real) = sqrt(chisqinvccdf(d.ν, p))
-invlogcdf(d::Chi, p::Real) = sqrt(chisqinvlogcdf(d.ν, p))
-invlogccdf(d::Chi, p::Real) = sqrt(chisqinvlogccdf(d.ν, p))
+for f in (:cdf, :ccdf, :logcdf, :logccdf)
+    @eval $f(d::Chi, x::Real) = $f(Chisq(d.ν; check_args=false), max(x, 0)^2)
+end
 
+for f in (:quantile, :cquantile, :invlogcdf, :invlogccdf)
+    @eval $f(d::Chi, p::Real) = sqrt($f(Chisq(d.ν; check_args=false), p))
+end
 
 #### Sampling
 

src/univariate/continuous/exponential.jl

@@ -62,19 +62,24 @@ entropy(d::Exponential{T}) where {T} = one(T) + log(d.θ)
 
 #### Evaluation
 
-zval(d::Exponential, x::Real) = x / d.θ
+zval(d::Exponential, x::Real) = max(x / d.θ, 0)
 xval(d::Exponential, z::Real) = z * d.θ
 
-pdf(d::Exponential, x::Real) = (λ = rate(d); x < 0 ? zero(λ) : λ * exp(-λ * x))
-function logpdf(d::Exponential{T}, x::Real) where T<:Real
+function pdf(d::Exponential, x::Real)
     λ = rate(d)
-    x < 0 ? -T(Inf) : log(λ) - λ * x
+    z = λ * exp(-λ * max(x, 0))
+    return x < 0 ? zero(z) : z
+end
+function logpdf(d::Exponential, x::Real)
+    λ = rate(d)
+    z = log(λ) - λ * x
+    return x < 0 ? oftype(z, -Inf) : z
 end
 
-cdf(d::Exponential{T}, x::Real) where {T<:Real} = x > 0 ? -expm1(-zval(d, x)) : zero(T)
-ccdf(d::Exponential{T}, x::Real) where {T<:Real} = x > 0 ? exp(-zval(d, x)) : zero(T)
-logcdf(d::Exponential{T}, x::Real) where {T<:Real} = x > 0 ? log1mexp(-zval(d, x)) : -T(Inf)
-logccdf(d::Exponential{T}, x::Real) where {T<:Real} = x > 0 ? -zval(d, x) : zero(T)
+cdf(d::Exponential, x::Real) = -expm1(-zval(d, x))
+ccdf(d::Exponential, x::Real) = exp(-zval(d, x))
+logcdf(d::Exponential, x::Real) = log1mexp(-zval(d, x))
+logccdf(d::Exponential, x::Real) = -zval(d, x)
 
 quantile(d::Exponential, p::Real) = -xval(d, log1p(-p))
 cquantile(d::Exponential, p::Real) = -xval(d, log(p))

src/univariate/continuous/frechet.jl

@@ -119,15 +119,18 @@ function logpdf(d::Frechet{T}, x::Real) where T<:Real
     end
 end
 
-cdf(d::Frechet{T}, x::Real) where {T<:Real} = x > 0 ? exp(-((d.θ / x) ^ d.α)) : zero(T)
-ccdf(d::Frechet{T}, x::Real) where {T<:Real} = x > 0 ? -expm1(-((d.θ / x) ^ d.α)) : one(T)
-logcdf(d::Frechet{T}, x::Real) where {T<:Real} = x > 0 ? -(d.θ / x) ^ d.α : -T(Inf)
-logccdf(d::Frechet{T}, x::Real) where {T<:Real} = x > 0 ? log1mexp(-((d.θ / x) ^ d.α)) : zero(T)
-
-quantile(d::Frechet, p::Real) = d.θ * (-log(p)) ^ (-1 / d.α)
-cquantile(d::Frechet, p::Real) = d.θ * (-log1p(-p)) ^ (-1 / d.α)
-invlogcdf(d::Frechet, lp::Real) = d.θ * (-lp)^(-1 / d.α)
-invlogccdf(d::Frechet, lp::Real) = d.θ * (-log1mexp(lp))^(-1 / d.α)
+zval(d::Frechet, x::Real) = (d.θ / max(x, 0))^d.α
+xval(d::Frechet, z::Real) = d.θ * z^(- 1 / d.α)
+
+cdf(d::Frechet, x::Real) = exp(- zval(d, x))
+ccdf(d::Frechet, x::Real) = -expm1(- zval(d, x))
+logcdf(d::Frechet, x::Real) = - zval(d, x)
+logccdf(d::Frechet, x::Real) = log1mexp(- zval(d, x))
+
+quantile(d::Frechet, p::Real) = xval(d, -log(p))
+cquantile(d::Frechet, p::Real) = xval(d, -log1p(-p))
+invlogcdf(d::Frechet, lp::Real) = xval(d, -lp)
+invlogccdf(d::Frechet, lp::Real) = xval(d, -log1mexp(lp))
 
 function gradlogpdf(d::Frechet{T}, x::Real) where T<:Real
     (α, θ) = params(d)

src/univariate/continuous/generalizedextremevalue.jl

@@ -219,43 +219,21 @@ function pdf(d::GeneralizedExtremeValue{T}, x::Real) where T<:Real
     end
 end
 
-function logcdf(d::GeneralizedExtremeValue{T}, x::Real) where T<:Real
-    if insupport(d, x)
-        (μ, σ, ξ) = params(d)
-
-        z = (x - μ) / σ # Normalise x.
-        if abs(ξ) < eps(one(ξ)) # ξ == 0
-            return -exp(- z)
-        else
-            return - (1 + z * ξ) ^ ( -1/ξ)
-        end
-    elseif x <= minimum(d)
-        return -T(Inf)
-    else
-        return zero(T)
-    end
-end
-
-function cdf(d::GeneralizedExtremeValue{T}, x::Real) where T<:Real
-    if insupport(d, x)
-        (μ, σ, ξ) = params(d)
-
-        z = (x - μ) / σ # Normalise x.
-        if abs(ξ) < eps(one(ξ)) # ξ == 0
-            t = exp(-z)
-        else
-            t = (1 + z * ξ) ^ (-1/ξ)
-        end
-        return exp(- t)
-    elseif x <= minimum(d)
-        return zero(T)
+function logcdf(d::GeneralizedExtremeValue, x::Real)
+    μ, σ, ξ = params(d)
+    z = (x - μ) / σ
+    return if abs(ξ) < eps(one(ξ)) # ξ == 0
+        -exp(- z)
     else
-        return one(T)
+        # y(x) = y(bound) = 0 if x is not in the support with lower/upper bound
+        y = max(1 + z * ξ, 0)
+        - y^(-1/ξ)
     end
 end
+cdf(d::GeneralizedExtremeValue, x::Real) = exp(logcdf(d, x))
 
-logccdf(d::GeneralizedExtremeValue, x::Real) = log1p(- cdf(d, x))
 ccdf(d::GeneralizedExtremeValue, x::Real) = - expm1(logcdf(d, x))
+logccdf(d::GeneralizedExtremeValue, x::Real) = log1mexp(logcdf(d, x))
 
 
 #### Sampling

src/univariate/continuous/generalizedpareto.jl

@@ -132,26 +132,22 @@ function logpdf(d::GeneralizedPareto{T}, x::Real) where T<:Real
     return p
 end
 
-function logccdf(d::GeneralizedPareto{T}, x::Real) where T<:Real
-    (μ, σ, ξ) = params(d)
-
-    # The logccdf is log(0) outside the support range.
-    p = -T(Inf)
-
-    if x >= μ
-        z = (x - μ) / σ
-        if abs(ξ) < eps()
-            p = -z
-        elseif ξ > 0 || (ξ < 0 && x < maximum(d))
-            p = (-1 / ξ) * log1p(z * ξ)
-        end
+function logccdf(d::GeneralizedPareto, x::Real)
+    μ, σ, ξ = params(d)
+    z = max((x - μ) / σ, 0) # z(x) = z(μ) = 0 if x < μ (lower bound)
+    return if abs(ξ) < eps(one(ξ)) # ξ == 0
+        -z
+    elseif ξ < 0
+        # y(x) = y(μ - σ / ξ) = -1 if x > μ - σ / ξ (upper bound)
+        -log1p(max(z * ξ, -1)) / ξ
+    else
+        -log1p(z * ξ) / ξ
     end
-
-    return p
 end
-
 ccdf(d::GeneralizedPareto, x::Real) = exp(logccdf(d, x))
+
 cdf(d::GeneralizedPareto, x::Real) = -expm1(logccdf(d, x))
+logcdf(d::GeneralizedPareto, x::Real) = log1mexp(logccdf(d, x))
 
 function quantile(d::GeneralizedPareto{T}, p::Real) where T<:Real
     (μ, σ, ξ) = params(d)