From 7ce0216cfffaecb4ba34f21f04c578550cbd817d Mon Sep 17 00:00:00 2001 From: Kris De Meyer Date: Sat, 30 Jan 2016 19:52:34 +0000 Subject: [PATCH] This fixes errors in the boundary conditions for univariate distributions. These errors were partly caused because for some distributions the boundary points of the support are not in the support itself. I solved this by allowing the support to be a closed, half-open or open interval. The default is a closed interval (which corresponds to the old behaviour), meaning that only those distributions with a half-open or open support interval needed to be ammended. The core of this functionality is in the file src/univariates.jl. For half-open or open intervals, the macro @distr_boundaries specifies which boundaries of the interval are open or closed. See e.g., in src/univariate/continuous/lognormal.jl for an example. Tests have been added to test/utils.jl and src/testutils.jl A second set of errors came from sloppy implementations in pdf, logpdf, cdf, ccdf, logcdf or logccdf functions of specific distributions that did not explicitly test if the presented values fall inside or outside the support. Tests that were added in src/testutils.jl flagged - I believe - all instances of such errors, and the errors have been fixed. I spotted 2 additional errors in the support() function in src/univariates.jl. There were no tests for this function so I added them to src/testutils.jl/test_support() --- src/testutils.jl | 50 ++++++++++++++++++- src/univariate/continuous/betaprime.jl | 9 +++- src/univariate/continuous/biweight.jl | 4 +- src/univariate/continuous/chi.jl | 20 +++++--- src/univariate/continuous/cosine.jl | 20 ++++++-- src/univariate/continuous/epanechnikov.jl | 4 +- src/univariate/continuous/exponential.jl | 2 +- .../continuous/generalizedpareto.jl | 15 +++--- src/univariate/continuous/inversegamma.jl | 17 ++++--- src/univariate/continuous/inversegaussian.jl | 13 ++--- src/univariate/continuous/levy.jl | 25 +++++++--- src/univariate/continuous/lognormal.jl | 18 +++---- src/univariate/continuous/rayleigh.jl | 18 +++++-- src/univariate/continuous/triweight.jl | 2 +- src/univariate/discrete/betabinomial.jl | 14 ++++-- src/univariates.jl | 39 ++++++++++++--- test/utils.jl | 19 +++++-- 17 files changed, 213 insertions(+), 76 deletions(-) diff --git a/src/testutils.jl b/src/testutils.jl index 0caf7451a..76ee03def 100644 --- a/src/testutils.jl +++ b/src/testutils.jl @@ -268,15 +268,61 @@ function test_support(d::UnivariateDistribution, vs::AbstractVector) end @test all(insupport(d, vs)) + @test lowerboundary(d) == :closed || lowerboundary(d) == :open + @test upperboundary(d) == :closed || upperboundary(d) == :open + + sp = support(d) + if isa(d,ContinuousUnivariateDistribution) + @test minimum(sp) == minimum(d) + @test maximum(sp) == maximum(d) + @test lowerboundary(sp) == lowerboundary(d) + @test upperboundary(sp) == upperboundary(d) + @test lowercomparator(sp) == lowercomparator(d) + @test uppercomparator(sp) == uppercomparator(d) + end + + if isa(d,DiscreteUnivariateDistribution) + @test (isfinite(minimum(d)) && minimum(d) == minimum(sp)) || (minimum(sp) == typemin(Int)) + @test (isfinite(minimum(d)) && maximum(d) == maximum(sp)) || (maximum(sp) == typemax(Int)) + end + if islowerbounded(d) @test isfinite(minimum(d)) - @test insupport(d, minimum(d)) + @test (lowerboundary(d) == :closed && insupport(d, minimum(d))) || (lowerboundary(d) == :open && !insupport(d, minimum(d))) @test !insupport(d, minimum(d)-1) + if lowerboundary(d) == :open + @test pdf(d,minimum(d)) == 0.0 + @test logpdf(d,minimum(d)) == -Inf + @test cdf(d,minimum(d)) == 0.0 + @test logcdf(d,minimum(d)) == -Inf + @test ccdf(d,minimum(d)) == 1.0 + @test logccdf(d,minimum(d)) == 0.0 + end + @test pdf(d,minimum(d)-1) == 0.0 + @test logpdf(d,minimum(d)-1) == -Inf + @test cdf(d,minimum(d)-1) == 0.0 + @test logcdf(d,minimum(d)-1) == -Inf + @test ccdf(d,minimum(d)-1) == 1.0 + @test logccdf(d,minimum(d)-1) == 0.0 end if isupperbounded(d) @test isfinite(maximum(d)) - @test insupport(d, maximum(d)) + @test (upperboundary(d) == :closed && insupport(d, maximum(d))) || (upperboundary(d) == :open && !insupport(d, maximum(d))) @test !insupport(d, maximum(d)+1) + if upperboundary(d) == :open + @test pdf(d,maximum(d)) == 0.0 + @test logpdf(d,maximum(d)) == -Inf + end + @test cdf(d,maximum(d)) == 1.0 + @test ccdf(d,maximum(d)) == 0.0 + @test logcdf(d,maximum(d)) == 0.0 + @test logccdf(d,maximum(d)) == -Inf + @test pdf(d,maximum(d)+1) == 0.0 + @test logpdf(d,maximum(d)+1) == -Inf + @test cdf(d,maximum(d)+1) == 1.0 + @test logcdf(d,maximum(d)+1) == 0.0 + @test ccdf(d,maximum(d)+1) == 0.0 + @test logccdf(d,maximum(d)+1) == -Inf end @test isbounded(d) == (isupperbounded(d) && islowerbounded(d)) diff --git a/src/univariate/continuous/betaprime.jl b/src/univariate/continuous/betaprime.jl index 47e0f19ae..756ca5d09 100644 --- a/src/univariate/continuous/betaprime.jl +++ b/src/univariate/continuous/betaprime.jl @@ -11,6 +11,7 @@ immutable BetaPrime <: ContinuousUnivariateDistribution end @distr_support BetaPrime 0.0 Inf +@distr_boundaries BetaPrime :open :closed #### Parameters @@ -42,8 +43,12 @@ end #### Evaluation function logpdf(d::BetaPrime, x::Float64) - (α, β) = params(d) - (α - 1.0) * log(x) - (α + β) * log1p(x) - lbeta(α, β) + if insupport(d, x) + (α, β) = params(d) + (α - 1.0) * log(x) - (α + β) * log1p(x) - lbeta(α, β) + else + return -Inf + end end pdf(d::BetaPrime, x::Float64) = exp(logpdf(d, x)) diff --git a/src/univariate/continuous/biweight.jl b/src/univariate/continuous/biweight.jl index a770e9d9d..757704d4d 100644 --- a/src/univariate/continuous/biweight.jl +++ b/src/univariate/continuous/biweight.jl @@ -36,8 +36,8 @@ end function ccdf(d::Biweight, x::Float64) u = (d.μ - x) / d.σ - u <= -1.0 ? 1.0 : - u >= 1.0 ? 0.0 : + u <= -1.0 ? 0.0 : + u >= 1.0 ? 1.0 : 0.0625 * (u + 1.0)^3 * @horner(u,8.0,-9.0,3.0) end diff --git a/src/univariate/continuous/chi.jl b/src/univariate/continuous/chi.jl index cdd3bd4bd..c2df345d2 100644 --- a/src/univariate/continuous/chi.jl +++ b/src/univariate/continuous/chi.jl @@ -5,6 +5,7 @@ immutable Chi <: ContinuousUnivariateDistribution end @distr_support Chi 0.0 Inf +@distr_boundaries Chi :open :closed #### Parameters @@ -45,16 +46,21 @@ end pdf(d::Chi, x::Float64) = exp(logpdf(d, x)) -logpdf(d::Chi, x::Float64) = (ν = d.ν; - (1.0 - 0.5 * ν) * logtwo + (ν - 1.0) * log(x) - 0.5 * x^2 - lgamma(0.5 * ν) -) +function logpdf(d::Chi, x::Float64) + if insupport(d, x) + ν = d.ν + return (1.0 - 0.5 * ν) * logtwo + (ν - 1.0) * log(x) - 0.5 * x^2 - lgamma(0.5 * ν) + else + return -Inf + end +end gradlogpdf(d::Chi, x::Float64) = x >= 0.0 ? (d.ν - 1.0) / x - x : 0.0 -cdf(d::Chi, x::Float64) = chisqcdf(d.ν, x^2) -ccdf(d::Chi, x::Float64) = chisqccdf(d.ν, x^2) -logcdf(d::Chi, x::Float64) = chisqlogcdf(d.ν, x^2) -logccdf(d::Chi, x::Float64) = chisqlogccdf(d.ν, x^2) +cdf(d::Chi, x::Float64) = insupport(d,x)?chisqcdf(d.ν, x^2):0.0 +ccdf(d::Chi, x::Float64) = insupport(d,x)?chisqccdf(d.ν, x^2):1.0 +logcdf(d::Chi, x::Float64) = insupport(d,x)?chisqlogcdf(d.ν, x^2):-Inf +logccdf(d::Chi, x::Float64) = insupport(d,x)?chisqlogccdf(d.ν, x^2):0.0 quantile(d::Chi, p::Float64) = sqrt(chisqinvcdf(d.ν, p)) cquantile(d::Chi, p::Float64) = sqrt(chisqinvccdf(d.ν, p)) diff --git a/src/univariate/continuous/cosine.jl b/src/univariate/continuous/cosine.jl index 9f3ee9364..56b33c8ec 100644 --- a/src/univariate/continuous/cosine.jl +++ b/src/univariate/continuous/cosine.jl @@ -52,13 +52,25 @@ end logpdf(d::Cosine, x::Float64) = insupport(d, x) ? log(pdf(d, x)) : -Inf function cdf(d::Cosine, x::Float64) - z = (x - d.μ) / d.σ - 0.5 * (1.0 + z + sinpi(z) * invπ) + if insupport(d, x) + z = (x - d.μ) / d.σ + return 0.5 * (1.0 + z + sinpi(z) * invπ) + elseif x < minimum(d) + return 0.0 + else + return 1.0 + end end function ccdf(d::Cosine, x::Float64) - nz = (d.μ - x) / d.σ - 0.5 * (1.0 + nz + sinpi(nz) * invπ) + if insupport(d, x) + nz = (d.μ - x) / d.σ + return 0.5 * (1.0 + nz + sinpi(nz) * invπ) + elseif x < minimum(d) + return 1.0 + else + return 0.0 + end end quantile(d::Cosine, p::Float64) = quantile_bisect(d, p) diff --git a/src/univariate/continuous/epanechnikov.jl b/src/univariate/continuous/epanechnikov.jl index 4e9b066bd..8da69e71b 100644 --- a/src/univariate/continuous/epanechnikov.jl +++ b/src/univariate/continuous/epanechnikov.jl @@ -39,8 +39,8 @@ end function ccdf(d::Epanechnikov, x::Float64) u = (d.μ - x) / d.σ - u <= -1 ? 1.0 : - u >= 1 ? 0.0 : + u <= -1 ? 0.0 : + u >= 1 ? 1.0 : 0.5 + u * (0.75 - 0.25 * u^2) end diff --git a/src/univariate/continuous/exponential.jl b/src/univariate/continuous/exponential.jl index 76d1ce788..9e2806321 100644 --- a/src/univariate/continuous/exponential.jl +++ b/src/univariate/continuous/exponential.jl @@ -38,7 +38,7 @@ pdf(d::Exponential, x::Float64) = (λ = rate(d); x < 0.0 ? 0.0 : λ * exp(-λ * logpdf(d::Exponential, x::Float64) = (λ = rate(d); x < 0.0 ? -Inf : log(λ) - λ * x) cdf(d::Exponential, x::Float64) = x > 0.0 ? -expm1(-zval(d, x)) : 0.0 -ccdf(d::Exponential, x::Float64) = x > 0.0 ? exp(-zval(d, x)) : 0.0 +ccdf(d::Exponential, x::Float64) = x > 0.0 ? exp(-zval(d, x)) : 1.0 logcdf(d::Exponential, x::Float64) = x > 0.0 ? log1mexp(-zval(d, x)) : -Inf logccdf(d::Exponential, x::Float64) = x > 0.0 ? -zval(d, x) : 0.0 diff --git a/src/univariate/continuous/generalizedpareto.jl b/src/univariate/continuous/generalizedpareto.jl index a1c861201..259d57c10 100644 --- a/src/univariate/continuous/generalizedpareto.jl +++ b/src/univariate/continuous/generalizedpareto.jl @@ -11,8 +11,7 @@ immutable GeneralizedPareto <: ContinuousUnivariateDistribution GeneralizedPareto() = new(1.0, 1.0, 0.0) end -minimum(d::GeneralizedPareto) = d.μ -maximum(d::GeneralizedPareto) = d.ξ < 0.0 ? d.μ - d.σ / d.ξ : Inf +@distr_support GeneralizedPareto d.μ d.ξ<0.0?d.μ-d.σ/d.ξ:Inf #### Parameters @@ -74,9 +73,9 @@ function logpdf(d::GeneralizedPareto, x::Float64) # The logpdf is log(0) outside the support range. p = -Inf - if x >= μ + if insupport(d,x) z = (x - μ) / σ - if abs(ξ) < eps() + if abs(ξ) < eps(x) p = -z - log(σ) elseif ξ > 0.0 || (ξ < 0.0 && x < maximum(d)) p = (-1.0 - 1.0 / ξ) * log1p(z * ξ) - log(σ) @@ -91,10 +90,9 @@ pdf(d::GeneralizedPareto, x::Float64) = exp(logpdf(d, x)) function logccdf(d::GeneralizedPareto, x::Float64) (ξ, σ, μ) = params(d) - # The logccdf is log(0) outside the support range. - p = -Inf + p = 0.0 - if x >= μ + if insupport(d,x) z = (x - μ) / σ if abs(ξ) < eps() p = -z @@ -102,6 +100,9 @@ function logccdf(d::GeneralizedPareto, x::Float64) p = (-1.0 / ξ) * log1p(z * ξ) end end + if x >= maximum(d) + p = -Inf + end return p end diff --git a/src/univariate/continuous/inversegamma.jl b/src/univariate/continuous/inversegamma.jl index 2d752876f..d052e3cb4 100644 --- a/src/univariate/continuous/inversegamma.jl +++ b/src/univariate/continuous/inversegamma.jl @@ -12,6 +12,7 @@ immutable InverseGamma <: ContinuousUnivariateDistribution end @distr_support InverseGamma 0.0 Inf +@distr_boundaries InverseGamma :open :closed #### Parameters @@ -55,14 +56,18 @@ end pdf(d::InverseGamma, x::Float64) = exp(logpdf(d, x)) function logpdf(d::InverseGamma, x::Float64) - (α, θ) = params(d) - α * log(θ) - lgamma(α) - (α + 1.0) * log(x) - θ / x + if insupport(d, x) + (α, θ) = params(d) + return α * log(θ) - lgamma(α) - (α + 1.0) * log(x) - θ / x + else + return -Inf + end end -cdf(d::InverseGamma, x::Float64) = ccdf(d.invd, 1.0 / x) -ccdf(d::InverseGamma, x::Float64) = cdf(d.invd, 1.0 / x) -logcdf(d::InverseGamma, x::Float64) = logccdf(d.invd, 1.0 / x) -logccdf(d::InverseGamma, x::Float64) = logcdf(d.invd, 1.0 / x) +cdf(d::InverseGamma, x::Float64) = insupport(d,x)?ccdf(d.invd, 1.0 / x):0.0 +ccdf(d::InverseGamma, x::Float64) = insupport(d,x)?cdf(d.invd, 1.0 / x):1.0 +logcdf(d::InverseGamma, x::Float64) = insupport(d,x)?logccdf(d.invd, 1.0 / x):-Inf +logccdf(d::InverseGamma, x::Float64) = insupport(d,x)?logcdf(d.invd, 1.0 / x):0.0 quantile(d::InverseGamma, p::Float64) = 1.0 / cquantile(d.invd, p) cquantile(d::InverseGamma, p::Float64) = 1.0 / quantile(d.invd, p) diff --git a/src/univariate/continuous/inversegaussian.jl b/src/univariate/continuous/inversegaussian.jl index 0eed7918b..0f4fdb399 100644 --- a/src/univariate/continuous/inversegaussian.jl +++ b/src/univariate/continuous/inversegaussian.jl @@ -11,6 +11,7 @@ immutable InverseGaussian <: ContinuousUnivariateDistribution end @distr_support InverseGaussian 0.0 Inf +@distr_boundaries InverseGaussian :open :closed #### Parameters @@ -39,7 +40,7 @@ end #### Evaluation function pdf(d::InverseGaussian, x::Float64) - if x > 0.0 + if insupport(d, x) μ, λ = params(d) return sqrt(λ / (twoπ * x^3)) * exp(-λ * (x - μ)^2 / (2.0 * μ^2 * x)) else @@ -48,7 +49,7 @@ function pdf(d::InverseGaussian, x::Float64) end function logpdf(d::InverseGaussian, x::Float64) - if x > 0.0 + if insupport(d, x) μ, λ = params(d) return 0.5 * (log(λ) - (log2π + 3.0 * log(x)) - λ * (x - μ)^2 / (μ^2 * x)) else @@ -57,7 +58,7 @@ function logpdf(d::InverseGaussian, x::Float64) end function cdf(d::InverseGaussian, x::Float64) - if x > 0.0 + if insupport(d, x) μ, λ = params(d) u = sqrt(λ / x) v = x / μ @@ -68,7 +69,7 @@ function cdf(d::InverseGaussian, x::Float64) end function ccdf(d::InverseGaussian, x::Float64) - if x > 0.0 + if insupport(d, x) μ, λ = params(d) u = sqrt(λ / x) v = x / μ @@ -79,7 +80,7 @@ function ccdf(d::InverseGaussian, x::Float64) end function logcdf(d::InverseGaussian, x::Float64) - if x > 0.0 + if insupport(d, x) μ, λ = params(d) u = sqrt(λ / x) v = x / μ @@ -92,7 +93,7 @@ function logcdf(d::InverseGaussian, x::Float64) end function logccdf(d::InverseGaussian, x::Float64) - if x > 0.0 + if insupport(d, x) μ, λ = params(d) u = sqrt(λ / x) v = x / μ diff --git a/src/univariate/continuous/levy.jl b/src/univariate/continuous/levy.jl index dab2519ac..04eb23463 100644 --- a/src/univariate/continuous/levy.jl +++ b/src/univariate/continuous/levy.jl @@ -8,6 +8,7 @@ immutable Levy <: ContinuousUnivariateDistribution end @distr_support Levy d.μ Inf +@distr_boundaries Levy :open :open #### Parameters @@ -33,19 +34,27 @@ median(d::Levy) = d.μ + d.σ / 0.4549364231195728 # 0.454... = (2.0 * erfcinv( #### Evaluation function pdf(d::Levy, x::Float64) - μ, σ = params(d) - z = x - μ - (sqrt(σ) / sqrt2π) * exp((-σ) / (2.0 * z)) / z^1.5 + if (insupport(d,x)) + μ, σ = params(d) + z = x - μ + return (sqrt(σ) / sqrt2π) * exp((-σ) / (2.0 * z)) / z^1.5 + else + return zero(x) + end end function logpdf(d::Levy, x::Float64) - μ, σ = params(d) - z = x - μ - 0.5 * (log(σ) - log2π - σ / z - 3.0 * log(z)) + if (insupport(d,x)) + μ, σ = params(d) + z = x - μ + return 0.5 * (log(σ) - log2π - σ / z - 3.0 * log(z)) + else + return -Inf + end end -cdf(d::Levy, x::Float64) = erfc(sqrt(d.σ / (2.0 * (x - d.μ)))) -ccdf(d::Levy, x::Float64) = erf(sqrt(d.σ / (2.0 * (x - d.μ)))) +cdf(d::Levy, x::Float64) = insupport(d,x)?erfc(sqrt(d.σ / (2.0 * (x - d.μ)))):0.0 +ccdf(d::Levy, x::Float64) = insupport(d,x)?erf(sqrt(d.σ / (2.0 * (x - d.μ)))):1.0 quantile(d::Levy, p::Float64) = d.μ + d.σ / (2.0 * erfcinv(p)^2) cquantile(d::Levy, p::Float64) = d.μ + d.σ / (2.0 * erfinv(p)^2) diff --git a/src/univariate/continuous/lognormal.jl b/src/univariate/continuous/lognormal.jl index a5f210380..8bfcf418c 100644 --- a/src/univariate/continuous/lognormal.jl +++ b/src/univariate/continuous/lognormal.jl @@ -8,7 +8,7 @@ immutable LogNormal <: ContinuousUnivariateDistribution end @distr_support LogNormal 0.0 Inf - +@distr_boundaries LogNormal :open :closed #### Parameters @@ -53,20 +53,20 @@ end #### Evalution -pdf(d::LogNormal, x::Float64) = normpdf(d.μ, d.σ, log(x)) / x +pdf(d::LogNormal, x::Float64) = insupport(d,x) ? normpdf(d.μ, d.σ, log(x))/x : zero(x) function logpdf(d::LogNormal, x::Float64) - if !insupport(d, x) - return -Inf - else + if insupport(d, x) lx = log(x) return normlogpdf(d.μ, d.σ, lx) - lx + else + return -Inf end end -cdf(d::LogNormal, x::Float64) = x > 0.0 ? normcdf(d.μ, d.σ, log(x)) : 0.0 -ccdf(d::LogNormal, x::Float64) = x > 0.0 ? normccdf(d.μ, d.σ, log(x)) : 1.0 -logcdf(d::LogNormal, x::Float64) = x > 0.0 ? normlogcdf(d.μ, d.σ, log(x)) : -Inf -logccdf(d::LogNormal, x::Float64) = x > 0.0 ? normlogccdf(d.μ, d.σ, log(x)) : 0.0 +cdf(d::LogNormal, x::Float64) = insupport(d, x) ? normcdf(d.μ, d.σ, log(x)) : 0.0 +ccdf(d::LogNormal, x::Float64) = insupport(d, x) ? normccdf(d.μ, d.σ, log(x)) : 1.0 +logcdf(d::LogNormal, x::Float64) = insupport(d, x) ? normlogcdf(d.μ, d.σ, log(x)) : -Inf +logccdf(d::LogNormal, x::Float64) = insupport(d, x) ? normlogccdf(d.μ, d.σ, log(x)) : 0.0 quantile(d::LogNormal, q::Float64) = exp(norminvcdf(d.μ, d.σ, q)) cquantile(d::LogNormal, q::Float64) = exp(norminvccdf(d.μ, d.σ, q)) diff --git a/src/univariate/continuous/rayleigh.jl b/src/univariate/continuous/rayleigh.jl index 9f7481e43..22fd49d55 100644 --- a/src/univariate/continuous/rayleigh.jl +++ b/src/univariate/continuous/rayleigh.jl @@ -32,16 +32,24 @@ entropy(d::Rayleigh) = 0.942034242170793776 + log(d.σ) #### Evaluation function pdf(d::Rayleigh, x::Float64) - σ2 = d.σ^2 - x > 0.0 ? (x / σ2) * exp(- (x^2) / (2.0 * σ2)) : 0.0 + if insupport(d,x) + σ2 = d.σ^2 + return (x / σ2) * exp(- (x^2) / (2.0 * σ2)) + else + return 0.0 + end end function logpdf(d::Rayleigh, x::Float64) - σ2 = d.σ^2 - x > 0.0 ? log(x / σ2) - (x^2) / (2.0 * σ2) : -Inf + if (insupport(d,x)) + σ2 = d.σ^2 + return log(x / σ2) - (x^2) / (2.0 * σ2) + else + return -Inf + end end -logccdf(d::Rayleigh, x::Float64) = - (x^2) / (2.0 * d.σ^2) +logccdf(d::Rayleigh, x::Float64) = insupport(d,x)?-(x^2)/(2.0*d.σ^2):0.0 ccdf(d::Rayleigh, x::Float64) = exp(logccdf(d, x)) cdf(d::Rayleigh, x::Float64) = 1.0 - ccdf(d, x) diff --git a/src/univariate/continuous/triweight.jl b/src/univariate/continuous/triweight.jl index 82edbf7fd..40ad1b637 100644 --- a/src/univariate/continuous/triweight.jl +++ b/src/univariate/continuous/triweight.jl @@ -39,7 +39,7 @@ end function ccdf(d::Triweight, x::Real) u = (d.μ - x)/d.σ - u <= -1 ? 1.0 : u >= 1 ? 0.0 : 0.03125*(1+u)^4*@horner(u,16.0,-29.0,20.0,-5.0) + u <= -1 ? 0.0 : u >= 1 ? 1.0 : 0.03125*(1+u)^4*@horner(u,16.0,-29.0,20.0,-5.0) end @quantile_newton Triweight diff --git a/src/univariate/discrete/betabinomial.jl b/src/univariate/discrete/betabinomial.jl index 8ab8e194b..0249d76ed 100644 --- a/src/univariate/discrete/betabinomial.jl +++ b/src/univariate/discrete/betabinomial.jl @@ -49,11 +49,15 @@ function kurtosis(d::BetaBinomial) end function pdf(d::BetaBinomial, k::Int) - n, α, β = d.n, d.α, d.β - @compat choose = Float64(binomial(n, k)) - numerator = beta(k + α, n - k + β) - denominator = beta(α, β) - return choose * (numerator / denominator) + if insupport(d,k) + n, α, β = d.n, d.α, d.β + @compat choose = Float64(binomial(n, k)) + numerator = beta(k + α, n - k + β) + denominator = beta(α, β) + return choose * (numerator / denominator) + else + return 0.0 + end end function pdf(d::BetaBinomial) diff --git a/src/univariates.jl b/src/univariates.jl index 9b237c15d..d27eadf63 100644 --- a/src/univariates.jl +++ b/src/univariates.jl @@ -1,15 +1,27 @@ #### Domain && Support +comparator(::Type{Val{:closed}}) = <= +comparator(::Type{Val{:open}}) = < +comparator(s::Symbol) = comparator(Val{s}) immutable RealInterval lb::Float64 ub::Float64 + lc::Symbol + uc::Symbol - @compat RealInterval(lb::Real, ub::Real) = new(Float64(lb), Float64(ub)) + @compat RealInterval(lb::Real, ub::Real, lc::Symbol = :closed, uc::Symbol = :closed) = new(Float64(lb), Float64(ub), lc, uc) end minimum(r::RealInterval) = r.lb maximum(r::RealInterval) = r.ub -@compat in(x::Real, r::RealInterval) = (r.lb <= Float64(x) <= r.ub) + +lowerboundary(r::RealInterval) = r.lc +upperboundary(r::RealInterval) = r.uc + +lowercomparator(r::RealInterval) = comparator(r.lc) +uppercomparator(r::RealInterval) = comparator(r.uc) + +@compat in(x::Real, r::RealInterval) = lowercomparator(r)(r.lb,Float64(x)) && uppercomparator(r)(Float64(x),r.ub) @compat isbounded{D<:UnivariateDistribution}(d::Union{D,Type{D}}) = isupperbounded(d) && islowerbounded(d) @@ -19,6 +31,12 @@ maximum(r::RealInterval) = r.ub @compat hasfinitesupport{D<:DiscreteUnivariateDistribution}(d::Union{D,Type{D}}) = isbounded(d) @compat hasfinitesupport{D<:ContinuousUnivariateDistribution}(d::Union{D,Type{D}}) = false +@compat lowerboundary{D<:UnivariateDistribution}(d::Union{D,Type{D}}) = :closed +@compat upperboundary{D<:UnivariateDistribution}(d::Union{D,Type{D}}) = :closed + +@compat lowercomparator{D<:UnivariateDistribution}(d::Union{D,Type{D}}) = comparator(lowerboundary(d)) +@compat uppercomparator{D<:UnivariateDistribution}(d::Union{D,Type{D}}) = comparator(upperboundary(d)) + @compat function insupport!{D<:UnivariateDistribution}(r::AbstractArray, d::Union{D,Type{D}}, X::AbstractArray) length(r) == length(X) || throw(DimensionMismatch("Inconsistent array dimensions.")) @@ -31,11 +49,11 @@ end @compat insupport{D<:UnivariateDistribution}(d::Union{D,Type{D}}, X::AbstractArray) = insupport!(BitArray(size(X)), d, X) -@compat insupport{D<:ContinuousUnivariateDistribution}(d::Union{D,Type{D}},x::Real) = minimum(d) <= x <= maximum(d) -@compat insupport{D<:DiscreteUnivariateDistribution}(d::Union{D,Type{D}},x::Real) = isinteger(x) && minimum(d) <= x <= maximum(d) +@compat insupport{D<:ContinuousUnivariateDistribution}(d::Union{D,Type{D}},x::Real) = lowercomparator(d)(minimum(d),x) && uppercomparator(d)(x,maximum(d)) +@compat insupport{D<:DiscreteUnivariateDistribution}(d::Union{D,Type{D}},x::Real) = isinteger(x) && lowercomparator(d)(minimum(d),x) && uppercomparator(d)(x,maximum(d)) -@compat support{D<:ContinuousUnivariateDistribution}(d::Union{D,Type{D}}) = RealInterval(minimum(d), maximum(d)) -@compat support{D<:DiscreteUnivariateDistribution}(d::Union{D,Type{D}}) = round(Int, minimum(d)):round(Int, maximum(d)) +@compat support{D<:ContinuousUnivariateDistribution}(d::Union{D,Type{D}}) = RealInterval(minimum(d), maximum(d), lowerboundary(d), upperboundary(d)) +@compat support{D<:DiscreteUnivariateDistribution}(d::Union{D,Type{D}}) = (isfinite(minimum(d))?round(Int,minimum(d)):typemin(Int)):(isfinite(maximum(d))?round(Int, maximum(d)):typemax(Int)) # Type used for dispatch on finite support # T = true or false @@ -56,6 +74,15 @@ macro distr_support(D, lb, ub) end) end +macro distr_boundaries(D, l, u) + + @compat paramdecl = :(d::Union{$D, Type{$D}}) + + esc(quote + lowerboundary($(paramdecl)) = $l + upperboundary($(paramdecl)) = $u + end) +end ##### generic methods (fallback) ##### diff --git a/test/utils.jl b/test/utils.jl index 8a1bac71b..bb966a145 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -2,7 +2,20 @@ using Distributions using Base.Test # RealInterval -r = RealInterval(1.5, 4.0) -@test minimum(r) == 1.5 -@test maximum(r) == 4.0 +for (l,u,lc,uc) in [(:closed,:closed,<=,<=), + (:open,:closed,<,<=), + (:closed,:open,<=,<), + (:open,:open,<,<)] + + r = RealInterval(1.5, 4.0, l, u) + @test minimum(r) == 1.5 + @test maximum(r) == 4.0 + @test Distributions.lowerboundary(r) == l + @test Distributions.upperboundary(r) == u + @test Distributions.lowercomparator(r) == lc + @test Distributions.uppercomparator(r) == uc + +end + +@test RealInterval(0.0,1.0) == RealInterval(0.0,1.0,:closed,:closed)