Skip to content

Commit

Permalink
Fix (log)(c)cdf with Inf, -Inf and NaN (#1348)
Browse files Browse the repository at this point in the history
  • Loading branch information
devmotion authored Jul 2, 2021
1 parent 63ca313 commit d608045
Show file tree
Hide file tree
Showing 38 changed files with 624 additions and 456 deletions.
13 changes: 4 additions & 9 deletions src/mixtures/mixturemodel.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down Expand Up @@ -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)
Expand All @@ -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
Expand Down
107 changes: 36 additions & 71 deletions src/truncate.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
33 changes: 18 additions & 15 deletions src/univariate/continuous/betaprime.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
15 changes: 6 additions & 9 deletions src/univariate/continuous/chi.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
21 changes: 13 additions & 8 deletions src/univariate/continuous/exponential.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down
21 changes: 12 additions & 9 deletions src/univariate/continuous/frechet.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
42 changes: 10 additions & 32 deletions src/univariate/continuous/generalizedextremevalue.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
28 changes: 12 additions & 16 deletions src/univariate/continuous/generalizedpareto.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
Loading

0 comments on commit d608045

Please sign in to comment.