diff --git a/README.md b/README.md index de3bc6c93a..a6ce1f77f0 100644 --- a/README.md +++ b/README.md @@ -11,7 +11,7 @@ Distributions.jl A Julia package for probability distributions and associated functions. Particularly, *Distributions* implements: * Moments (e.g mean, variance, skewness, and kurtosis), entropy, and other properties -* Probability density/mass functions (pdf) and their logarithm (logpdf) +* Probability density/mass functions (pdf/pmf) and their logarithm (logpdf/logpmf) * Moment generating functions and characteristic functions * Sampling from population or from a distribution * Maximum likelihood estimation diff --git a/docs/src/types.md b/docs/src/types.md index c00e172133..81ef6f2e48 100644 --- a/docs/src/types.md +++ b/docs/src/types.md @@ -41,8 +41,10 @@ The `ValueSupport` sub-types defined in `Distributions.jl` are: **Type** | **Element type** | **Descriptions** --- | --- | --- -`Discrete` | `Int` | Samples take discrete values -`Continuous` | `Float64` | Samples take continuous real values +`DiscreteSupport{T}` | `T` | Samples take any discrete values +`Discrete = DiscreteSupport{Int}` | `Int` | Samples take `Int` values +`ContinuousSupport{T <: Number}` | `T` | Samples take continuous values +`Continuous = ContinuousSupport{Float64}` | `Float64` | Samples take continuous `Float64` values Multiple samples are often organized into an array, depending on the variate form. @@ -69,22 +71,28 @@ abstract type Distribution{F<:VariateForm,S<:ValueSupport} <: Sampleable{F,S} en Distributions.Distribution ``` -To simplify the use in practice, we introduce a series of type alias as follows: +To simplify the use in practice, we introduce a series of type aliases as follows: ```julia const UnivariateDistribution{S<:ValueSupport} = Distribution{Univariate,S} const MultivariateDistribution{S<:ValueSupport} = Distribution{Multivariate,S} const MatrixDistribution{S<:ValueSupport} = Distribution{Matrixvariate,S} const NonMatrixDistribution = Union{UnivariateDistribution, MultivariateDistribution} -const DiscreteDistribution{F<:VariateForm} = Distribution{F,Discrete} -const ContinuousDistribution{F<:VariateForm} = Distribution{F,Continuous} +const CountableDistribution{F<:VariateForm, C<:CountableSupport} = Distribution{F,C} +const DiscreteDistribution{F<:VariateForm} = CountableDistribution{F,Discrete} +const ContinuousDistribution{F<:VariateForm} = Distribution{F,Continuous} -const DiscreteUnivariateDistribution = Distribution{Univariate, Discrete} -const ContinuousUnivariateDistribution = Distribution{Univariate, Continuous} -const DiscreteMultivariateDistribution = Distribution{Multivariate, Discrete} -const ContinuousMultivariateDistribution = Distribution{Multivariate, Continuous} -const DiscreteMatrixDistribution = Distribution{Matrixvariate, Discrete} -const ContinuousMatrixDistribution = Distribution{Matrixvariate, Continuous} +const CountableUnivariateDistribution{C<:CountableSupport} = UnivariateDistribution{C} +const DiscreteUnivariateDistribution = CountableUnivariateDistribution{Discrete} +const ContinuousUnivariateDistribution = UnivariateDistribution{Continuous} + +const CountableMultivariateDistribution{C<:CountableSupport} = MultivariateDistribution{C} +const DiscreteMultivariateDistribution = CountableMultivariateDistribution{Discrete} +const ContinuousMultivariateDistribution = MultivariateDistribution{Continuous} + +const CountableMatrixDistribution{C<:CountableSupport} = MatrixDistribution{C} +const DiscreteMatrixDistribution = CountableMatrixDistribution{Discrete} +const ContinuousMatrixDistribution = MatrixDistribution{Continuous} ``` All methods applicable to `Sampleable` also applies to `Distribution`. The API for distributions of different variate forms are different (refer to [univariates](@ref univariates), [multivariates](@ref multivariates), and [matrix](@ref matrix-variates) for details). diff --git a/src/Distributions.jl b/src/Distributions.jl index 242ab505ce..a4d57df046 100644 --- a/src/Distributions.jl +++ b/src/Distributions.jl @@ -30,20 +30,28 @@ export # generic types VariateForm, ValueSupport, + CountableSupport, + ContiguousSupport, + ContinuousSupport, + DiscontinuousSupport, + UnionSupport, Univariate, Multivariate, Matrixvariate, Discrete, Continuous, + Discontinuous, Sampleable, Distribution, UnivariateDistribution, MultivariateDistribution, MatrixDistribution, NoncentralHypergeometric, - NonMatrixDistribution, DiscreteDistribution, ContinuousDistribution, + CountableUnivariateDistribution, + CountableMultivariateDistribution, + CountableMatrixDistribution, DiscreteUnivariateDistribution, DiscreteMultivariateDistribution, DiscreteMatrixDistribution, @@ -69,9 +77,11 @@ export Chernoff, Chi, Chisq, + CompoundDistribution, Cosine, DiagNormal, DiagNormalCanon, + Dirac, Dirichlet, DirichletMultinomial, DiscreteUniform, @@ -93,6 +103,7 @@ export GeneralizedExtremeValue, Geometric, Gumbel, + Hurdle, Hypergeometric, InverseWishart, InverseGamma, @@ -138,6 +149,7 @@ export Rayleigh, Semicircle, Skellam, + SpikeSlab, StudentizedRange, SymTriangularDist, TDist, @@ -152,6 +164,7 @@ export WalleniusNoncentralHypergeometric, Weibull, Wishart, + ZeroInflated, ZeroMeanIsoNormal, ZeroMeanIsoNormalCanon, ZeroMeanDiagNormal, @@ -172,6 +185,8 @@ export components, # get components from a mixture model componentwise_pdf, # component-wise pdf for mixture models componentwise_logpdf, # component-wise logpdf for mixture models + componentwise_pmf, # component-wise pmf for mixture models + componentwise_logpmf, # component-wise logpmf for mixture models concentration, # the concentration parameter convolve, # convolve distributions of the same type dim, # sample dimension of multivariate distribution @@ -199,6 +214,8 @@ export loglikelihood, # log probability of array of IID draws logpdf, # log probability density logpdf!, # evaluate log pdf to provided storage + logpmf, # log probability mass + logpmf!, # evaluate log pmf to provided storage invscale, # Inverse scale parameter sqmahal, # squared Mahalanobis distance to Gaussian center @@ -221,7 +238,8 @@ export params, # get the tuple of parameters params!, # provide storage space to calculate the tuple of parameters for a multivariate distribution like mvlognormal partype, # returns a type large enough to hold all of a distribution's parameters' element types - pdf, # probability density function (ContinuousDistribution) + pdf, # probability density function (non-CountableSupport) + pmf, # probability mass function (non-ContinuousSupport) probs, # Get the vector of probabilities probval, # The pdf/pmf value for a uniform distribution quantile, # inverse of cdf (defined for p in (0,1)) @@ -277,6 +295,7 @@ include("qq.jl") include("estimators.jl") # mixture distributions (TODO: moveout) +include("mixtures/mixturedist.jl") include("mixtures/mixturemodel.jl") include("mixtures/unigmm.jl") @@ -290,6 +309,7 @@ API overview (major features): - `d = Dist(parameters...)` creates a distribution instance `d` for some distribution `Dist` (see choices below) with the specified `parameters` - `rand(d, sz)` samples from the distribution - `pdf(d, x)` and `logpdf(d, x)` compute the probability density or log-probability density of `d` at `x` +- `pmf(d, x)` and `logpmf(d, x)` compute the probability mass or log-probability mass of `d` at `x` - `cdf(d, x)` and `ccdf(d, x)` compute the (complementary) cumulative distribution function at `x` - `quantile(d, p)` is the inverse `cdf` (see also `cquantile`) - `mean(d)`, `var(d)`, `std(d)`, `skewness(d)`, `kurtosis(d)` compute moments of `d` @@ -303,23 +323,26 @@ information. Supported distributions: Arcsine, Bernoulli, Beta, BetaBinomial, BetaPrime, Binomial, Biweight, - Categorical, Cauchy, Chi, Chisq, Cosine, DiagNormal, DiagNormalCanon, + Categorical, Cauchy, Chi, Chisq, CompoundDistribution, + Cosine, DiagNormal, DiagNormalCanon, Dirichlet, DiscreteUniform, DoubleExponential, EdgeworthMean, EdgeworthSum, EdgeworthZ, Erlang, Epanechnikov, Exponential, FDist, FisherNoncentralHypergeometric, Frechet, FullNormal, FullNormalCanon, Gamma, GeneralizedPareto, - GeneralizedExtremeValue, Geometric, Gumbel, Hypergeometric, + GeneralizedExtremeValue, Geometric, Gumbel, Hurdle, Hypergeometric, InverseWishart, InverseGamma, InverseGaussian, IsoNormal, IsoNormalCanon, Kolmogorov, KSDist, KSOneSided, Laplace, Levy, - Logistic, LogNormal, MatrixBeta, MatrixFDist, MatrixNormal, MatrixTDist, MixtureModel, - Multinomial, MultivariateNormal, MvLogNormal, MvNormal, MvNormalCanon, - MvNormalKnownCov, MvTDist, NegativeBinomial, NoncentralBeta, NoncentralChisq, + Logistic, LogNormal, MatrixBeta, MatrixFDist, MatrixNormal, + MatrixTDist, MixtureModel, Multinomial, MultivariateNormal, MvLogNormal, + MvNormal, MvNormalCanon, MvNormalKnownCov, MvTDist, NegativeBinomial, + NoncentralBeta, NoncentralChisq, NoncentralF, NoncentralHypergeometric, NoncentralT, Normal, NormalCanon, - NormalInverseGaussian, Pareto, PGeneralizedGaussian, Poisson, PoissonBinomial, - QQPair, Rayleigh, Skellam, StudentizedRange, SymTriangularDist, TDist, TriangularDist, + NormalInverseGaussian, Pareto, PGeneralizedGaussian, + Poisson, PoissonBinomial, QQPair, Rayleigh, Skellam, SpikeSlab, + StudentizedRange, SymTriangularDist, TDist, TriangularDist, Triweight, Truncated, TruncatedNormal, Uniform, UnivariateGMM, VonMises, VonMisesFisher, WalleniusNoncentralHypergeometric, Weibull, - Wishart, ZeroMeanIsoNormal, ZeroMeanIsoNormalCanon, + Wishart, ZeroInflated, ZeroMeanIsoNormal, ZeroMeanIsoNormalCanon, ZeroMeanDiagNormal, ZeroMeanDiagNormalCanon, ZeroMeanFullNormal, ZeroMeanFullNormalCanon diff --git a/src/common.jl b/src/common.jl index 4f40d2abfb..37e3a78b2f 100644 --- a/src/common.jl +++ b/src/common.jl @@ -13,9 +13,21 @@ struct Matrixvariate <: VariateForm end `S <: ValueSupport` specifies the support of sample elements, either discrete or continuous. """ -abstract type ValueSupport end -struct Discrete <: ValueSupport end -struct Continuous <: ValueSupport end +abstract type ValueSupport{N} end +struct ContinuousSupport{N <: Number} <: ValueSupport{N} end +abstract type CountableSupport{C} <: ValueSupport{C} end +struct ContiguousSupport{C <: Integer} <: CountableSupport{C} end +struct UnionSupport{N1, N2, + S1 <: ValueSupport{N1}, + S2 <: ValueSupport{N2}} <: + ValueSupport{Union{N1, N2}} end + +const Discrete = ContiguousSupport{Int} +const Continuous = ContinuousSupport{Float64} +const DiscontinuousSupport{I, F} = + UnionSupport{I, F, <: CountableSupport{I}, + ContinuousSupport{F}} where {I <: Number, F <: Number} +const Discontinuous = DiscontinuousSupport{Int, Float64} ## Sampleable @@ -50,13 +62,14 @@ Base.size(s::Sampleable{Multivariate}) = (length(s),) """ eltype(s::Sampleable) + eltype(::ValueSupport) The default element type of a sample. This is the type of elements of the samples generated by the `rand` method. However, one can provide an array of different element types to store the samples using `rand!`. """ -Base.eltype(s::Sampleable{F,Discrete}) where {F} = Int -Base.eltype(s::Sampleable{F,Continuous}) where {F} = Float64 +Base.eltype(::Sampleable{F, <: ValueSupport{N}}) where {F, N} = N +Base.eltype(::ValueSupport{N}) where {N} = N """ nsamples(s::Sampleable) @@ -67,10 +80,11 @@ into an array, depending on the variate form. nsamples(t::Type{Sampleable}, x::Any) nsamples(::Type{D}, x::Number) where {D<:Sampleable{Univariate}} = 1 nsamples(::Type{D}, x::AbstractArray) where {D<:Sampleable{Univariate}} = length(x) -nsamples(::Type{D}, x::AbstractVector) where {D<:Sampleable{Multivariate}} = 1 +nsamples(::Type{D}, x::AbstractArray{<:AbstractVector}) where {D<:Sampleable{Multivariate}} = length(x) +nsamples(::Type{D}, x::AbstractVector{<:Number}) where {D<:Sampleable{Multivariate}} = 1 nsamples(::Type{D}, x::AbstractMatrix) where {D<:Sampleable{Multivariate}} = size(x, 2) -nsamples(::Type{D}, x::Number) where {D<:Sampleable{Matrixvariate}} = 1 -nsamples(::Type{D}, x::Array{Matrix{T}}) where {D<:Sampleable{Matrixvariate},T<:Number} = length(x) +nsamples(::Type{D}, x::AbstractMatrix{<:Number}) where {D<:Sampleable{Matrixvariate}} = 1 +nsamples(::Type{D}, x::AbstractArray{<:AbstractMatrix{T}}) where {D<:Sampleable{Matrixvariate},T<:Number} = length(x) """ Distribution{F<:VariateForm,S<:ValueSupport} <: Sampleable{F,S} @@ -85,23 +99,45 @@ abstract type Distribution{F<:VariateForm,S<:ValueSupport} <: Sampleable{F,S} en const UnivariateDistribution{S<:ValueSupport} = Distribution{Univariate,S} const MultivariateDistribution{S<:ValueSupport} = Distribution{Multivariate,S} const MatrixDistribution{S<:ValueSupport} = Distribution{Matrixvariate,S} -const NonMatrixDistribution = Union{UnivariateDistribution, MultivariateDistribution} -const DiscreteDistribution{F<:VariateForm} = Distribution{F,Discrete} +const CountableDistribution{F<:VariateForm, + C<:CountableSupport} = Distribution{F,C} +const DiscreteDistribution{F<:VariateForm} = CountableDistribution{F,Discrete} const ContinuousDistribution{F<:VariateForm} = Distribution{F,Continuous} -const DiscreteUnivariateDistribution = Distribution{Univariate, Discrete} -const ContinuousUnivariateDistribution = Distribution{Univariate, Continuous} -const DiscreteMultivariateDistribution = Distribution{Multivariate, Discrete} -const ContinuousMultivariateDistribution = Distribution{Multivariate, Continuous} -const DiscreteMatrixDistribution = Distribution{Matrixvariate, Discrete} -const ContinuousMatrixDistribution = Distribution{Matrixvariate, Continuous} - -variate_form(::Type{Distribution{VF,VS}}) where {VF<:VariateForm,VS<:ValueSupport} = VF -variate_form(::Type{T}) where {T<:Distribution} = variate_form(supertype(T)) - -value_support(::Type{Distribution{VF,VS}}) where {VF<:VariateForm,VS<:ValueSupport} = VS -value_support(::Type{T}) where {T<:Distribution} = value_support(supertype(T)) +const CountableUnivariateDistribution{C<:CountableSupport} = + UnivariateDistribution{C} +const DiscreteUnivariateDistribution = + CountableUnivariateDistribution{Discrete} +const ContinuousUnivariateDistribution = UnivariateDistribution{Continuous} +const CountableMultivariateDistribution{C<:CountableSupport} = + MultivariateDistribution{C} +const DiscreteMultivariateDistribution = + CountableMultivariateDistribution{Discrete} +const ContinuousMultivariateDistribution = MultivariateDistribution{Continuous} + +const CountableMatrixDistribution{C<:CountableSupport} = MatrixDistribution{C} +const DiscreteMatrixDistribution = CountableMatrixDistribution{Discrete} +const ContinuousMatrixDistribution = MatrixDistribution{Continuous} + +pdf(d::CountableDistribution) = pmf(d) +pdf(d::CountableDistribution, x) = pmf(d, x) +logpdf(d::CountableDistribution) = logpmf(d) +logpdf(d::CountableDistribution, x) = logpmf(d, x) + +const CountableMultivariateDistribution{C<:CountableSupport} = + MultivariateDistribution{C} +const DiscreteMultivariateDistribution = + CountableMultivariateDistribution{Discrete} +const ContinuousMultivariateDistribution = MultivariateDistribution{Continuous} + +const CountableMatrixDistribution{C<:CountableSupport} = MatrixDistribution{C} +const DiscreteMatrixDistribution = CountableMatrixDistribution{Discrete} +const ContinuousMatrixDistribution = MatrixDistribution{Continuous} + + +variate_form(::Type{<:Sampleable{VF, <:ValueSupport}}) where {VF<:VariateForm} = VF +value_support(::Type{<:Sampleable{<:VariateForm,VS}}) where {VS<:ValueSupport} = VS # allow broadcasting over distribution objects # to be decided: how to handle multivariate/matrixvariate distributions? diff --git a/src/edgeworth.jl b/src/edgeworth.jl index 1f23e568ee..2e85bd997e 100644 --- a/src/edgeworth.jl +++ b/src/edgeworth.jl @@ -55,7 +55,7 @@ end # Cornish-Fisher expansion. -function quantile(d::EdgeworthZ, p::Float64) +function quantile(d::EdgeworthZ, p::Real) s = skewness(d) k = kurtosis(d) z = quantile(Normal(0,1),p) @@ -63,7 +63,7 @@ function quantile(d::EdgeworthZ, p::Float64) z + s*(z2-1)/6.0 + k*z*(z2-3)/24.0 - s*s/36.0*z*(2.0*z2-5.0) end -function cquantile(d::EdgeworthZ, p::Float64) +function cquantile(d::EdgeworthZ, p::Real) s = skewness(d) k = kurtosis(d) z = cquantile(Normal(0,1),p) @@ -112,5 +112,5 @@ cdf(d::EdgeworthAbstract, x::Float64) = cdf(EdgeworthZ(d.dist,d.n), (x-mean(d))/ ccdf(d::EdgeworthAbstract, x::Float64) = ccdf(EdgeworthZ(d.dist,d.n), (x-mean(d))/std(d)) -quantile(d::EdgeworthAbstract, p::Float64) = mean(d) + std(d)*quantile(EdgeworthZ(d.dist,d.n), p) -cquantile(d::EdgeworthAbstract, p::Float64) = mean(d) + std(d)*cquantile(EdgeworthZ(d.dist,d.n), p) +quantile(d::EdgeworthAbstract, p::Real) = mean(d) + std(d)*quantile(EdgeworthZ(d.dist,d.n), p) +cquantile(d::EdgeworthAbstract, p::Real) = mean(d) + std(d)*cquantile(EdgeworthZ(d.dist,d.n), p) diff --git a/src/functionals.jl b/src/functionals.jl index 7f19b75b14..ef6626604e 100644 --- a/src/functionals.jl +++ b/src/functionals.jl @@ -18,6 +18,12 @@ function expectation(distr::DiscreteUnivariateDistribution, g::Function, epsilon sum(x -> f(x)*g(x), leftEnd:rightEnd) end +function expectation(distr::CountableUnivariateDistribution, + g::Function, epsilon::Real) + f = x->pdf(distr,x) + sum(x -> f(x)*g(x), support(distr)) +end + function expectation(distr::UnivariateDistribution, g::Function) expectation(distr, g, 1e-10) end diff --git a/src/mixtures/mixturedist.jl b/src/mixtures/mixturedist.jl new file mode 100644 index 0000000000..b56e9dfbe9 --- /dev/null +++ b/src/mixtures/mixturedist.jl @@ -0,0 +1,345 @@ +# General mixture distributions +""" + + All subtypes of `AbstractMixtureDistribution` must implement the following method: + + - components(d): return the distribution components + + - probs(d): return a vector of prior probabilities over components. + + And may implement these: + + - ncomponents(d): the number of components + + - component(d, k): return the k-th component + +Mixture distributions can contain any combination of distributions with countable and uncountable support. All of the distributions must have the same VariateForm, however. +""" +abstract type AbstractMixtureDistribution{VF<:VariateForm,VS<:ValueSupport} <: Distribution{VF, VS} end + +""" + components(d::AbstractMixtureDistribution) + +Get a list of components of the mixture distribution `d`. +""" +components(d::AbstractMixtureDistribution) + +""" + components(d::AbstractMixtureDistribution) + +Get the number of components of the mixture distribution `d`. +""" +ncomponents(d::AbstractMixtureDistribution) = length(components(d)) + +""" + component(d::AbstractMixtureDistribution, k::Integer) + +Get the `k'th component of the mixture distribution `d`. +""" +component(d::AbstractMixtureDistribution, k::Integer) = components(d)[k] + +""" + probs(d::AbstractMixtureDistribution) + +Get the vector of prior probabilities of all components of `d`. +""" +probs(d::AbstractMixtureDistribution) + +const UnivariateMixtureDistribution{S<:ValueSupport} = + AbstractMixtureDistribution{Univariate,S} +const MultivariateMixtureDistribution{S<:ValueSupport} = + AbstractMixtureDistribution{Multivariate,S} +const MatrixvariateMixtureDistribution{S<:ValueSupport} = + AbstractMixtureDistribution{Matrixvariate,S} + +minimum(md::AbstractMixtureDistribution) = minimum(minimum.(components(md))) +maximum(md::AbstractMixtureDistribution) = maximum(maximum.(components(md))) + +function mean(d::UnivariateMixtureDistribution) + K = ncomponents(d) + p = probs(d) + m = 0.0 + for i = 1:K + pi = p[i] + if pi > 0.0 + c = component(d, i) + m += mean(c) * pi + end + end + return m +end + +function mean(d::MultivariateMixtureDistribution) + K = ncomponents(d) + p = probs(d) + m = zeros(length(d)) + for i = 1:K + pi = p[i] + if pi > 0.0 + c = component(d, i) + BLAS.axpy!(pi, mean(c), m) + end + end + return m +end + +function var(d::UnivariateMixtureDistribution) + K = ncomponents(d) + p = probs(d) + means = Vector{Float64}(undef, K) + m = 0.0 + v = 0.0 + for i = 1:K + pi = p[i] + if pi > 0.0 + ci = component(d, i) + means[i] = mi = mean(ci) + m += pi * mi + v += pi * var(ci) + end + end + for i = 1:K + pi = p[i] + if pi > 0.0 + v += pi * abs2(means[i] - m) + end + end + return v +end + +function cov(d::MultivariateMixtureDistribution) + K = ncomponents(d) + p = probs(d) + m = zeros(length(d)) + md = zeros(length(d)) + V = zeros(length(d),length(d)) + + for i = 1:K + pi = p[i] + if pi > 0.0 + c = component(d, i) + BLAS.axpy!(pi, mean(c), m) + BLAS.axpy!(pi, cov(c), V) + end + end + for i = 1:K + pi = p[i] + if pi > 0.0 + c = component(d, i) + # todo: use more in-place operations + md = mean(c) - m + BLAS.axpy!(pi, md*md', V) + end + end + return V +end + +#### show + +function show(io::IO, d::D) where D <: AbstractMixtureDistribution + K = ncomponents(d) + pr = probs(d) + println(io, "$D(K = $K)") + Ks = min(K, 8) + for i = 1:Ks + @printf(io, "components[%d] (prior = %.4f): ", i, pr[i]) + println(io, component(d, i)) + end + if Ks < K + println(io, "The rest are omitted ...") + end +end + +""" + + A SpikeSlab distribution contains a Dirac delta and a continuous distribution. +""" +struct SpikeSlab{N<:Number, + C<:Distribution{Univariate, ContinuousSupport{N}}} <: + AbstractMixtureDistribution{Univariate, + DiscontinuousSupport{N, N}} + pSpike::Float64 + spike::N + slab::C +end + +ncomponents(d::SpikeSlab) = 2 +components(ss::SpikeSlab) = [Dirac(ss.spike), ss.slab] +component(ss::SpikeSlab, k::Integer) = k == 1 ? Dirac(ss.spike) : + (k==2 ? ss.slab : throw(BoundsError())) +probs(ss::SpikeSlab) = [ss.pSpike, one(ss.pSpike) - ss.pSpike] +minimum(ss::SpikeSlab) = min(ss.spike, minimum(ss.slab)) +maximum(ss::SpikeSlab) = max(ss.spike, maximum(ss.slab)) +mean(ss::SpikeSlab) = ss.spike * ss.pSpike + + mean(ss.slab) * (one(ss.pSpike) - ss.pSpike) +var(ss::SpikeSlab) = + (one(ss.pSpike) - ss.pSpike) * (var(ss.slab) + mean(ss.slab)^2) + + ss.pSpike * ss.spike^2 - mean(ss)^2 +# Probability mass always beats density! +mode(ss::SpikeSlab) = iszero(ss.pSpike) ? mode(ss.slab) : ss.spike +pdf(ss::SpikeSlab, x) = (x ≈ ss.spike) ? Inf : pdf(ss.slab, x) +pmf(ss::SpikeSlab, x) = + (iszero(ss.pSpike) || x ≉ ss.spike) ? 0.0 : ss.pSpike +logpmf(ss::SpikeSlab, x) = + iszero(ss.pSpike) || x ≉ ss.spike ? -Inf : log(ss.pSpike) +logpdf(ss::SpikeSlab, x) = (x ≈ ss.spike) ? Inf : logpdf(ss.slab, x) +cdf(ss::SpikeSlab, x) = + cdf(ss.slab, x) + (x < ss.spike ? zero(ss.pSpike) : ss.pSpike) +rand(rng::AbstractRNG, ss::SpikeSlab) = + rand(rng) < ss.sSpike ? ss.spike : rand(rng, ss.slab) +function quantile(ss::SpikeSlab, p::Real) + toSpike = cdf(ss.slab, ss.spike) * (one(ss.pSpike) - ss.pSpike) + if p ≤ toSpike + quantile(ss.slab, p / (one(ss.pSpike) - ss.pSpike)) + elseif p ≤ toSpike + ss.pSpike + ss.spike + else + quantile(ss.slab, (p - ss.pSpike) / (one(ss.pSpike) - ss.pSpike)) + end +end +insupport(ss::SpikeSlab, x::Number) = insupport(ss.slab, x) || x ≈ ss.spike + +""" + + A `CompoundDistribution` contains multiple discrete and continuous distributions. +""" +struct CompoundDistribution{TR<:Tuple, TD<:Tuple, NR<:Number, ND<:Number, + Counts<:Union{CountableUnivariateDistribution{ND}, + Nothing}} <: + AbstractMixtureDistribution{Univariate, + UnionSupport{NR, ND, + ContinuousSupport{NR}, + CountableSupport{ND}}} + choice::Categorical{Float64, Vector{Float64}, Int} + countable::Counts + nContinuous::Int + pContinuous::Float64 + continuous::TR + discrete::TD + + function CompoundDistribution{TR, Tuple{}}(ps::Vector{Float64}, + continuous::TR, ::Tuple{}, + ::NoArgCheck) where TR + return new{TR, Tuple{}, eltype(first(TR)), eltype(first(TR)), + Nothing}(Categorical(ps), nothing, + length(continuous), 1.0, continuous, ()) + end + + function CompoundDistribution{TR, Tuple{Counts}}(ps::Vector{Float64}, + continuous::TR, + discrete::Tuple{Counts}, + ::NoArgCheck) where + {TR, Counts} + n = length(continuous) + return new{TR, Tuple{Counts}, eltype(first(TR)), Int, + Counts}(Categorical(ps), first(discrete), + n, sum(probs(choice)[1:n]), + continuous, discrete) + end + + function CompoundDistribution{TR, TD}(ps::Vector{Float64}, + continuous::TR, + discrete::TD, + ::NoArgCheck) where {TR, TD} + n = length(continuous) + dict = Dict{ND, Float64}() + for d in discrete + for s in support(d) + dict[s] = get(dict, s, 0.0) + end + end + choice = DiscreteNonParametric(dict) + DR = length(continuous) == 0 ? Float64 : eltype(first(TR)) + return new{TR, TD, DR, eltype(first(TD)), + typeof(choice)}(Categorical(ps), choice, + n, sum(probs(choice)[1:n]), + continuous, discrete) + end + + function CompoundDistribution{TR, TD}(ps::Vector{Float64}, + continuous::TR, + discrete::TD) where {TR <: Tuple, + TD <: Tuple} + @check_args(CompoundDistribution, + length(support(choice)) == + length(continuous) + length(discrete)) + if length(continuous) > 0 + vs = Set(value_support.(continuous)) + @check_args(CompoundDistribution, length(vs) == 1) + @check_args(CompoundDistribution, eltype(first(vs)) ≡ ND) + @check_args(CompoundDistribution, + length(Set(variate_form.(continuous))) == 1) + end + if length(discrete) > 0 + @check_args(CompoundDistribution, all(hasfinitesupport.(discrete))) + ds = Set(map(value_support, discrete)) + @check_args(CompoundDistribution, length(ds) == 1) + @check_args(CompoundDistribution, eltype(first(ds)) ≡ NR) + @check_args(CompoundDistribution, + length(Set(variate_form.(discrete))) == 1) + end + return CompoundDistribution{TR, TD}(ps, continuous, discrete, + NoArgCheck()) + end +end + +CompoundDistribution(ps::Vector{Float64}, continuous::TR, discrete::TD) where +{TR <: Tuple, TD <: Tuple} = + CompoundDistribution{TR, TD}(ps, continuous, discrete) +function CompoundDistribution(ps::Vector{Float64}, continuous, discrete) + tcont = tuple(continuous...) + tdisc = tuple(discrete...) + return CompoundDistribution{typeof(tcont), typeof(tdisc)}(ps, tcont, tdisc) +end + +function CompoundDistribution(ps::Vector{Float64}, distributions) + continuous = filter(d -> value_support(d) <: ContinuousSupport, + collect(distributions)) + discrete = filter(d -> value_support(d) <: CountableSupport, + collect(distributions)) + return CompoundDistribution(ps, continuous, discrete) +end + +CompoundDistribution(dict::Dict) = CompoundDistribution(values(dict), + keys(dict)) + +function ZeroInflated(pZero::Float64, dist::D) where D <: UnivariateDistribution + z = zero(eltype(D)) + insupport(dist, z) || + throw(DomainError(z, + "ZeroInflated distribution must have $z in support.")) + return CompoundDistribution([pZero, 1.0 - pZero], + [Dirac(z), dist]) +end + +function Hurdle(pVal::Float64, dist::D; val = zero(eltype(D))) where + D <: UnivariateDistribution + insupport(dist, val) && + throw(DomainError(val, + "Hurdle distribution must not have $val in support.")) + return CompoundDistribution([pVal, 1.0 - pVal], + [Dirac(val), dist]) +end + +components(cd::CompoundDistribution) = [cd.continuous..., cd.discrete...] +probs(cd::CompoundDistribution) = probs(cd.choice) +# Probability mass always beats density! +pdf(cd::CompoundDistribution, x) = + (cd.countable ≢ nothing && any(x .≈ support(cd.countable))) ? Inf : + reduce(+, pdf.(cd.continuous, x)) * cd.pContinuous + +pmf(cd::CompoundDistribution, x) = isnothing(cd.countable) ? 0.0 : + (pmf(cd.countable) * (1.0 - cd.pContinuous)) +function cdf(cd::CompoundDistribution, x) + cdf_point = cdf(cd.countable, x) * (1.0 - cd.pContinuous) + cs = cd.continuous + ps = probs(cd.choice) + for idx in Base.OneTo(cd.nContinuous) + cdf_point += ps[idx] * cdf(cs[idx], x) + end + return cdf_point +end +rand(rng::AbstractRNG, cd::CompoundDistribution) = + rand(rng, component(cd, rand(rng, cd.choice))) +insupport(cd::CompoundDistribution, x::Number) = + any(insupport.(cd.continuous, x)) || any(insupport.(cd.discrete, x)) diff --git a/src/mixtures/mixturemodel.jl b/src/mixtures/mixturemodel.jl index 5efdb01eaf..fe187081c9 100644 --- a/src/mixtures/mixturemodel.jl +++ b/src/mixtures/mixturemodel.jl @@ -2,15 +2,21 @@ """ - All subtypes of `AbstractMixtureModel` should implement the following methods: + All subtypes of `AbstractMixtureModel` must implement the following methods: + + - components(d): return the distribution components + + - probs(d): return a vector of prior probabilities over components. + + And may implement these: - ncomponents(d): the number of components - component(d, k): return the k-th component - - probs(d): return a vector of prior probabilities over components. +Mixture models are subtypes of AbstractMixtureDistributions, where all of the distributions being mixed are of the same form (though their parameters may be different). """ -abstract type AbstractMixtureModel{VF<:VariateForm,VS<:ValueSupport,C<:Distribution} <: Distribution{VF, VS} end +abstract type AbstractMixtureModel{VF<:VariateForm,VS<:ValueSupport,C<:Distribution} <: AbstractMixtureDistribution{VF, VS} end """ MixtureModel{VF<:VariateForm,VS<:ValueSupport,C<:Distribution,CT<:Real} @@ -43,20 +49,6 @@ The type of the components of `d`. """ component_type(d::AbstractMixtureModel{VF,VS,C}) where {VF,VS,C} = C -""" - components(d::AbstractMixtureModel) - -Get a list of components of the mixture model `d`. -""" -components(d::AbstractMixtureModel) - -""" - probs(d::AbstractMixtureModel) - -Get the vector of prior probabilities of all components of `d`. -""" -probs(d::AbstractMixtureModel) - """ mean(d::Union{UnivariateMixture, MultivariateMixture}) @@ -125,7 +117,7 @@ the components given by ``params``, and a prior probability vector. If no `prior` is provided then all components will have the same prior probabilities. """ function MixtureModel(::Type{C}, params::AbstractArray) where C<:Distribution - components = C[_construct_component(C, a) for a in params] + components = [_construct_component(C, a) for a in params] MixtureModel(components) end @@ -142,7 +134,7 @@ _construct_component(::Type{C}, arg) where {C<:Distribution} = C(arg) _construct_component(::Type{C}, args::Tuple) where {C<:Distribution} = C(args...) function MixtureModel(::Type{C}, params::AbstractArray, p::Vector{T}) where {C<:Distribution,T<:Real} - components = C[_construct_component(C, a) for a in params] + components = [_construct_component(C, a) for a in params] MixtureModel(components, p) end @@ -167,94 +159,6 @@ probs(d::MixtureModel) = probs(d.prior) params(d::MixtureModel) = ([params(c) for c in d.components], params(d.prior)[1]) partype(d::MixtureModel) = promote_type(partype(d.prior), map(partype, d.components)...) -minimum(d::MixtureModel) = minimum([minimum(dci) for dci in d.components]) -maximum(d::MixtureModel) = maximum([maximum(dci) for dci in d.components]) - -function mean(d::UnivariateMixture) - K = ncomponents(d) - p = probs(d) - m = 0.0 - for i = 1:K - pi = p[i] - if pi > 0.0 - c = component(d, i) - m += mean(c) * pi - end - end - return m -end - -function mean(d::MultivariateMixture) - K = ncomponents(d) - p = probs(d) - m = zeros(length(d)) - for i = 1:K - pi = p[i] - if pi > 0.0 - c = component(d, i) - BLAS.axpy!(pi, mean(c), m) - end - end - return m -end - -""" - var(d::UnivariateMixture) - -Compute the overall variance (only for ``UnivariateMixture``). -""" -function var(d::UnivariateMixture) - K = ncomponents(d) - p = probs(d) - means = Vector{Float64}(undef, K) - m = 0.0 - v = 0.0 - for i = 1:K - pi = p[i] - if pi > 0.0 - ci = component(d, i) - means[i] = mi = mean(ci) - m += pi * mi - v += pi * var(ci) - end - end - for i = 1:K - pi = p[i] - if pi > 0.0 - v += pi * abs2(means[i] - m) - end - end - return v -end - -function cov(d::MultivariateMixture) - K = ncomponents(d) - p = probs(d) - m = zeros(length(d)) - md = zeros(length(d)) - V = zeros(length(d),length(d)) - - for i = 1:K - pi = p[i] - if pi > 0.0 - c = component(d, i) - BLAS.axpy!(pi, mean(c), m) - BLAS.axpy!(pi, cov(c), V) - end - end - for i = 1:K - pi = p[i] - if pi > 0.0 - c = component(d, i) - # todo: use more in-place operations - md = mean(c) - m - BLAS.axpy!(pi, md*md', V) - end - end - return V -end - - #### show @@ -272,8 +176,6 @@ function show(io::IO, d::MixtureModel) end end - - #### Evaluation function insupport(d::AbstractMixtureModel, x::AbstractVector) diff --git a/src/multivariate/dirichlet.jl b/src/multivariate/dirichlet.jl index 9052d04e1f..07d74522c6 100644 --- a/src/multivariate/dirichlet.jl +++ b/src/multivariate/dirichlet.jl @@ -20,7 +20,7 @@ Dirichlet(alpha) # Dirichlet distribution with parameter vector alpha Dirichlet(k, a) # Dirichlet distribution with parameter a * ones(k) ``` """ -struct Dirichlet{T<:Real} <: ContinuousMultivariateDistribution +struct Dirichlet{T<:Real} <: MultivariateDistribution{ContinuousSupport{T}} alpha::Vector{T} alpha0::T lmnB::T @@ -57,7 +57,6 @@ end length(d::DirichletCanon) = length(d.alpha) -eltype(::Dirichlet{T}) where {T} = T #### Conversions convert(::Type{Dirichlet{Float64}}, cf::DirichletCanon) = Dirichlet(cf.alpha) convert(::Type{Dirichlet{T}}, alpha::Vector{S}) where {T<:Real, S<:Real} = diff --git a/src/multivariate/dirichletmultinomial.jl b/src/multivariate/dirichletmultinomial.jl index caaadd80af..2167e4e529 100644 --- a/src/multivariate/dirichletmultinomial.jl +++ b/src/multivariate/dirichletmultinomial.jl @@ -54,7 +54,7 @@ function insupport(d::DirichletMultinomial, x::AbstractVector{T}) where T<:Real end return sum(x) == ntrials(d) end -function _logpdf(d::DirichletMultinomial{S}, x::AbstractVector{T}) where {T<:Real, S<:Real} +function _logpmf(d::DirichletMultinomial{S}, x::AbstractVector{T}) where {T<:Real, S<:Real} c = lgamma(S(d.n + 1)) + lgamma(d.α0) - lgamma(d.n + d.α0) for j in eachindex(x) @inbounds xj, αj = x[j], d.α[j] diff --git a/src/multivariate/multinomial.jl b/src/multivariate/multinomial.jl index d69659c587..208559dcdc 100644 --- a/src/multivariate/multinomial.jl +++ b/src/multivariate/multinomial.jl @@ -117,7 +117,7 @@ function entropy(d::Multinomial) for pr in p b = Binomial(n, pr) for x in 0:n - s += pdf(b, x) * lgamma(x+1) + s += pmf(b, x) * lgamma(x+1) end end return s @@ -140,7 +140,7 @@ function insupport(d::Multinomial, x::AbstractVector{T}) where T<:Real return s == ntrials(d) # integer computation would not yield truncation errors end -function _logpdf(d::Multinomial, x::AbstractVector{T}) where T<:Real +function _logpmf(d::Multinomial, x::AbstractVector{T}) where T<:Real p = probs(d) n = ntrials(d) S = eltype(p) diff --git a/src/multivariate/mvlognormal.jl b/src/multivariate/mvlognormal.jl index bfb4b1a7d5..ce3179c21f 100644 --- a/src/multivariate/mvlognormal.jl +++ b/src/multivariate/mvlognormal.jl @@ -23,7 +23,7 @@ # ########################################################### -abstract type AbstractMvLogNormal <: ContinuousMultivariateDistribution end +abstract type AbstractMvLogNormal{T} <: MultivariateDistribution{ContinuousSupport{T}} end function insupport(::Type{D},x::AbstractVector{T}) where {T<:Real,D<:AbstractMvLogNormal} for i=1:length(x) @@ -161,7 +161,7 @@ Mean vector ``\\boldsymbol{\\mu}`` and covariance matrix ``\\boldsymbol{\\Sigma} underlying normal distribution are known as the *location* and *scale* parameters of the corresponding lognormal distribution. """ -struct MvLogNormal{T<:Real,Cov<:AbstractPDMat,Mean<:AbstractVector} <: AbstractMvLogNormal +struct MvLogNormal{T<:Real,Cov<:AbstractPDMat,Mean<:AbstractVector} <: AbstractMvLogNormal{T} normal::MvNormal{T,Cov,Mean} end @@ -176,7 +176,6 @@ MvLogNormal(σ::AbstractVector) = MvLogNormal(MvNormal(σ)) MvLogNormal(d::Int,s::Real) = MvLogNormal(MvNormal(d,s)) -eltype(::MvLogNormal{T}) where {T} = T ### Conversion function convert(::Type{MvLogNormal{T}}, d::MvLogNormal) where T<:Real MvLogNormal(convert(MvNormal{T}, d.normal)) diff --git a/src/multivariate/mvnormal.jl b/src/multivariate/mvnormal.jl index 6effd498f8..67f1ee36c1 100644 --- a/src/multivariate/mvnormal.jl +++ b/src/multivariate/mvnormal.jl @@ -67,7 +67,7 @@ const ZeroMeanDiagNormal = MvNormal{PDiagMat, ZeroVector{Float64}} const ZeroMeanFullNormal = MvNormal{PDMat, ZeroVector{Float64}} ``` """ -abstract type AbstractMvNormal <: ContinuousMultivariateDistribution end +abstract type AbstractMvNormal{T<:Real} <: MultivariateDistribution{ContinuousSupport{T}} end ### Generic methods (for all AbstractMvNormal subtypes) @@ -170,7 +170,7 @@ isotropic covariance as `abs2(sig) * eye(d)`. **Note:** The constructor will choose an appropriate covariance form internally, so that special structure of the covariance can be exploited. """ -struct MvNormal{T<:Real,Cov<:AbstractPDMat,Mean<:AbstractVector} <: AbstractMvNormal +struct MvNormal{T<:Real,Cov<:AbstractPDMat,Mean<:AbstractVector} <: AbstractMvNormal{T} μ::Mean Σ::Cov end @@ -219,8 +219,6 @@ MvNormal(Σ::Matrix{<:Real}) = MvNormal(PDMat(Σ)) MvNormal(σ::Vector{<:Real}) = MvNormal(PDiagMat(abs2.(σ))) MvNormal(d::Int, σ::Real) = MvNormal(ScalMat(d, abs2(σ))) - -eltype(::MvNormal{T}) where {T} = T ### Conversion function convert(::Type{MvNormal{T}}, d::MvNormal) where T<:Real MvNormal(convert(AbstractArray{T}, d.μ), convert(AbstractArray{T}, d.Σ)) diff --git a/src/multivariate/mvnormalcanon.jl b/src/multivariate/mvnormalcanon.jl index 23b3b00d93..70933a186c 100644 --- a/src/multivariate/mvnormalcanon.jl +++ b/src/multivariate/mvnormalcanon.jl @@ -63,7 +63,7 @@ Construct a multivariate normal distribution of dimension `d`, with zero mean an **Note:** `MvNormalCanon` share the same set of methods as `MvNormal`. """ -struct MvNormalCanon{T<:Real,P<:AbstractPDMat,V<:AbstractVector} <: AbstractMvNormal +struct MvNormalCanon{T<:Real,P<:AbstractPDMat,V<:AbstractVector} <: AbstractMvNormal{T} μ::V # the mean vector h::V # potential vector, i.e. inv(Σ) * μ J::P # precision matrix, i.e. inv(Σ) @@ -156,7 +156,6 @@ length(d::MvNormalCanon) = length(d.μ) mean(d::MvNormalCanon) = convert(Vector{eltype(d.μ)}, d.μ) params(d::MvNormalCanon) = (d.μ, d.h, d.J) @inline partype(d::MvNormalCanon{T}) where {T<:Real} = T -eltype(::MvNormalCanon{T}) where {T} = T var(d::MvNormalCanon) = diag(inv(d.J)) cov(d::MvNormalCanon) = Matrix(inv(d.J)) diff --git a/src/multivariate/mvtdist.jl b/src/multivariate/mvtdist.jl index dc031c8627..45ceff17e4 100644 --- a/src/multivariate/mvtdist.jl +++ b/src/multivariate/mvtdist.jl @@ -2,9 +2,9 @@ ## Generic multivariate t-distribution class -abstract type AbstractMvTDist <: ContinuousMultivariateDistribution end +abstract type AbstractMvTDist{T} <: MultivariateDistribution{ContinuousSupport{T}} end -struct GenericMvTDist{T<:Real, Cov<:AbstractPDMat, Mean<:AbstractVector} <: AbstractMvTDist +struct GenericMvTDist{T<:Real, Cov<:AbstractPDMat, Mean<:AbstractVector} <: AbstractMvTDist{T} df::T # non-integer degrees of freedom allowed dim::Int zeromean::Bool @@ -103,7 +103,6 @@ logdet_cov(d::GenericMvTDist) = d.df>2 ? logdet((d.df/(d.df-2))*d.Σ) : NaN params(d::GenericMvTDist) = (d.df, d.μ, d.Σ) @inline partype(d::GenericMvTDist{T}) where {T} = T -eltype(::GenericMvTDist{T}) where {T} = T # For entropy calculations see "Multivariate t Distributions and their Applications", S. Kotz & S. Nadarajah function entropy(d::GenericMvTDist) diff --git a/src/multivariates.jl b/src/multivariates.jl index d2bd273440..babec8384e 100644 --- a/src/multivariates.jl +++ b/src/multivariates.jl @@ -170,7 +170,8 @@ end # pdf and logpdf """ - pdf(d::MultivariateDistribution, x::AbstractArray) + pdf(d::ContinuousMultivariateDistribution, x::AbstractArray) + pmf(d::CountableMultivariateDistribution, x::AbstractArray) Return the probability density of distribution `d` evaluated at `x`. @@ -180,10 +181,12 @@ to `x[:,i]` (i.e. treating each column as a sample). `pdf!(r, d, x)` will write the results to a pre-allocated array `r`. """ -pdf(d::MultivariateDistribution, x::AbstractArray) +pdf(d::ContinuousMultivariateDistribution, x::AbstractArray) +# pmf(d::CountableMultivariateDistribution, x::AbstractArray) """ - logpdf(d::MultivariateDistribution, x::AbstractArray) + logpdf(d::ContinuousMultivariateDistribution, x::AbstractArray) + logpmf(d::CountableMultivariateDistribution, x::AbstractArray) Return the logarithm of probability density evaluated at `x`. @@ -192,69 +195,134 @@ Return the logarithm of probability density evaluated at `x`. `logpdf!(r, d, x)` will write the results to a pre-allocated array `r`. """ -logpdf(d::MultivariateDistribution, x::AbstractArray) +logpdf(d::ContinuousMultivariateDistribution, x::AbstractArray) +# logpmf(d::CountableMultivariateDistribution, x::AbstractArray) -_pdf(d::MultivariateDistribution, X::AbstractVector) = exp(_logpdf(d, X)) +_pdf(d::ContinuousMultivariateDistribution, X::AbstractVector) = + exp(_logpdf(d, X)) +_pmf(d::CountableMultivariateDistribution, X::AbstractVector) = + exp(_logpmf(d, X)) -function logpdf(d::MultivariateDistribution, X::AbstractVector) +function logpdf(d::ContinuousMultivariateDistribution, X::AbstractVector) length(X) == length(d) || throw(DimensionMismatch("Inconsistent array dimensions.")) _logpdf(d, X) end +function logpmf(d::CountableMultivariateDistribution, X::AbstractVector) + length(X) == length(d) || + throw(DimensionMismatch("Inconsistent array dimensions.")) + _logpmf(d, X) +end -function pdf(d::MultivariateDistribution, X::AbstractVector) +function pdf(d::ContinuousMultivariateDistribution, X::AbstractVector) length(X) == length(d) || throw(DimensionMismatch("Inconsistent array dimensions.")) _pdf(d, X) end -function _logpdf!(r::AbstractArray, d::MultivariateDistribution, X::AbstractMatrix) +function pmf(d::CountableMultivariateDistribution, X::AbstractVector) + length(X) == length(d) || + throw(DimensionMismatch("Inconsistent array dimensions.")) + _pmf(d, X) +end + +function _logpdf!(r::AbstractArray, d::ContinuousMultivariateDistribution, + X::AbstractMatrix) for i in 1 : size(X,2) @inbounds r[i] = logpdf(d, view(X,:,i)) end return r end -function _pdf!(r::AbstractArray, d::MultivariateDistribution, X::AbstractMatrix) +function _logpmf!(r::AbstractArray, d::CountableMultivariateDistribution, + X::AbstractMatrix) + for i in 1 : size(X,2) + @inbounds r[i] = logpmf(d, view(X,:,i)) + end + return r +end + +function _pdf!(r::AbstractArray, d::ContinuousMultivariateDistribution, + X::AbstractMatrix) for i in 1 : size(X,2) @inbounds r[i] = pdf(d, view(X,:,i)) end return r end -function logpdf!(r::AbstractArray, d::MultivariateDistribution, X::AbstractMatrix) +function _pmf!(r::AbstractArray, d::CountableMultivariateDistribution, + X::AbstractMatrix) + for i in 1 : size(X,2) + @inbounds r[i] = pmf(d, view(X,:,i)) + end + return r +end + +function logpdf!(r::AbstractArray, d::ContinuousMultivariateDistribution, + X::AbstractMatrix) size(X) == (length(d), length(r)) || throw(DimensionMismatch("Inconsistent array dimensions.")) _logpdf!(r, d, X) end -function pdf!(r::AbstractArray, d::MultivariateDistribution, X::AbstractMatrix) +function logpmf!(r::AbstractArray, d::CountableMultivariateDistribution, + X::AbstractMatrix) + size(X) == (length(d), length(r)) || + throw(DimensionMismatch("Inconsistent array dimensions.")) + _logpmf!(r, d, X) +end + +function pdf!(r::AbstractArray, d::ContinuousMultivariateDistribution, + X::AbstractMatrix) size(X) == (length(d), length(r)) || throw(DimensionMismatch("Inconsistent array dimensions.")) _pdf!(r, d, X) end -function logpdf(d::MultivariateDistribution, X::AbstractMatrix) +function pmf!(r::AbstractArray, d::CountableMultivariateDistribution, + X::AbstractMatrix) + size(X) == (length(d), length(r)) || + throw(DimensionMismatch("Inconsistent array dimensions.")) + _pmf!(r, d, X) +end + +function logpdf(d::ContinuousMultivariateDistribution, X::AbstractMatrix) size(X, 1) == length(d) || throw(DimensionMismatch("Inconsistent array dimensions.")) T = promote_type(partype(d), eltype(X)) _logpdf!(Vector{T}(undef, size(X,2)), d, X) end -function pdf(d::MultivariateDistribution, X::AbstractMatrix) +function logpmf(d::CountableMultivariateDistribution, X::AbstractMatrix) + size(X, 1) == length(d) || + throw(DimensionMismatch("Inconsistent array dimensions.")) + T = promote_type(partype(d), eltype(X)) + _logpmf!(Vector{T}(undef, size(X,2)), d, X) +end + +function pdf(d::ContinuousMultivariateDistribution, X::AbstractMatrix) size(X, 1) == length(d) || throw(DimensionMismatch("Inconsistent array dimensions.")) T = promote_type(partype(d), eltype(X)) _pdf!(Vector{T}(undef, size(X,2)), d, X) end +function pmf(d::CountableMultivariateDistribution, X::AbstractMatrix) + size(X, 1) == length(d) || + throw(DimensionMismatch("Inconsistent array dimensions.")) + T = promote_type(partype(d), eltype(X)) + _pmf!(Vector{T}(undef, size(X,2)), d, X) +end + """ - _logpdf{T<:Real}(d::MultivariateDistribution, x::AbstractArray) + _logpdf{T<:Real}(d::ContinuousMultivariateDistribution, x::AbstractArray) + _logpmf{T<:Real}(d::CountableMultivariateDistribution, x::AbstractArray) -Evaluate logarithm of pdf value for a given vector `x`. This function need not perform dimension checking. +Evaluate logarithm of pdf/pmf value for a given vector `x`. This function need not perform dimension checking. Generally, one does not need to implement `pdf` (or `_pdf`) as fallback methods are provided in `src/multivariates.jl`. """ -_logpdf(d::MultivariateDistribution, x::AbstractArray) +_logpdf(d::ContinuousMultivariateDistribution, x::AbstractArray) +# _logpmf(d::CountableMultivariateDistribution, x::AbstractArray) """ loglikelihood(d::MultivariateDistribution, x::AbstractMatrix) diff --git a/src/samplers/discretenonparametric.jl b/src/samplers/discretenonparametric.jl index ce964a3ef8..7829cfaa3c 100644 --- a/src/samplers/discretenonparametric.jl +++ b/src/samplers/discretenonparametric.jl @@ -4,19 +4,19 @@ Data structure for efficiently sampling from an arbitrary probability mass function defined by support `xs` and probabilities `ps`. """ -struct DiscreteNonParametricSampler{T<:Real, S<:AbstractVector{T}, A<:AliasTable} <: Sampleable{Univariate,Discrete} +struct DiscreteNonParametricSampler{T, S<:AbstractVector{T}, A<:AliasTable} <: Sampleable{Univariate,Discrete} support::S aliastable::A function DiscreteNonParametricSampler{T,S}(support::S, probs::AbstractVector{<:Real} - ) where {T<:Real,S<:AbstractVector{T}} + ) where {T,S<:AbstractVector{T}} aliastable = AliasTable(probs) new{T,S,typeof(aliastable)}(support, aliastable) end end DiscreteNonParametricSampler(support::S, probs::AbstractVector{<:Real} - ) where {T<:Real,S<:AbstractVector{T}} = + ) where {T,S<:AbstractVector{T}} = DiscreteNonParametricSampler{T,S}(support, probs) rand(rng::AbstractRNG, s::DiscreteNonParametricSampler) = diff --git a/src/truncate.jl b/src/truncate.jl index 714e1fc6bf..9c40415c12 100644 --- a/src/truncate.jl +++ b/src/truncate.jl @@ -41,7 +41,7 @@ isupperbounded(d::Truncated) = isupperbounded(d.untruncated) || isfinite(d.upper minimum(d::Truncated) = max(minimum(d.untruncated), d.lower) maximum(d::Truncated) = min(maximum(d.untruncated), d.upper) -insupport(d::Truncated{D,Union{Discrete,Continuous}}, x::Real) where {D<:UnivariateDistribution} = +insupport(d::Truncated{D,VS}, x::Real) where {VS<:ValueSupport, D<:UnivariateDistribution{VS}} = d.lower <= x <= d.upper && insupport(d.untruncated, x) diff --git a/src/univariate/continuous/normal.jl b/src/univariate/continuous/normal.jl index 4de2f017a3..bb08e08fa2 100644 --- a/src/univariate/continuous/normal.jl +++ b/src/univariate/continuous/normal.jl @@ -27,7 +27,7 @@ External links * [Normal distribution on Wikipedia](http://en.wikipedia.org/wiki/Normal_distribution) """ -struct Normal{T<:Real} <: ContinuousUnivariateDistribution +struct Normal{T<:Real} <: UnivariateDistribution{ContinuousSupport{T}} μ::T σ::T Normal{T}(µ::T, σ::T) where {T<:Real} = new{T}(µ, σ) @@ -61,8 +61,6 @@ params(d::Normal) = (d.μ, d.σ) location(d::Normal) = d.μ scale(d::Normal) = d.σ -eltype(::Normal{T}) where {T} = T - #### Statistics mean(d::Normal) = d.μ @@ -169,11 +167,11 @@ invlogccdf(d::Normal, lp::Real) = xval(d, -norminvlogcdf(lp)) function quantile(d::Normal, p::Real) if iszero(d.σ) if iszero(p) - -Inf + return -Inf elseif isone(p) - Inf + return Inf else - 0.5 + return d.μ end end xval(d, -erfcinv(2p) * sqrt2) @@ -182,11 +180,11 @@ end function cquantile(d::Normal, q::Real) if iszero(d.σ) if iszero(q) - Inf + return Inf elseif isone(q) - -Inf + return -Inf else - 0.5 + return d.μ end end xval(d, erfcinv(2q) * sqrt2) diff --git a/src/univariate/discrete/bernoulli.jl b/src/univariate/discrete/bernoulli.jl index f663f04630..886a138209 100644 --- a/src/univariate/discrete/bernoulli.jl +++ b/src/univariate/discrete/bernoulli.jl @@ -80,8 +80,8 @@ end #### Evaluation -pdf(d::Bernoulli, x::Bool) = x ? succprob(d) : failprob(d) -pdf(d::Bernoulli, x::Int) = x == 0 ? failprob(d) : +pmf(d::Bernoulli, x::Bool) = x ? succprob(d) : failprob(d) +pmf(d::Bernoulli, x::Int) = x == 0 ? failprob(d) : x == 1 ? succprob(d) : zero(d.p) cdf(d::Bernoulli, x::Bool) = x ? failprob(d) : one(d.p) diff --git a/src/univariate/discrete/betabinomial.jl b/src/univariate/discrete/betabinomial.jl index 5d7e678554..5b96d212f1 100644 --- a/src/univariate/discrete/betabinomial.jl +++ b/src/univariate/discrete/betabinomial.jl @@ -83,7 +83,7 @@ function kurtosis(d::BetaBinomial) return (left * right) - 3 end -function pdf(d::BetaBinomial{T}, k::Int) where T +function pmf(d::BetaBinomial{T}, k::Int) where T n, α, β = d.n, d.α, d.β insupport(d, k) || return zero(T) chooseinv = (n + 1) * beta(k + 1, n - k + 1) @@ -92,7 +92,7 @@ function pdf(d::BetaBinomial{T}, k::Int) where T return numerator / (denominator * chooseinv) end -function logpdf(d::BetaBinomial{T}, k::Int) where T +function logpmf(d::BetaBinomial{T}, k::Int) where T n, α, β = d.n, d.α, d.β logbinom = - log1p(n) - lbeta(k + 1, n - k + 1) lognum = lbeta(k + α, n - k + β) @@ -100,12 +100,13 @@ function logpdf(d::BetaBinomial{T}, k::Int) where T logbinom + lognum - logdenom end -entropy(d::BetaBinomial) = entropy(Categorical(pdf.(Ref(d),support(d)))) -median(d::BetaBinomial) = median(Categorical(pdf.(Ref(d),support(d)))) - 1 -mode(d::BetaBinomial) = argmax(pdf.(Ref(d),support(d))) - 1 -modes(d::BetaBinomial) = modes(Categorical(pdf.(Ref(d),support(d)))) .- 1 +entropy(d::BetaBinomial) = entropy(Categorical(pmf.(Ref(d),support(d)))) +median(d::BetaBinomial) = median(Categorical(pmf.(Ref(d),support(d)))) - 1 +mode(d::BetaBinomial) = argmax(pmf.(Ref(d),support(d))) - 1 +modes(d::BetaBinomial) = modes(Categorical(pmf.(Ref(d),support(d)))) .- 1 -quantile(d::BetaBinomial, p::Float64) = quantile(Categorical(pdf.(Ref(d), support(d))), p) - 1 +quantile(d::BetaBinomial, p::Real) = + quantile(Categorical(pmf.(Ref(d), support(d))), p) - 1 #### Sampling diff --git a/src/univariate/discrete/binomial.jl b/src/univariate/discrete/binomial.jl index 7a96987d93..f17765bb5a 100644 --- a/src/univariate/discrete/binomial.jl +++ b/src/univariate/discrete/binomial.jl @@ -130,42 +130,42 @@ struct RecursiveBinomProbEvaluator{T<:Real} <: RecursiveProbabilityEvaluator end RecursiveBinomProbEvaluator(d::Binomial) = RecursiveBinomProbEvaluator(d.n, d.p / (1 - d.p)) -nextpdf(s::RecursiveBinomProbEvaluator, pv::T, x::Integer) where T <: Real = ((s.n - x + 1) / x) * s.coef * pv +nextpmf(s::RecursiveBinomProbEvaluator, pv::T, x::Integer) where T <: Real = ((s.n - x + 1) / x) * s.coef * pv -function Base.broadcast!(::typeof(pdf), r::AbstractArray, d::Binomial, X::UnitRange) +function Base.broadcast!(::typeof(pmf), r::AbstractArray, d::Binomial, X::UnitRange) axes(r) == axes(X) || throw(DimensionMismatch("Inconsistent array dimensions.")) - vl,vr, vfirst, vlast = _pdf_fill_outside!(r, d, X) + vl,vr, vfirst, vlast = _pmf_fill_outside!(r, d, X) if succprob(d) <= 1/2 # fill normal rpe = RecursiveBinomProbEvaluator(d::Binomial) - # fill central part: with non-zero pdf + # fill central part: with non-zero pmf if vl <= vr fm1 = vfirst - 1 - r[vl - fm1] = pv = pdf(d, vl) + r[vl - fm1] = pv = pmf(d, vl) for v = (vl + 1):vr - r[v - fm1] = pv = nextpdf(rpe, pv, v) + r[v - fm1] = pv = nextpmf(rpe, pv, v) end end else # fill reversed to avoid 1/0 for d.p==1. rpe = RecursiveBinomProbEvaluator(d.n, (1 - d.p) / d.p) - # fill central part: with non-zero pdf + # fill central part: with non-zero pmf if vl <= vr fm1 = vfirst - 1 - r[vr - fm1] = pv = pdf(d, vr) + r[vr - fm1] = pv = pmf(d, vr) for v = (vr-1):-1:vl - r[v - fm1] = pv = nextpdf(rpe, pv, d.n-v) + r[v - fm1] = pv = nextpmf(rpe, pv, d.n-v) end end end return r end -function Base.broadcast(::typeof(pdf), d::Binomial, X::UnitRange) +function Base.broadcast(::typeof(pmf), d::Binomial, X::UnitRange) r = similar(Array{promote_type(partype(d), eltype(X))}, axes(X)) - r .= pdf.(Ref(d),X) + r .= pmf.(Ref(d),X) end diff --git a/src/univariate/discrete/categorical.jl b/src/univariate/discrete/categorical.jl index bda56a04aa..24186afdb7 100644 --- a/src/univariate/discrete/categorical.jl +++ b/src/univariate/discrete/categorical.jl @@ -24,36 +24,41 @@ External links: * [Categorical distribution on Wikipedia](http://en.wikipedia.org/wiki/Categorical_distribution) """ -const Categorical{P,Ps} = DiscreteNonParametric{Int,P,Base.OneTo{Int},Ps} +const Categorical{P,Ps,I <: Integer} = + DiscreteNonParametric{I,P,Base.OneTo{I},Ps} -Categorical{P,Ps}(p::Ps, ::NoArgCheck) where {P<:Real, Ps<:AbstractVector{P}} = - Categorical{P,Ps}(Base.OneTo(length(p)), p, NoArgCheck()) +Categorical{P,Ps,I}(p::Ps, ::NoArgCheck) where {P<:Real, Ps<:AbstractVector{P}, I} = + Categorical{P,Ps,I}(Base.OneTo(I(length(p))), p, NoArgCheck()) Categorical(p::Ps, ::NoArgCheck) where {P<:Real, Ps<:AbstractVector{P}} = - Categorical{P,Ps}(p, NoArgCheck()) + Categorical{P,Ps,Int}(p, NoArgCheck()) -function Categorical{P,Ps}(p::Ps) where {P<:Real, Ps<:AbstractVector{P}} +function Categorical{P,Ps,I}(p::Ps) where {P<:Real, Ps<:AbstractVector{P},I} @check_args(Categorical, isprobvec(p)) - Categorical{P,Ps}(Base.OneTo(length(p)), p, NoArgCheck()) + Categorical{P,Ps,I}(Base.OneTo(I(length(p))), p, NoArgCheck()) end Categorical(p::Ps) where {P<:Real, Ps<:AbstractVector{P}} = - Categorical{P,Ps}(p) + Categorical{P,Ps,Int}(p) +Categorical{I}(p::Ps) where {P<:Real, Ps<:AbstractVector{P}, I<:Integer} = + Categorical{P,Ps,I}(p) -function Categorical(k::Integer) +function Categorical(k::I) where I <: Integer @check_args(Categorical, k >= 1) - Categorical{Float64,Vector{Float64}}(Base.OneTo(k), fill(1/k, k), NoArgCheck()) + Categorical{Float64,Vector{Float64},I}(Base.OneTo(k), fill(1/k, k), + NoArgCheck()) end ### Conversions -convert(::Type{Categorical{P,Ps}}, x::AbstractVector{<:Real}) where { - P<:Real,Ps<:AbstractVector{P}} = Categorical{P,Ps}(Ps(x)) +convert(::Type{Categorical{P,Ps,I}}, x::AbstractVector{<:Real}) where { + P<:Real,Ps<:AbstractVector{P},I} = Categorical{P,Ps,I}(Ps(x)) ### Parameters ncategories(d::Categorical) = support(d).stop -params(d::Categorical{P,Ps}) where {P<:Real, Ps<:AbstractVector{P}} = (probs(d),) +params(d::Categorical{P,Ps,I}) where {P<:Real, Ps<:AbstractVector{P}, I} = + (probs(d),) @inline partype(d::Categorical{T}) where {T<:Real} = T ### Statistics @@ -72,7 +77,7 @@ end ### Evaluation -function cdf(d::Categorical{T}, x::Int) where T<:Real +function _cdf(d::Categorical{T}, x::T) where T<:Real k = ncategories(d) p = probs(d) x < 1 && return zero(T) @@ -114,7 +119,7 @@ end # sampling -sampler(d::Categorical{P,Ps}) where {P<:Real,Ps<:AbstractVector{P}} = +sampler(d::Categorical{P,Ps,I}) where {P<:Real,Ps<:AbstractVector{P},I} = AliasTable(probs(d)) diff --git a/src/univariate/discrete/dirac.jl b/src/univariate/discrete/dirac.jl new file mode 100644 index 0000000000..fd9a8dfd7b --- /dev/null +++ b/src/univariate/discrete/dirac.jl @@ -0,0 +1,31 @@ +""" + Dirac(value) + +A *Dirac distribution* is parametrized by its only value, and takes its value with probability 1. + +```julia +d = Dirac(2.0) # Dirac distribution with value = 2. +rand(d) # Always returns the same value +``` +""" +struct Dirac{T <: Number} <: + CountableUnivariateDistribution{CountableSupport{T}} + value::T +end + +rand(::AbstractRNG, d::Dirac) = d.value +pmf(d::Dirac, x) = x ≈ d.value ? 1.0 : 0.0 +cmf(d::Dirac, x) = x < d.value ? 0.0 : 1.0 +quantile(d::Dirac{T}, p::Real) where T = 0 ≤ p ≤ 1 ? d.value : T(NaN) +minimum(d::Dirac) = d.value +maximum(d::Dirac) = d.value +support(d::Dirac) = [d.value] +insupport(d::Dirac, x) = x == d.value +mean(d::Dirac) = d.value +var(d::Dirac) = 0.0 +mode(d::Dirac) = d.value +skewness(d::Dirac) = 0.0 +kurtosis(d::Dirac) = 0.0 +entropy(d::Dirac) = 0.0 +mgf(d::Dirac, t) = exp(t * d.value) +cf(d::Dirac, t) = exp(im * t * d.value) diff --git a/src/univariate/discrete/discretenonparametric.jl b/src/univariate/discrete/discretenonparametric.jl index c36b5e76b3..284dee4b7d 100644 --- a/src/univariate/discrete/discretenonparametric.jl +++ b/src/univariate/discrete/discretenonparametric.jl @@ -2,7 +2,7 @@ DiscreteNonParametric(xs, ps) A *Discrete nonparametric distribution* explicitly defines an arbitrary -probability mass function in terms of a list of real support values and their +probability mass function in terms of a list of support values and their corresponding probabilities ```julia @@ -17,16 +17,18 @@ External links * [Probability mass function on Wikipedia](http://en.wikipedia.org/wiki/Probability_mass_function) """ -struct DiscreteNonParametric{T<:Real,P<:Real,Ts<:AbstractVector{T},Ps<:AbstractVector{P}} <: DiscreteUnivariateDistribution +struct DiscreteNonParametric{T,P<:Real, + Ts<:AbstractVector{T},Ps<:AbstractVector{P}} <: + UnivariateDistribution{CountableSupport{T}} support::Ts p::Ps DiscreteNonParametric{T,P,Ts,Ps}(vs::Ts, ps::Ps, ::NoArgCheck) where { - T<:Real,P<:Real,Ts<:AbstractVector{T},Ps<:AbstractVector{P}} = + T,P<:Real,Ts<:AbstractVector{T},Ps<:AbstractVector{P}} = new{T,P,Ts,Ps}(vs, ps) function DiscreteNonParametric{T,P,Ts,Ps}(vs::Ts, ps::Ps) where { - T<:Real,P<:Real,Ts<:AbstractVector{T},Ps<:AbstractVector{P}} + T,P<:Real,Ts<:AbstractVector{T},Ps<:AbstractVector{P}} @check_args(DiscreteNonParametric, length(vs) == length(ps)) @check_args(DiscreteNonParametric, isprobvec(ps)) @check_args(DiscreteNonParametric, allunique(vs)) @@ -36,14 +38,20 @@ struct DiscreteNonParametric{T<:Real,P<:Real,Ts<:AbstractVector{T},Ps<:AbstractV end DiscreteNonParametric(vs::Ts, ps::Ps) where { - T<:Real,P<:Real,Ts<:AbstractVector{T},Ps<:AbstractVector{P}} = + T,P<:Real,Ts<:AbstractVector{T},Ps<:AbstractVector{P}} = DiscreteNonParametric{T,P,Ts,Ps}(vs, ps) DiscreteNonParametric(vs::Ts, ps::Ps, a::NoArgCheck) where { - T<:Real,P<:Real,Ts<:AbstractVector{T},Ps<:AbstractVector{P}} = + T,P<:Real,Ts<:AbstractVector{T},Ps<:AbstractVector{P}} = DiscreteNonParametric{T,P,Ts,Ps}(vs, ps, a) -eltype(d::DiscreteNonParametric{T}) where T = T +DiscreteNonParametric(dict::Dict{K, V}) where {K,V<:Real} = + DiscreteNonParametric{K,V,Vector{K},Vector{V}}(collect(keys(dict)), + collect(values(dict))) + +DiscreteNonParametric(dict::Dict{K, V}, a::NoArgCheck) where {K,V<:Real} = + DiscreteNonParametric{K,V,Vector{K},Vector{V}}(collect(keys(dict)), + collect(values(dict)), a) # Conversion convert(::Type{DiscreteNonParametric{T,P,Ts,Ps}}, d::DiscreteNonParametric) where {T,P,Ts,Ps} = @@ -88,8 +96,6 @@ function rand(rng::AbstractRNG, d::DiscreteNonParametric{T,P}) where {T,P} x[max(i,1)] end -rand(d::DiscreteNonParametric) = rand(GLOBAL_RNG, d) - sampler(d::DiscreteNonParametric) = DiscreteNonParametricSampler(support(d), probs(d)) @@ -99,11 +105,11 @@ get_evalsamples(d::DiscreteNonParametric, ::Float64) = support(d) # Evaluation -pdf(d::DiscreteNonParametric) = copy(probs(d)) +pmf(d::DiscreteNonParametric) = copy(probs(d)) # Helper functions for pdf and cdf required to fix ambiguous method -# error involving [pc]df(::DisceteUnivariateDistribution, ::Int) -function _pdf(d::DiscreteNonParametric{T,P}, x::T) where {T,P} +# error involving [pc]mf(::DisceteUnivariateDistribution, ::Int) +function _pmf(d::DiscreteNonParametric{T,P}, x::T) where {T,P} idx_range = searchsorted(support(d), x) if length(idx_range) > 0 return probs(d)[first(idx_range)] @@ -111,8 +117,9 @@ function _pdf(d::DiscreteNonParametric{T,P}, x::T) where {T,P} return zero(P) end end -pdf(d::DiscreteNonParametric{T}, x::Int) where T = _pdf(d, convert(T, x)) -pdf(d::DiscreteNonParametric{T}, x::Real) where T = _pdf(d, convert(T, x)) +pmf(d::DiscreteNonParametric{T}, x::Int) where T = _pmf(d, convert(T, x)) +pmf(d::DiscreteNonParametric{T}, x::Real) where T = _pmf(d, convert(T, x)) +pmf(d::DiscreteNonParametric{T}, x) where T = _pmf(d, convert(T, x)) function _cdf(d::DiscreteNonParametric{T,P}, x::T) where {T,P} x > maximum(d) && return 1.0 @@ -126,6 +133,7 @@ function _cdf(d::DiscreteNonParametric{T,P}, x::T) where {T,P} end cdf(d::DiscreteNonParametric{T}, x::Integer) where T = _cdf(d, convert(T, x)) cdf(d::DiscreteNonParametric{T}, x::Real) where T = _cdf(d, convert(T, x)) +cdf(d::DiscreteNonParametric{T}, x) where T = _cdf(d, convert(T, x)) function _ccdf(d::DiscreteNonParametric{T,P}, x::T) where {T,P} x < minimum(d) && return 1.0 @@ -139,6 +147,7 @@ function _ccdf(d::DiscreteNonParametric{T,P}, x::T) where {T,P} end ccdf(d::DiscreteNonParametric{T}, x::Integer) where T = _ccdf(d, convert(T, x)) ccdf(d::DiscreteNonParametric{T}, x::Real) where T = _ccdf(d, convert(T, x)) +ccdf(d::DiscreteNonParametric{T}, x) where T = _ccdf(d, convert(T, x)) function quantile(d::DiscreteNonParametric, q::Real) 0 <= q <= 1 || throw(DomainError()) @@ -156,8 +165,10 @@ end minimum(d::DiscreteNonParametric) = first(support(d)) maximum(d::DiscreteNonParametric) = last(support(d)) -insupport(d::DiscreteNonParametric, x::Real) = +insupport(d::DiscreteNonParametric{T}, x::T) where {T} = length(searchsorted(support(d), x)) > 0 +insupport(d::DiscreteNonParametric{T}, x::Number) where {T<:Number} = + length(searchsorted(support(d), convert(T, x))) > 0 mean(d::DiscreteNonParametric) = dot(probs(d), support(d)) diff --git a/src/univariate/discrete/discreteuniform.jl b/src/univariate/discrete/discreteuniform.jl index a8c0f5a647..163ac1ed11 100644 --- a/src/univariate/discrete/discreteuniform.jl +++ b/src/univariate/discrete/discreteuniform.jl @@ -74,11 +74,11 @@ cdf(d::DiscreteUniform, x::Int) = (x < d.a ? 0.0 : x > d.b ? 1.0 : (floor(Int,x) - d.a + 1.0) * d.pv) -pdf(d::DiscreteUniform, x::Int) = insupport(d, x) ? d.pv : 0.0 +pmf(d::DiscreteUniform, x::Int) = insupport(d, x) ? d.pv : 0.0 -logpdf(d::DiscreteUniform, x::Int) = insupport(d, x) ? log(d.pv) : -Inf +logpmf(d::DiscreteUniform, x::Int) = insupport(d, x) ? log(d.pv) : -Inf -quantile(d::DiscreteUniform, p::Float64) = d.a + floor(Int,p * span(d)) +quantile(d::DiscreteUniform, p::Real) = d.a + floor(Int,p * span(d)) function mgf(d::DiscreteUniform, t::T) where {T <: Real} a, b = d.a, d.b diff --git a/src/univariate/discrete/geometric.jl b/src/univariate/discrete/geometric.jl index e26c14d7d8..547a22cb95 100644 --- a/src/univariate/discrete/geometric.jl +++ b/src/univariate/discrete/geometric.jl @@ -70,7 +70,7 @@ entropy(d::Geometric) = (-xlogx(succprob(d)) - xlogx(failprob(d))) / d.p ### Evaluations -function pdf(d::Geometric{T}, x::Int) where T<:Real +function pmf(d::Geometric{T}, x::Int) where T<:Real if x >= 0 p = d.p return p < one(p) / 10 ? p * exp(log1p(-p) * x) : d.p * (one(p) - p)^x @@ -79,7 +79,7 @@ function pdf(d::Geometric{T}, x::Int) where T<:Real end end -function logpdf(d::Geometric{T}, x::Int) where T<:Real +function logpmf(d::Geometric{T}, x::Int) where T<:Real x >= 0 ? log(d.p) + log1p(-d.p) * x : -T(Inf) end @@ -88,13 +88,13 @@ struct RecursiveGeomProbEvaluator <: RecursiveProbabilityEvaluator end RecursiveGeomProbEvaluator(d::Geometric) = RecursiveGeomProbEvaluator(failprob(d)) -nextpdf(s::RecursiveGeomProbEvaluator, p::Real, x::Integer) = p * s.p0 +nextpmf(s::RecursiveGeomProbEvaluator, p::Real, x::Integer) = p * s.p0 -Base.broadcast!(::typeof(pdf), r::AbstractArray, d::Geometric, rgn::UnitRange) = - _pdf!(r, d, rgn, RecursiveGeomProbEvaluator(d)) -function Base.broadcast(::typeof(pdf), d::Geometric, X::UnitRange) +Base.broadcast!(::typeof(pmf), r::AbstractArray, d::Geometric, rgn::UnitRange) = + _pmf!(r, d, rgn, RecursiveGeomProbEvaluator(d)) +function Base.broadcast(::typeof(pmf), d::Geometric, X::UnitRange) r = similar(Array{promote_type(partype(d), eltype(X))}, axes(X)) - r .= pdf.(Ref(d),X) + r .= pmf.(Ref(d),X) end diff --git a/src/univariate/discrete/hypergeometric.jl b/src/univariate/discrete/hypergeometric.jl index 72e6f2a14e..4b093bee51 100644 --- a/src/univariate/discrete/hypergeometric.jl +++ b/src/univariate/discrete/hypergeometric.jl @@ -74,7 +74,7 @@ function kurtosis(d::Hypergeometric) a/b end -entropy(d::Hypergeometric) = entropy(pdf.(Ref(d), support(d))) +entropy(d::Hypergeometric) = entropy(pmf.(Ref(d), support(d))) ### Evaluation & Sampling @@ -99,13 +99,13 @@ end RecursiveHypergeomProbEvaluator(d::Hypergeometric) = RecursiveHypergeomProbEvaluator(d.ns, d.nf, d.n) -nextpdf(s::RecursiveHypergeomProbEvaluator, p::Float64, x::Integer) = +nextpmf(s::RecursiveHypergeomProbEvaluator, p::Float64, x::Integer) = ((s.ns - x + 1) / x) * ((s.n - x + 1) / (s.nf - s.n + x)) * p -Base.broadcast!(::typeof(pdf), r::AbstractArray, d::Hypergeometric, rgn::UnitRange) = - _pdf!(r, d, rgn, RecursiveHypergeomProbEvaluator(d)) +Base.broadcast!(::typeof(pmf), r::AbstractArray, d::Hypergeometric, rgn::UnitRange) = + _pmf!(r, d, rgn, RecursiveHypergeomProbEvaluator(d)) -function Base.broadcast(::typeof(pdf), d::Hypergeometric, X::UnitRange) +function Base.broadcast(::typeof(pmf), d::Hypergeometric, X::UnitRange) r = similar(Array{promote_type(partype(d), eltype(X))}, axes(X)) - r .= pdf.(Ref(d),X) + r .= pmf.(Ref(d),X) end diff --git a/src/univariate/discrete/negativebinomial.jl b/src/univariate/discrete/negativebinomial.jl index 682ebe171a..512163ebda 100644 --- a/src/univariate/discrete/negativebinomial.jl +++ b/src/univariate/discrete/negativebinomial.jl @@ -105,13 +105,13 @@ struct RecursiveNegBinomProbEvaluator <: RecursiveProbabilityEvaluator end RecursiveNegBinomProbEvaluator(d::NegativeBinomial) = RecursiveNegBinomProbEvaluator(d.r, failprob(d)) -nextpdf(s::RecursiveNegBinomProbEvaluator, p::Float64, x::Integer) = ((x + s.r - 1) / x) * s.p0 * p +nextpmf(s::RecursiveNegBinomProbEvaluator, p::Float64, x::Integer) = ((x + s.r - 1) / x) * s.p0 * p -Base.broadcast!(::typeof(pdf), r::AbstractArray, d::NegativeBinomial, rgn::UnitRange) = - _pdf!(r, d, rgn, RecursiveNegBinomProbEvaluator(d)) -function Base.broadcast(::typeof(pdf), d::NegativeBinomial, X::UnitRange) +Base.broadcast!(::typeof(pmf), r::AbstractArray, d::NegativeBinomial, rgn::UnitRange) = + _pmf!(r, d, rgn, RecursiveNegBinomProbEvaluator(d)) +function Base.broadcast(::typeof(pmf), d::NegativeBinomial, X::UnitRange) r = similar(Array{promote_type(partype(d), eltype(X))}, axes(X)) - r .= pdf.(Ref(d),X) + r .= pmf.(Ref(d),X) end function mgf(d::NegativeBinomial, t::Real) diff --git a/src/univariate/discrete/noncentralhypergeometric.jl b/src/univariate/discrete/noncentralhypergeometric.jl index ac6d8b0fde..2639cf99e3 100644 --- a/src/univariate/discrete/noncentralhypergeometric.jl +++ b/src/univariate/discrete/noncentralhypergeometric.jl @@ -26,7 +26,7 @@ function quantile(d::NoncentralHypergeometric{T}, q::Real) where T<:Real qsum, i = zero(T), 0 while qsum < q i += 1 - qsum += pdf(d, range[i]) + qsum += pmf(d, range[i]) end range[i] end @@ -81,12 +81,12 @@ mode(d::FisherNoncentralHypergeometric) = floor(Int, _mode(d)) testfd(d::FisherNoncentralHypergeometric) = d.ω^3 -logpdf(d::FisherNoncentralHypergeometric, k::Int) = +logpmf(d::FisherNoncentralHypergeometric, k::Int) = -log(d.ns + 1) - lbeta(d.ns - k + 1, k + 1) - log(d.nf + 1) - lbeta(d.nf - d.n + k + 1, d.n - k + 1) + xlogy(k, d.ω) - _P(d, 0) -pdf(d::FisherNoncentralHypergeometric, k::Int) = exp(logpdf(d, k)) +pmf(d::FisherNoncentralHypergeometric, k::Int) = exp(logpmf(d, k)) ## Wallenius' noncentral hypergeometric distribution @@ -114,15 +114,15 @@ convert(::Type{WalleniusNoncentralHypergeometric{T}}, ns::Real, nf::Real, n::Rea convert(::Type{WalleniusNoncentralHypergeometric{T}}, d::WalleniusNoncentralHypergeometric{S}) where {T<:Real, S<:Real} = WalleniusNoncentralHypergeometric(d.ns, d.nf, d.n, T(d.ω)) # Properties -mean(d::WalleniusNoncentralHypergeometric) = sum(support(d) .* pdf.(Ref(d), support(d))) -var(d::WalleniusNoncentralHypergeometric) = sum((support(d) .- mean(d)).^2 .* pdf.(Ref(d), support(d))) -mode(d::WalleniusNoncentralHypergeometric) = support(d)[argmax(pdf.(Ref(d), support(d)))] +mean(d::WalleniusNoncentralHypergeometric) = sum(support(d) .* pmf.(Ref(d), support(d))) +var(d::WalleniusNoncentralHypergeometric) = sum((support(d) .- mean(d)).^2 .* pmf.(Ref(d), support(d))) +mode(d::WalleniusNoncentralHypergeometric) = support(d)[argmax(pmf.(Ref(d), support(d)))] entropy(d::WalleniusNoncentralHypergeometric) = 1 testfd(d::WalleniusNoncentralHypergeometric) = d.ω^3 -function logpdf(d::WalleniusNoncentralHypergeometric, k::Int) +function logpmf(d::WalleniusNoncentralHypergeometric, k::Int) D = d.ω * (d.ns - k) + (d.nf - d.n + k) f(t) = (1 - t^(d.ω / D))^k * (1 - t^(1 / D))^(d.n - k) I, _ = quadgk(f, 0, 1) @@ -130,4 +130,4 @@ function logpdf(d::WalleniusNoncentralHypergeometric, k::Int) log(d.nf + 1) - lbeta(d.nf - d.n + k + 1, d.n - k + 1) + log(I) end -pdf(d::WalleniusNoncentralHypergeometric, k::Int) = exp(logpdf(d, k)) +pmf(d::WalleniusNoncentralHypergeometric, k::Int) = exp(logpmf(d, k)) diff --git a/src/univariate/discrete/poisson.jl b/src/univariate/discrete/poisson.jl index 8f3ec3bb3f..60adb36eb2 100644 --- a/src/univariate/discrete/poisson.jl +++ b/src/univariate/discrete/poisson.jl @@ -95,13 +95,13 @@ struct RecursivePoissonProbEvaluator <: RecursiveProbabilityEvaluator end RecursivePoissonProbEvaluator(d::Poisson) = RecursivePoissonProbEvaluator(rate(d)) -nextpdf(s::RecursivePoissonProbEvaluator, p::Float64, x::Integer) = p * s.λ / x +nextpmf(s::RecursivePoissonProbEvaluator, p::Float64, x::Integer) = p * s.λ / x -Base.broadcast!(::typeof(pdf), r::AbstractArray, d::Poisson, rgn::UnitRange) = - _pdf!(r, d, rgn, RecursivePoissonProbEvaluator(d)) -function Base.broadcast(::typeof(pdf), d::Poisson, X::UnitRange) +Base.broadcast!(::typeof(pmf), r::AbstractArray, d::Poisson, rgn::UnitRange) = + _pmf!(r, d, rgn, RecursivePoissonProbEvaluator(d)) +function Base.broadcast(::typeof(pmf), d::Poisson, X::UnitRange) r = similar(Array{promote_type(partype(d), eltype(X))}, axes(X)) - r .= pdf.(Ref(d),X) + r .= pmf.(Ref(d),X) end diff --git a/src/univariate/discrete/poissonbinomial.jl b/src/univariate/discrete/poissonbinomial.jl index 37a7bcfff3..20233c444e 100644 --- a/src/univariate/discrete/poissonbinomial.jl +++ b/src/univariate/discrete/poissonbinomial.jl @@ -27,9 +27,8 @@ External links: struct PoissonBinomial{T<:Real} <: DiscreteUnivariateDistribution p::Vector{T} pmf::Vector{T} - function PoissonBinomial{T}(p::AbstractArray) where {T <: Real} - pb = poissonbinomial_pdf_fft(p) + pb = poissonbinomial_pmf_fft(p) @assert isprobvec(pb) new{T}(p, pb) end @@ -98,7 +97,7 @@ modes(d::PoissonBinomial) = [x - 1 for x in modes(Categorical(d.pmf))] #### Evaluation -quantile(d::PoissonBinomial, x::Float64) = quantile(Categorical(d.pmf), x) - 1 +quantile(d::PoissonBinomial, x::Real) = quantile(Categorical(d.pmf), x) - 1 function mgf(d::PoissonBinomial{T}, t::Real) where {T} p, = params(d) @@ -110,20 +109,20 @@ function cf(d::PoissonBinomial{T}, t::Real) where {T} prod(one(T) .- p .+ p .* cis(t)) end -pdf(d::PoissonBinomial, k::Integer) = insupport(d, k) ? d.pmf[k+1] : 0 -function logpdf(d::PoissonBinomial{T}, k::Int) where T<:Real +pmf(d::PoissonBinomial, k::Integer) = insupport(d, k) ? d.pmf[k+1] : 0 +function logpmf(d::PoissonBinomial{T}, k::Int) where T<:Real insupport(d, k) ? log(d.pmf[k + 1]) : -T(Inf) end -# Computes the pdf of a poisson-binomial random variable using +# Computes the pmf of a poisson-binomial random variable using # fast fourier transform # # Hong, Y. (2013). # On computing the distribution function for the Poisson binomial # distribution. Computational Statistics and Data Analysis, 59, 41–51. # -function poissonbinomial_pdf_fft(p::AbstractArray) +function poissonbinomial_pmf_fft(p::AbstractArray) n = length(p) ω = 2 / (n + 1) diff --git a/src/univariate/discrete/skellam.jl b/src/univariate/discrete/skellam.jl index beb538e832..ec6b946666 100644 --- a/src/univariate/discrete/skellam.jl +++ b/src/univariate/discrete/skellam.jl @@ -69,12 +69,12 @@ kurtosis(d::Skellam) = 1 / var(d) #### Evaluation -function logpdf(d::Skellam, x::Integer) +function logpmf(d::Skellam, x::Integer) μ1, μ2 = params(d) - (μ1 + μ2) + (x/2) * log(μ1/μ2) + log(besseli(x, 2*sqrt(μ1)*sqrt(μ2))) end -pdf(d::Skellam, x::Integer) = exp(logpdf(d, x)) +pmf(d::Skellam, x::Integer) = exp(logpmf(d, x)) function mgf(d::Skellam, t::Real) μ1, μ2 = params(d) @@ -92,8 +92,8 @@ end Implementation based on SciPy: https://github.com/scipy/scipy/blob/v0.15.1/scipy/stats/_discrete_distns.py Refer to Eqn (5) in On an Extension of the Connexion Between Poisson and χ2 Distributions, N.L Johnson(1959) -Vol 46, No 3/4, doi:10.2307/2333532 -It relates the Skellam and Non-central chisquare PDFs, which is very similar to their CDFs computation as well. +Vol 46, No 3/4, doi:10.2307/2333532 +It relates the Skellam and Non-central chisquare PMFs, which is very similar to their CDFs computation as well. Computing cdf of the Skellam distribution. """ diff --git a/src/univariates.jl b/src/univariates.jl index fe97d1f6c9..59ac9f1360 100644 --- a/src/univariates.jl +++ b/src/univariates.jl @@ -17,7 +17,7 @@ isbounded(d::Union{D,Type{D}}) where {D<:UnivariateDistribution} = isupperbounde islowerbounded(d::Union{D,Type{D}}) where {D<:UnivariateDistribution} = minimum(d) > -Inf isupperbounded(d::Union{D,Type{D}}) where {D<:UnivariateDistribution} = maximum(d) < +Inf -hasfinitesupport(d::Union{D,Type{D}}) where {D<:DiscreteUnivariateDistribution} = isbounded(d) +hasfinitesupport(d::Union{D,Type{D}}) where {D<:CountableUnivariateDistribution} = isbounded(d) hasfinitesupport(d::Union{D,Type{D}}) where {D<:ContinuousUnivariateDistribution} = false """ @@ -124,10 +124,11 @@ end insupport(d::Union{D,Type{D}}, X::AbstractArray) where {D<:UnivariateDistribution} = insupport!(BitArray(undef, size(X)), d, X) -insupport(d::Union{D,Type{D}},x::Real) where {D<:ContinuousUnivariateDistribution} = minimum(d) <= x <= maximum(d) +insupport(d::Union{D,Type{D}},x::Real) where {D<:UnivariateDistribution{ContinuousSupport{T}}} where {T} = minimum(d) <= x <= maximum(d) +insupport(d::D,x::T) where {T, D<:UnivariateDistribution{CountableSupport{T}}} = x ∈ support(d) insupport(d::Union{D,Type{D}},x::Real) where {D<:DiscreteUnivariateDistribution} = isinteger(x) && minimum(d) <= x <= maximum(d) -support(d::Union{D,Type{D}}) where {D<:ContinuousUnivariateDistribution} = RealInterval(minimum(d), maximum(d)) +support(d::Union{D,Type{D}}) where {D<:UnivariateDistribution{ContinuousSupport{T}}} where {T} = RealInterval(minimum(d), maximum(d)) support(d::Union{D,Type{D}}) where {D<:DiscreteUnivariateDistribution} = round(Int, minimum(d)):round(Int, maximum(d)) # Type used for dispatch on finite support @@ -308,7 +309,7 @@ cf(d::UnivariateDistribution, t) #### pdf, cdf, and friends -# pdf +# pmf """ pdf(d::UnivariateDistribution, x::Real) @@ -320,7 +321,7 @@ See also: [`logpdf`](@ref). pdf(d::UnivariateDistribution, x::Real) """ - pdf(d::DiscreteUnivariateDistribution, x::T) where {T<:Real} + pmf(d::CountableUnivariateDistribution, x::T) where {T<:Real} Evaluate the probability density (mass) at `x`. If `T` is not an `Integer` type but `x` is integer, the value is converted to `Int`. @@ -328,9 +329,11 @@ type but `x` is integer, the value is converted to `Int`. The version with `x::Integer` must be implemented by discrete distributions. -See also: [`logpdf`](@ref). +See also: [`logpmf`](@ref). """ -pdf(d::DiscreteUnivariateDistribution, x::Real) = isinteger(x) ? pdf(d, round(Int, x)) : 0.0 +pmf(d::UnivariateDistribution{C}, x::Real) where {I <: Integer, + C <: CountableSupport{I}} = + insupport(d, x) ? pmf(d, round(Int, x)) : 0.0 """ logpdf(d::UnivariateDistribution, x::Real) @@ -339,9 +342,10 @@ Evaluate the logarithm of probability density (mass) at `x`. Whereas there is a fallback implemented `logpdf(d, x) = log(pdf(d, x))`. Relying on this fallback is not recommended in general, as it is prone to overflow or underflow. """ -logpdf(d::UnivariateDistribution, x::Real) = log(pdf(d, x)) -logpdf(d::DiscreteUnivariateDistribution, x::Integer) = log(pdf(d, x)) -logpdf(d::DiscreteUnivariateDistribution, x::Real) = isinteger(x) ? logpdf(d, round(Int, x)) : -Inf +logpdf(d::ContinuousUnivariateDistribution, x::Real) = log(pdf(d, x)) +logpmf(d::DiscreteUnivariateDistribution, x::Integer) = log(pmf(d, x)) +logpmf(d::CountableUnivariateDistribution, x) = log(pmf(d, x)) +logpmf(d::DiscreteUnivariateDistribution, x::Real) = isinteger(x) ? logpmf(d, round(Int, x)) : -Inf """ cdf(d::UnivariateDistribution, x::Real) @@ -350,20 +354,27 @@ Evaluate the cumulative probability at `x`. See also [`ccdf`](@ref), [`logcdf`](@ref), and [`logccdf`](@ref). """ -cdf(d::UnivariateDistribution, x::Real) -cdf(d::DiscreteUnivariateDistribution, x::Integer) = cdf(d, x, FiniteSupport{hasfinitesupport(d)}) +cdf(d::UnivariateDistribution, x) + +cdf(d::CountableUnivariateDistribution, x) = + cdf(d, x, FiniteSupport{hasfinitesupport(d)}) # Discrete univariate with infinite support -function cdf(d::DiscreteUnivariateDistribution, x::Integer, ::Type{FiniteSupport{false}}) +function cdf(d::DiscreteUnivariateDistribution, x::Integer, + ::Type{FiniteSupport{false}}) c = 0.0 for y = minimum(d):x - c += pdf(d, y) + c += pmf(d, y) end return c end +cdf(d::DiscreteUnivariateDistribution, x::Real, ::Type{FiniteSupport{false}}) = + cdf(d, floor(Int,x), FiniteSupport{false}) + # Discrete univariate with finite support -function cdf(d::DiscreteUnivariateDistribution, x::Integer, ::Type{FiniteSupport{true}}) +function cdf(d::DiscreteUnivariateDistribution, x::Integer, + ::Type{FiniteSupport{true}}) # calculate from left if x < (min + max)/2 # (same as infinite support version) x <= div(minimum(d) + maximum(d),2) && return cdf(d, x, FiniteSupport{false}) @@ -376,7 +387,18 @@ function cdf(d::DiscreteUnivariateDistribution, x::Integer, ::Type{FiniteSupport return c end -cdf(d::DiscreteUnivariateDistribution, x::Real) = cdf(d, floor(Int,x)) +function cdf(d::CountableUnivariateDistribution, x, ::Type{FiniteSupport{true}}) + # calculate from left if in first half of support + s = support(d) + val = s[floor(Int, length(s) / 2 + 1)] + if x <= val + return reduce(+, pmf(d, y) for y in s if y <= x; init = 0.0) + else + # otherwise, calculate from the right + return 1.0 - reduce(+, pmf(d, y) for y in s if y > x; init = 0.0) + end +end + cdf(d::ContinuousUnivariateDistribution, x::Real) = throw(MethodError(cdf, (d, x))) @@ -385,7 +407,7 @@ cdf(d::ContinuousUnivariateDistribution, x::Real) = throw(MethodError(cdf, (d, x The complementary cumulative function evaluated at `x`, i.e. `1 - cdf(d, x)`. """ -ccdf(d::UnivariateDistribution, x::Real) = 1.0 - cdf(d, x) +ccdf(d::UnivariateDistribution, x) = 1.0 - cdf(d, x) ccdf(d::DiscreteUnivariateDistribution, x::Integer) = 1.0 - cdf(d, x) ccdf(d::DiscreteUnivariateDistribution, x::Real) = ccdf(d, floor(Int,x)) @@ -394,18 +416,20 @@ ccdf(d::DiscreteUnivariateDistribution, x::Real) = ccdf(d, floor(Int,x)) The logarithm of the cumulative function value(s) evaluated at `x`, i.e. `log(cdf(x))`. """ -logcdf(d::UnivariateDistribution, x::Real) = log(cdf(d, x)) +logcdf(d::UnivariateDistribution, x) = log(cdf(d, x)) logcdf(d::DiscreteUnivariateDistribution, x::Integer) = log(cdf(d, x)) logcdf(d::DiscreteUnivariateDistribution, x::Real) = logcdf(d, floor(Int,x)) +logcdf(d::CountableUnivariateDistribution, x) = log(cdf(d, x)) """ logccdf(d::UnivariateDistribution, x::Real) The logarithm of the complementary cumulative function values evaluated at x, i.e. `log(ccdf(x))`. """ -logccdf(d::UnivariateDistribution, x::Real) = log(ccdf(d, x)) +logccdf(d::UnivariateDistribution, x) = log(ccdf(d, x)) logccdf(d::DiscreteUnivariateDistribution, x::Integer) = log(ccdf(d, x)) logccdf(d::DiscreteUnivariateDistribution, x::Real) = logccdf(d, floor(Int,x)) +logccdf(d::CountableUnivariateDistribution, x) = log(ccdf(d, x)) """ quantile(d::UnivariateDistribution, q::Real) @@ -421,7 +445,7 @@ 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) +cquantile(d::UnivariateDistribution, p::Real) = quantile(d, one(p) - p) """ invlogcdf(d::UnivariateDistribution, lp::Real) @@ -495,15 +519,15 @@ end abstract type RecursiveProbabilityEvaluator end -function _pdf!(r::AbstractArray, d::DiscreteUnivariateDistribution, X::UnitRange, rpe::RecursiveProbabilityEvaluator) - vl,vr, vfirst, vlast = _pdf_fill_outside!(r, d, X) +function _pmf!(r::AbstractArray, d::DiscreteUnivariateDistribution, X::UnitRange, rpe::RecursiveProbabilityEvaluator) + vl,vr, vfirst, vlast = _pmf_fill_outside!(r, d, X) - # fill central part: with non-zero pdf + # fill central part: with non-zero pmf if vl <= vr fm1 = vfirst - 1 - r[vl - fm1] = pv = pdf(d, vl) + r[vl - fm1] = pv = pmf(d, vl) for v = (vl+1):vr - r[v - fm1] = pv = nextpdf(rpe, pv, v) + r[v - fm1] = pv = nextpmf(rpe, pv, v) end end @@ -527,6 +551,8 @@ macro _delegate_statsfuns(D, fpre, psyms...) # function names from StatsFuns fpdf = Symbol(fpre, "pdf") flogpdf = Symbol(fpre, "logpdf") + fpmf = Symbol(fpre, "pdf") + flogpmf = Symbol(fpre, "logpdf") fcdf = Symbol(fpre, "cdf") fccdf = Symbol(fpre, "ccdf") flogcdf = Symbol(fpre, "logcdf") @@ -542,6 +568,8 @@ macro _delegate_statsfuns(D, fpre, psyms...) esc(quote pdf(d::$D, x::$T) = $(fpdf)($(pargs...), x) logpdf(d::$D, x::$T) = $(flogpdf)($(pargs...), x) + pmf(d::$D, x::$T) = $(fpmf)($(pargs...), x) + logpmf(d::$D, x::$T) = $(flogpmf)($(pargs...), x) cdf(d::$D, x::$T) = $(fcdf)($(pargs...), x) ccdf(d::$D, x::$T) = $(fccdf)($(pargs...), x) @@ -562,6 +590,7 @@ const discrete_distributions = [ "bernoulli", "betabinomial", "binomial", + "dirac", "discreteuniform", "discretenonparametric", "categorical", diff --git a/test/categorical.jl b/test/categorical.jl index e2691c9a2e..0b1ca810e4 100644 --- a/test/categorical.jl +++ b/test/categorical.jl @@ -52,5 +52,5 @@ println(" testing $d as Categorical") p = ones(10^6) * 1.0e-6 @test Distributions.isprobvec(p) -@test typeof(convert(Categorical{Float32,Vector{Float32}}, d)) == Categorical{Float32,Vector{Float32}} -@test typeof(convert(Categorical{Float32,Vector{Float32}}, d.p)) == Categorical{Float32,Vector{Float32}} +@test typeof(convert(Categorical{Float32,Vector{Float32},Int}, d)) == Categorical{Float32,Vector{Float32},Int} +@test typeof(convert(Categorical{Float32,Vector{Float32},Int32}, d.p)) == Categorical{Float32,Vector{Float32},Int32} diff --git a/test/dirac.jl b/test/dirac.jl new file mode 100644 index 0000000000..73a6fbced9 --- /dev/null +++ b/test/dirac.jl @@ -0,0 +1,20 @@ +using Distributions +using Test + +@testset "Dirac tests" begin + for val in [3, 3.0] + d = Dirac(val) + @test minimum(d) == val + @test maximum(d) == val + @test pdf(d, val - 1) == 0 + @test pdf(d, val) == 1 + @test pdf(d, val + 1) == 0 + @test cdf(d, val - 1) == 0 + @test cdf(d, val) == 1 + @test cdf(d, val + 1) == 1 + @test quantile(d, 0) == val + @test quantile(d, 0.5) == val + @test quantile(d, 1) == val + @test rand(d) == val + end +end diff --git a/test/discretenonparametric.jl b/test/discretenonparametric.jl index 7a0e76bd40..fcf121929b 100644 --- a/test/discretenonparametric.jl +++ b/test/discretenonparametric.jl @@ -9,31 +9,47 @@ rng = MersenneTwister(123) "rand(rng, ...)" => [dist -> rand(rng, dist), (dist, n) -> rand(rng, dist, n)]) d = DiscreteNonParametric([40., 80., 120., -60.], [.4, .3, .1, .2]) +d_string = DiscreteNonParametric(["forty", "eighty", "one hundred and twenty", "negative sixty"], [.4, .3, .1, .2]) @test !(d ≈ DiscreteNonParametric([40., 80, 120, -60], [.4, .3, .1, .2], Distributions.NoArgCheck())) @test d ≈ DiscreteNonParametric([-60., 40., 80, 120], [.2, .4, .3, .1], Distributions.NoArgCheck()) # Invalid probability @test_throws ArgumentError DiscreteNonParametric([40., 80, 120, -60], [.5, .3, .1, .2]) +@test_throws ArgumentError DiscreteNonParametric(["forty", "eighty", "one hundred and twenty", "negative sixty"], [.5, .3, .1, .2]) # Invalid probability, but no arg check DiscreteNonParametric([40., 80, 120, -60], [.5, .3, .1, .2], Distributions.NoArgCheck()) +DiscreteNonParametric(["forty", "eighty", "one hundred and twenty", "negative sixty"], [.5, .3, .1, .2], Distributions.NoArgCheck()) test_range(d) vs = Distributions.get_evalsamples(d, 0.00001) +vs_string = Distributions.get_evalsamples(d_string, 0.00001) test_evaluation(d, vs, true) +test_evaluation(d_string, vs_string, true) test_stats(d, vs) test_params(d) +test_params(d_string) @test func[1](d) ∈ [40., 80., 120., -60.] +@test func[1](d_string) ∈ ["forty", "eighty", "one hundred and twenty", "negative sixty"] +@test func[1](d) ∈ [40., 80., 120., -60.] +@test func[1](d_string) ∈ ["forty", "eighty", "one hundred and twenty", "negative sixty"] @test func[1](sampler(d)) ∈ [40., 80., 120., -60.] +@test func[1](sampler(d_string)) ∈ ["forty", "eighty", "one hundred and twenty", "negative sixty"] @test pdf(d, -100.) == 0. +@test pdf(d_string, "negative one hundred") == 0. @test pdf(d, -100) == 0. +@test pdf(d_string, "negative one hundred") == 0. @test pdf(d, -60.) == .2 +@test pdf(d_string, "negative sixty") == .2 @test pdf(d, -60) == .2 +@test pdf(d_string, "negative sixty") == .2 @test pdf(d, 100.) == 0. +@test pdf(d_string, "one hundred") == 0. @test pdf(d, 120.) == .1 +@test pdf(d_string, "one hundred and twenty") == .1 @test cdf(d, -100.) == 0. @test cdf(d, -100) == 0. @@ -59,20 +75,28 @@ test_params(d) @test insupport(d, -60) @test insupport(d, -60.) +@test in("negative sixty", support(d_string)) @test !insupport(d, 20) @test !insupport(d, 20.) +@test !in("twenty", support(d_string)) @test insupport(d, 80) +@test in("eighty", support(d_string)) @test !insupport(d, 150) +@test !in("one hundred and fifty", support(d_string)) xs = support(d) +xs_string = support(d_string) ws = ProbabilityWeights(probs(d)) +ws_string = ProbabilityWeights(probs(d_string)) @test mean(d) ≈ mean(xs, ws) @test var(d) ≈ var(xs, ws, corrected=false) @test skewness(d) ≈ skewness(xs, ws) @test kurtosis(d) ≈ kurtosis(xs, ws) @test entropy(d) ≈ 1.2798542258336676 @test mode(d) == 40 +@test mode(d_string) == "forty" @test modes(d) == [40] +@test modes(d_string) == ["forty"] @test mgf(d, 0) ≈ 1.0 @test mgf(d, 0.17) ≈ 7.262034e7 @test cf(d, 0) ≈ 1.0 diff --git a/test/locationscale.jl b/test/locationscale.jl index d0e702c371..54aea31be1 100644 --- a/test/locationscale.jl +++ b/test/locationscale.jl @@ -39,6 +39,7 @@ function test_location_scale_normal(μ::Real, σ::Real, μD::Real, σD::Real, #### Evaluation & Sampling insupport(d,0.4) == insupport(dref,0.4) + insupport(d,0.4) == insupport(dref,Dual(0.4)) @test pdf(d,0.1) ≈ pdf(dref,0.1) @test pdf.(d,2:4) ≈ pdf.(dref,2:4) @test logpdf(d,0.4) ≈ logpdf(dref,0.4) diff --git a/test/matrixnormal.jl b/test/matrixnormal.jl index f321636e8d..4c32fc7978 100644 --- a/test/matrixnormal.jl +++ b/test/matrixnormal.jl @@ -10,6 +10,7 @@ p = 6 M = randn(n, p) U = rand(InverseWishart(n + 2, Matrix(1.0I, n, n))) +@test nsamples(InverseWishart, U) == 1 V = rand(InverseWishart(p + 2, Matrix(1.0I, p, p))) UF32 = Matrix{Float32}(U) VBF = Matrix{BigFloat}(V) @@ -138,6 +139,7 @@ end @testset "MatrixNormal sample moments" begin @test isapprox(mean(rand(D, 100000)), mean(D) , atol = 0.1) + @test nsamples(typeof(D), rand(D, 100)) == 100 @test isapprox(cov(hcat(vec.(rand(D, 100000))...)'), kron(V, U) , atol = 0.1) end diff --git a/test/mixture.jl b/test/mixture.jl index 06b69b857c..fc5ef2082a 100644 --- a/test/mixture.jl +++ b/test/mixture.jl @@ -169,8 +169,18 @@ end "rand(rng, ...)" => MersenneTwister(123)) @testset "Testing UnivariateMixture" begin + g_u32 = MixtureModel(Normal, [(0.0f0, 1.0f0), (2.0f0, 1.0f0), (-4.0f0, 1.5f0)]) + @test isa(g_u32, MixtureModel{Univariate, + <: ContinuousSupport, + <: Normal}) + @test isa(g_u32, MixtureModel{Univariate, ContinuousSupport{Float32}, + Normal{Float32}}) g_u = MixtureModel(Normal, [(0.0, 1.0), (2.0, 1.0), (-4.0, 1.5)], [0.2, 0.5, 0.3]) - @test isa(g_u, MixtureModel{Univariate, Continuous, Normal}) + @test isa(g_u, MixtureModel{Univariate, + <: ContinuousSupport, + <: Normal}) + @test isa(g_u, MixtureModel{Univariate, <: ContinuousSupport{Float64}, + Normal{Float64}}) @test ncomponents(g_u) == 3 test_mixture(g_u, 1000, 10^6, rng) test_params(g_u) diff --git a/test/multinomial.jl b/test/multinomial.jl index 6cc574567a..7791921dae 100644 --- a/test/multinomial.jl +++ b/test/multinomial.jl @@ -142,6 +142,7 @@ end # random sampling X = Matrix{Int}(undef, length(p), 100) x = func(d, X) + @test nsamples(typeof(d), x) == 100 @test x ≡ X @test isa(x, Matrix{Int}) @test all(sum(x, dims=1) .== nt) @@ -163,6 +164,7 @@ end @test x ≡ X @test all(sum.(x) .== nt) @test all(insupport(d, a) for a in x) + @test nsamples(typeof(d), x) == 100 end @testset "Testing Multinomial with $key" for (key, func) in @@ -175,6 +177,7 @@ end @test x1 ≡ X[1] @test all(sum.(x) .== nt) @test all(insupport(d, a) for a in x) + @test nsamples(typeof(d), x) == 100 end repeats = 10 @@ -195,3 +198,5 @@ nt = 10 d = Multinomial(nt, p) @test_throws DimensionMismatch rand!(d, m, false) @test_nowarn rand!(d, m) +@test Distributions.variate_form(typeof(d)) ≡ Multivariate +@test Distributions.value_support(typeof(d)) ≡ Discrete diff --git a/test/multivariate_stats.jl b/test/multivariate_stats.jl index f5bd9d8eb6..1ffcde9a1a 100644 --- a/test/multivariate_stats.jl +++ b/test/multivariate_stats.jl @@ -26,6 +26,7 @@ for d in [ dent = entropy(d) x = rand(d, n_samples) + @test nsamples(d, x) == n_samples xmean = vec(mean(x, 2)) z = x .- xmean xcov = (z * z') * (1 / n_samples) diff --git a/test/runtests.jl b/test/runtests.jl index 914f48eca9..3dcb72962f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,6 +7,7 @@ using StatsBase using LinearAlgebra const tests = [ + "dirac", "truncate", "truncnormal", "truncated_exponential", diff --git a/test/testutils.jl b/test/testutils.jl index c5901f6d5d..75c36fcb9c 100644 --- a/test/testutils.jl +++ b/test/testutils.jl @@ -25,7 +25,7 @@ end # testing the implementation of a discrete univariate distribution # -function test_distr(distr::DiscreteUnivariateDistribution, n::Int; +function test_distr(distr::CountableUnivariateDistribution, n::Int; testquan::Bool=true) test_range(distr) @@ -44,7 +44,7 @@ end # testing the implementation of a continuous univariate distribution # -function test_distr(distr::ContinuousUnivariateDistribution, n::Int; +function test_distr(distr::UnivariateDistribution{<:ContinuousSupport}, n::Int; testquan::Bool=true, rng::AbstractRNG=MersenneTwister(123)) test_range(distr) vs = get_evalsamples(distr, 0.01, 2000) @@ -73,12 +73,12 @@ end # for discrete samplers # -function test_samples(s::Sampleable{Univariate, Discrete}, # the sampleable instance - distr::DiscreteUnivariateDistribution, # corresponding distribution - n::Int; # number of samples to generate - q::Float64=1.0e-7, # confidence interval, 1 - q as confidence - verbose::Bool=false, # show intermediate info (for debugging) - rng::Union{AbstractRNG, Missing}=missing) # add an rng? +function test_samples(s::Sampleable{Univariate, <:CountableSupport}, # the sampleable instance + distr::CountableUnivariateDistribution, # corresponding distribution + n::Int; # number of samples to generate + q::Float64=1.0e-7, # confidence interval, 1 - q as confidence + verbose::Bool=false, # show intermediate info (for debugging) + rng::Union{AbstractRNG, Missing}=missing) # add an rng? # The basic idea # ------------------ @@ -105,7 +105,7 @@ function test_samples(s::Sampleable{Univariate, Discrete}, # the sampleable rmin = floor(Int,quantile(distr, 0.00001))::Int rmax = floor(Int,quantile(distr, 0.99999))::Int m = rmax - rmin + 1 # length of the range - p0 = pdf.(Ref(distr), rmin:rmax) # reference probability masses + p0 = pmf.(Ref(distr), rmin:rmax) # reference probability masses @assert length(p0) == m # determine confidence intervals for counts: @@ -145,13 +145,13 @@ function test_samples(s::Sampleable{Univariate, Discrete}, # the sampleable return samples end -test_samples(distr::DiscreteUnivariateDistribution, n::Int; +test_samples(distr::CountableUnivariateDistribution, n::Int; q::Float64=1.0e-6, verbose::Bool=false, rng=missing) = test_samples(distr, distr, n; q=q, verbose=verbose, rng=rng) # for continuous samplers # -function test_samples(s::Sampleable{Univariate, Continuous}, # the sampleable instance +function test_samples(s::Sampleable{Univariate, <: ContinuousSupport}, # the sampleable instance distr::ContinuousUnivariateDistribution, # corresponding distribution n::Int; # number of samples to generate nbins::Int=50, # divide the main interval into nbins @@ -239,7 +239,9 @@ function test_samples(s::Sampleable{Univariate, Continuous}, # the sampleable return samples end -test_samples(distr::ContinuousUnivariateDistribution, n::Int; nbins::Int=50, q::Float64=1.0e-6, verbose::Bool=false, rng=missing) = +test_samples(distr::UnivariateDistribution{<:ContinuousSupport}, + n::Int; nbins::Int=50, q::Float64=1.0e-6, verbose::Bool=false, + rng=missing) = test_samples(distr, distr, n; nbins=nbins, q=q, verbose=verbose, rng=rng) @@ -258,7 +260,7 @@ function test_range(d::UnivariateDistribution) @test isbounded(d) == (is_lb && is_ub) end -function get_evalsamples(d::DiscreteUnivariateDistribution, q::Float64) +function get_evalsamples(d::CountableUnivariateDistribution, q::Float64) # samples for testing evaluation functions (even spacing) T = eltype(d) @@ -268,7 +270,8 @@ function get_evalsamples(d::DiscreteUnivariateDistribution, q::Float64) return lv:hv end -function get_evalsamples(d::ContinuousUnivariateDistribution, q::Float64, n::Int) +function get_evalsamples(d::UnivariateDistribution{<:ContinuousSupport}, + q::Float64, n::Int) # samples for testing evaluation functions (even spacing) lv = quantile(d, q/2) @@ -280,8 +283,11 @@ end function test_support(d::UnivariateDistribution, vs::AbstractVector) for v in vs @test insupport(d, v) + @test nsamples(typeof(d), v) == 1 end + @test length(d) == 1 @test all(insupport(d, vs)) + @test nsamples(typeof(d), collect(vs)) == length(vs) if islowerbounded(d) @test isfinite(minimum(d)) @@ -297,20 +303,23 @@ function test_support(d::UnivariateDistribution, vs::AbstractVector) @test isbounded(d) == (isupperbounded(d) && islowerbounded(d)) if isbounded(d) - if isa(d, DiscreteUnivariateDistribution) + if isa(d, CountableUnivariateDistribution) s = support(d) @test isa(s, AbstractUnitRange) @test first(s) == minimum(d) @test last(s) == maximum(d) end end + if isa(d, UnivariateDistribution{<:ContinuousSupport}) + @test support(d) == RealInterval(minimum(d), maximum(d)) + end end #### Testing evaluation methods -function test_range_evaluation(d::DiscreteUnivariateDistribution) - # check the consistency between range-based and ordinary pdf +function test_range_evaluation(d::CountableUnivariateDistribution) + # check the consistency between range-based and ordinary pmf vmin = minimum(d) vmax = maximum(d) @test vmin <= vmax @@ -324,24 +333,25 @@ function test_range_evaluation(d::DiscreteUnivariateDistribution) rmin = round(Int, islowerbounded(d) ? vmin : quantile(d, 0.001))::Int rmax = round(Int, isupperbounded(d) ? vmax : quantile(d, 0.999))::Int - p0 = pdf.(Ref(d), collect(rmin:rmax)) - @test pdf.(Ref(d), rmin:rmax) ≈ p0 + p0 = pmf.(Ref(d), collect(rmin:rmax)) + @test pmf.(Ref(d), rmin:rmax) ≈ p0 if rmin + 2 <= rmax - @test pdf.(Ref(d), rmin+1:rmax-1) ≈ p0[2:end-1] + @test pmf.(Ref(d), rmin+1:rmax-1) ≈ p0[2:end-1] end if isbounded(d) - @test pdf.(Ref(d), support(d)) ≈ p0 - @test pdf.(Ref(d), rmin-2:rmax) ≈ vcat(0.0, 0.0, p0) - @test pdf.(Ref(d), rmin:rmax+3) ≈ vcat(p0, 0.0, 0.0, 0.0) - @test pdf.(Ref(d), rmin-2:rmax+3) ≈ vcat(0.0, 0.0, p0, 0.0, 0.0, 0.0) + @test pmf.(Ref(d), support(d)) ≈ p0 + @test pmf.(Ref(d), rmin-2:rmax) ≈ vcat(0.0, 0.0, p0) + @test pmf.(Ref(d), rmin:rmax+3) ≈ vcat(p0, 0.0, 0.0, 0.0) + @test pmf.(Ref(d), rmin-2:rmax+3) ≈ vcat(0.0, 0.0, p0, 0.0, 0.0, 0.0) elseif islowerbounded(d) - @test pdf.(Ref(d), rmin-2:rmax) ≈ vcat(0.0, 0.0, p0) + @test pmf.(Ref(d), rmin-2:rmax) ≈ vcat(0.0, 0.0, p0) end end -function test_evaluation(d::DiscreteUnivariateDistribution, vs::AbstractVector, testquan::Bool=true) +function test_evaluation(d::CountableUnivariateDistribution, + vs::AbstractVector, testquan::Bool=true) nv = length(vs) p = Vector{Float64}(undef, nv) c = Vector{Float64}(undef, nv) @@ -352,10 +362,10 @@ function test_evaluation(d::DiscreteUnivariateDistribution, vs::AbstractVector, ci = 0. for (i, v) in enumerate(vs) - p[i] = pdf(d, v) + p[i] = pmf(d, v) c[i] = cdf(d, v) cc[i] = ccdf(d, v) - lp[i] = logpdf(d, v) + lp[i] = logpmf(d, v) lc[i] = logcdf(d, v) lcc[i] = logccdf(d, v) @@ -383,17 +393,18 @@ function test_evaluation(d::DiscreteUnivariateDistribution, vs::AbstractVector, end # consistency of scalar-based and vectorized evaluation - @test pdf.(Ref(d), vs) ≈ p + @test pmf.(Ref(d), vs) ≈ p @test cdf.(Ref(d), vs) ≈ c @test ccdf.(Ref(d), vs) ≈ cc - @test logpdf.(Ref(d), vs) ≈ lp + @test logpmf.(Ref(d), vs) ≈ lp @test logcdf.(Ref(d), vs) ≈ lc @test logccdf.(Ref(d), vs) ≈ lcc end -function test_evaluation(d::ContinuousUnivariateDistribution, vs::AbstractVector, testquan::Bool=true) +function test_evaluation(d::UnivariateDistribution{<:ContinuousSupport}, + vs::AbstractVector, testquan::Bool=true) nv = length(vs) p = Vector{Float64}(undef, nv) c = Vector{Float64}(undef, nv) @@ -429,10 +440,14 @@ 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)) + @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)) end end end @@ -464,11 +479,11 @@ end #### Testing statistics methods -function test_stats(d::DiscreteUnivariateDistribution, vs::AbstractVector) +function test_stats(d::CountableUnivariateDistribution, vs::AbstractVector) # using definition (or an approximation) vf = Float64[v for v in vs] - p = pdf.(Ref(d), vf) + p = pmf.(Ref(d), vf) xmean = dot(p, vf) xvar = dot(p, abs2.(vf .- xmean)) xstd = sqrt(xvar) @@ -502,7 +517,8 @@ allow_test_stats(d::UnivariateDistribution) = true allow_test_stats(d::NoncentralBeta) = false allow_test_stats(::StudentizedRange) = false -function test_stats(d::ContinuousUnivariateDistribution, xs::AbstractVector{Float64}) +function test_stats(d::UnivariateDistribution{<:ContinuousSupport}, + xs::AbstractVector{Float64}) # using Monte Carlo methods if !(isfinite(mean(d)) && isfinite(var(d))) @@ -511,6 +527,7 @@ function test_stats(d::ContinuousUnivariateDistribution, xs::AbstractVector{Floa vd = var(d) n = length(xs) + @test length(d) == 1 xmean = mean(xs) xvar = var(xs) xstd = sqrt(xvar)