Add fit_mle methods

docs/Project.toml

@@ -1,4 +1,5 @@
 [deps]
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
 

docs/src/fit.md

@@ -42,6 +42,7 @@ The `fit_mle` method has been implemented for the following distributions:
 - [`Binomial`](@ref)
 - [`Categorical`](@ref)
 - [`DiscreteUniform`](@ref)
+- [`Dirac`](@ref)
 - [`Exponential`](@ref)
 - [`LogNormal`](@ref)
 - [`Normal`](@ref)
@@ -60,6 +61,7 @@ The `fit_mle` method has been implemented for the following distributions:
 - [`Multinomial`](@ref)
 - [`MvNormal`](@ref)
 - [`Dirichlet`](@ref)
+- [`Product`](@ref)
 
 For most of these distributions, the usage is as described above. For a few special distributions that require additional information for estimation, we have to use a modified interface:
 
@@ -74,6 +76,21 @@ fit_mle(Categorical, x)        # equivalent to fit_mle(Categorical, max(x), x)
 fit_mle(Categorical, x, w)
 ```
 
+It is also possible to directly input a distribution `fit_mle(d::Distribution, x[, w])`. This form avoids the extra arguments:
+```julia
+fit_mle(Binomial(n, 0.1), x) 
+# equivalent to fit_mle(Binomial, ntrials(Binomial(n, 0.1)), x), here the parameter 0.1 is not used
+
+fit_mle(Categorical(p), x) 
+# equivalent to fit_mle(Categorical, ncategories(Categorical(p)), x), here the only the length of p is used not its values
+
+d = product_distribution([Exponential(0.5), Normal(11.3, 3.2)])
+fit_mle(d, x) 
+# equivalent to product_distribution([fit_mle(Exponential, x[1,:]), fit_mle(Normal, x[2, :])]). Again parameters of d are not used.
+```
+Note that for standard distributions, the values of the distribution parameters `d` are not used in `fit_mle` only the “structure” of `d` is passed into `fit_mle`.
+However, for complex Maximum Likelihood estimation requiring optimization, e.g., EM algorithm, one could use `D` as an initial guess.
+
 ## Sufficient Statistics
 
 For many distributions, the estimation can be based on (sum of) sufficient statistics computed from a dataset. To simplify implementation, for such distributions, we implement `suffstats` method instead of `fit_mle` directly:

docs/src/mixture.md

@@ -106,4 +106,4 @@ rand!(::AbstractMixtureModel, ::AbstractArray)
 ## Estimation
 
 There are several methods for the estimation of mixture models from data, and this problem remains an open research topic.
-This package does not provide facilities for estimating mixture models. One can resort to other packages, *e.g.* [*GaussianMixtures.jl*](https://github.com/davidavdav/GaussianMixtures.jl), for this purpose.
+This package does not provide facilities for estimating mixture models. One can resort to other packages, *e.g.* [*GaussianMixtures.jl*](https://github.com/davidavdav/GaussianMixtures.jl) or [*ExpectationMaximization.jl*](https://github.com/dmetivie/ExpectationMaximization.jl), for this purpose.

src/genericfit.jl

@@ -1,6 +1,6 @@
 # generic functions for distribution fitting
 
-function suffstats(dt::Type{D}, xs...) where D<:Distribution
+function suffstats(dt::Type{D}, xs...) where {D<:Distribution}
     argtypes = tuple(D, map(typeof, xs)...)
     error("suffstats is not implemented for $argtypes.")
 end
@@ -30,5 +30,7 @@ fit_mle(dt::Type{D}, x::AbstractArray, w::AbstractArray) where {D<:UnivariateDis
 fit_mle(dt::Type{D}, x::AbstractMatrix) where {D<:MultivariateDistribution} = fit_mle(D, suffstats(D, x))
 fit_mle(dt::Type{D}, x::AbstractMatrix, w::AbstractArray) where {D<:MultivariateDistribution} = fit_mle(D, suffstats(D, x, w))
 
+fit_mle(dt::D, args...) where {D<:Distribution} = fit_mle(typeof(dt), args...)
+
 fit(dt::Type{D}, x) where {D<:Distribution} = fit_mle(D, x)
 fit(dt::Type{D}, args...) where {D<:Distribution} = fit_mle(D, args...)

src/multivariate/product.jl

@@ -27,15 +27,17 @@ function Product(v::V) where {S<:ValueSupport,T<:UnivariateDistribution{S},V<:Ab
         "`Product(v)` is deprecated, please use `product_distribution(v)`",
         :Product,
     )
-    return Product{S, T, V}(v)
+    return Product{S,T,V}(v)
 end
 
 length(d::Product) = length(d.v)
 function Base.eltype(::Type{<:Product{S,T}}) where {S<:ValueSupport,
-                                                    T<:UnivariateDistribution{S}}
+    T<:UnivariateDistribution{S}}
     return eltype(T)
 end
 
+params(g::Product) = params.(g.v)
+
 _rand!(rng::AbstractRNG, d::Product, x::AbstractVector{<:Real}) =
     map!(Base.Fix1(rand, rng), x, d.v)
 function _logpdf(d::Product, x::AbstractVector{<:Real})
@@ -60,3 +62,17 @@ maximum(d::Product) = map(maximum, d.v)
 function product_distribution(dists::V) where {S<:ValueSupport,T<:UnivariateDistribution{S},V<:AbstractVector{T}}
     return Product{S,T,V}(dists)
 end
+
+#### Fitting
+
+"""
+    fit_mle(g::Product, x::AbstractMatrix)
+    fit_mle(g::Product, x::AbstractMatrix, γ::AbstractVector)
+
+The `fit_mle` for a multivariate Product distributions `g` is the `product_distribution` of `fit_mle` of each components of `g`.
+"""
+function fit_mle(g::Product, x::AbstractMatrix, args...)
+    d = size(x, 1)
+    length(g) == d || throw(DimensionMismatch("The dimensions of g and x are inconsistent."))
+    return product_distribution([fit_mle(g.v[s], y, args...) for (s, y) in enumerate(eachrow(x))])
+end

src/univariate/continuous/laplace.jl

@@ -136,3 +136,15 @@ function fit_mle(::Type{<:Laplace}, x::AbstractArray{<:Real})
     xc .= abs.(x .- m)
     return Laplace(m, mean(xc))
 end
+
+function fit_mle(::Type{<:Laplace}, x::AbstractArray{<:Real}, w::AbstractArray{<:Real})
+    n = length(x)
+    if n != length(w)
+        throw(DimensionMismatch("Inconsistent array lengths."))
+    end
+    xc = similar(x)
+    copyto!(xc, x)
+    m = median(xc, weights(w)) # https://github.com/JuliaStats/StatsBase.jl/blob/0b64a9c8a16355da16469d0fe5016e0fdf099af5/src/weights.jl#L788
+    xc .= abs.(x .- m)
+    return Laplace(m, mean(xc, weights(w)))
+end

src/univariate/discrete/binomial.jl

@@ -196,6 +196,7 @@ suffstats(::Type{T}, data::BinomData, w::AbstractArray{<:Real}) where {T<:Binomi
 
 fit_mle(::Type{T}, ss::BinomialStats) where {T<:Binomial} = T(ss.n, ss.ns / (ss.ne * ss.n))
 
+fit_mle(d::T, x::AbstractArray{<:Integer}) where {T<:Binomial} = fit_mle(T, suffstats(T, ntrials(d), x))
 fit_mle(::Type{T}, n::Integer, x::AbstractArray{<:Integer}) where {T<:Binomial}= fit_mle(T, suffstats(T, n, x))
 fit_mle(::Type{T}, n::Integer, x::AbstractArray{<:Integer}, w::AbstractArray{<:Real}) where {T<:Binomial} = fit_mle(T, suffstats(T, n, x, w))
 fit_mle(::Type{T}, data::BinomData) where {T<:Binomial} = fit_mle(T, suffstats(T, data))

src/univariate/discrete/categorical.jl

@@ -26,7 +26,7 @@ External links:
 """
 const Categorical{P<:Real,Ps<:AbstractVector{P}} = DiscreteNonParametric{Int,P,Base.OneTo{Int},Ps}
 
-function Categorical{P,Ps}(p::Ps; check_args::Bool=true) where {P<:Real, Ps<:AbstractVector{P}}
+function Categorical{P,Ps}(p::Ps; check_args::Bool=true) where {P<:Real,Ps<:AbstractVector{P}}
     @check_args Categorical (p, isprobvec(p), "vector p is not a probability vector")
     return Categorical{P,Ps}(Base.OneTo(length(p)), p; check_args=check_args)
 end
@@ -36,7 +36,7 @@ Categorical(p::AbstractVector{P}; check_args::Bool=true) where {P<:Real} =
 
 function Categorical(k::Integer; check_args::Bool=true)
     @check_args Categorical (k, k >= 1, "at least one category is required")
-    return Categorical{Float64,Vector{Float64}}(Base.OneTo(k), fill(1/k, k); check_args=false)
+    return Categorical{Float64,Vector{Float64}}(Base.OneTo(k), fill(1 / k, k); check_args=false)
 end
 
 Categorical(probabilities::Real...; check_args::Bool=true) = Categorical([probabilities...]; check_args=check_args)
@@ -49,7 +49,7 @@ convert(::Type{Categorical{P,Ps}}, x::AbstractVector{<:Real}) where {
 ### 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}) where {P<:Real,Ps<:AbstractVector{P}} = (probs(d),)
 partype(::Categorical{T}) where {T<:Real} = T
 
 ### Statistics
@@ -59,7 +59,7 @@ function median(d::Categorical{T}) where {T<:Real}
     p = probs(d)
     cp = zero(T)
     i = 0
-    while cp < 1/2 && i <= k
+    while cp < 1 / 2 && i <= k
         i += 1
         @inbounds cp += p[i]
     end
@@ -87,16 +87,16 @@ function _pdf!(r::AbstractArray, d::Categorical{T}, rgn::UnitRange) where {T<:Re
     vr = min(vlast, ncategories(d))
     p = probs(d)
     if vl > vfirst
-        for i = 1:(vl - vfirst)
+        for i = 1:(vl-vfirst)
             r[i] = zero(T)
         end
     end
     fm1 = vfirst - 1
     for v = vl:vr
-        r[v - fm1] = p[v]
+        r[v-fm1] = p[v]
     end
     if vr < vlast
-        for i = (vr - vfirst + 2):length(rgn)
+        for i = (vr-vfirst+2):length(rgn)
             r[i] = zero(T)
         end
     end
@@ -107,7 +107,7 @@ end
 # sampling
 
 sampler(d::Categorical{P,Ps}) where {P<:Real,Ps<:AbstractVector{P}} =
-   AliasTable(probs(d))
+    AliasTable(probs(d))
 
 
 ### sufficient statistics
@@ -116,36 +116,36 @@ struct CategoricalStats <: SufficientStats
     h::Vector{Float64}
 end
 
-function add_categorical_counts!(h::Vector{Float64}, x::AbstractArray{T}) where T<:Integer
-    for i = 1 : length(x)
+function add_categorical_counts!(h::Vector{Float64}, x::AbstractArray{T}) where {T<:Integer}
+    for i = 1:length(x)
         @inbounds xi = x[i]
-        h[xi] += 1.   # cannot use @inbounds, as no guarantee that x[i] is in bound
+        h[xi] += 1.0   # cannot use @inbounds, as no guarantee that x[i] is in bound
     end
     h
 end
 
-function add_categorical_counts!(h::Vector{Float64}, x::AbstractArray{T}, w::AbstractArray{Float64}) where T<:Integer
+function add_categorical_counts!(h::Vector{Float64}, x::AbstractArray{T}, w::AbstractArray{Float64}) where {T<:Integer}
     n = length(x)
     if n != length(w)
         throw(DimensionMismatch("Inconsistent array lengths."))
     end
-    for i = 1 : n
+    for i = 1:n
         @inbounds xi = x[i]
         @inbounds wi = w[i]
         h[xi] += wi   # cannot use @inbounds, as no guarantee that x[i] is in bound
     end
     h
 end
 
-function suffstats(::Type{<:Categorical}, k::Int, x::AbstractArray{T}) where T<:Integer
+function suffstats(::Type{<:Categorical}, k::Int, x::AbstractArray{T}) where {T<:Integer}
     CategoricalStats(add_categorical_counts!(zeros(k), x))
 end
 
-function suffstats(::Type{<:Categorical}, k::Int, x::AbstractArray{T}, w::AbstractArray{Float64}) where T<:Integer
+function suffstats(::Type{<:Categorical}, k::Int, x::AbstractArray{T}, w::AbstractArray{Float64}) where {T<:Integer}
     CategoricalStats(add_categorical_counts!(zeros(k), x, w))
 end
 
-const CategoricalData = Tuple{Int, AbstractArray}
+const CategoricalData = Tuple{Int,AbstractArray}
 
 suffstats(::Type{<:Categorical}, data::CategoricalData) = suffstats(Categorical, data...)
 suffstats(::Type{<:Categorical}, data::CategoricalData, w::AbstractArray{Float64}) = suffstats(Categorical, data..., w)
@@ -156,11 +156,12 @@ function fit_mle(::Type{<:Categorical}, ss::CategoricalStats)
     Categorical(normalize!(ss.h, 1))
 end
 
-function fit_mle(::Type{<:Categorical}, k::Integer, x::AbstractArray{T}) where T<:Integer
+fit_mle(d::T, x::AbstractArray{<:Integer}) where {T<:Categorical} = fit_mle(T, ncategories(d), x)
+function fit_mle(::Type{<:Categorical}, k::Integer, x::AbstractArray{T}) where {T<:Integer}
     Categorical(normalize!(add_categorical_counts!(zeros(k), x), 1), check_args=false)
 end
 
-function fit_mle(::Type{<:Categorical}, k::Integer, x::AbstractArray{T}, w::AbstractArray{Float64}) where T<:Integer
+function fit_mle(::Type{<:Categorical}, k::Integer, x::AbstractArray{T}, w::AbstractArray{Float64}) where {T<:Integer}
     Categorical(normalize!(add_categorical_counts!(zeros(k), x, w), 1), check_args=false)
 end
 

src/univariate/discrete/dirac.jl

@@ -50,9 +50,21 @@ logccdf(d::Dirac, x::Real) = x < d.value ? 0.0 : isnan(x) ? NaN : -Inf
 quantile(d::Dirac{T}, p::Real) where {T} = 0 <= p <= 1 ? d.value : T(NaN)
 
 mgf(d::Dirac, t) = exp(t * d.value)
-cgf(d::Dirac, t) = t*d.value
+cgf(d::Dirac, t) = t * d.value
 cf(d::Dirac, t) = cis(t * d.value)
 
 #### Sampling
 
 rand(rng::AbstractRNG, d::Dirac) = d.value
+
+#### Fitting
+
+fit_mle(::Type{<:Dirac}, x::AbstractArray{T}) where {T<:Real} = length(unique(x)) == 1 ? Dirac(first(x)) : Dirac(NaN)
+
+function fit_mle(::Type{<:Dirac}, x::AbstractArray{T}, w::AbstractArray{Float64}) where {T<:Real}
+    n = length(x)
+    if n != length(w)
+        throw(DimensionMismatch("Inconsistent array lengths."))
+    end
+    return length(unique(x[findall(!iszero, w)])) == 1 ? Dirac(first(x)) : Dirac(NaN)
+end