Algorithms & Methods

ClassicEM<:AbstractEM

The EM algorithm was introduced by A. P. Dempster, N. M. Laird and D. B. Rubin in 1977 in the reference paper Maximum Likelihood from Incomplete Data Via the EM Algorithm.

Base.@kwdef struct StochasticEM<:AbstractEM
    rng::AbstractRNG = Random.GLOBAL_RNG
end

The Stochastic EM algorithm was introduced by G. Celeux, and J. Diebolt. in 1985 in The SEM Algorithm: A probabilistic teacher algorithm derived from the EM algorithm for the mixture problem.

The default random seed is Random.GLOBAL_RNG but it can be changed via StochasticEM(seed).

fit_mle(mix::MixtureModel, y::AbstractVecOrMat, weights...; method = ClassicEM(), display=:none, maxiter=1000, atol=1e-3, rtol=nothing, robust=false, infos=false)

Use the an Expectation Maximization (EM) algorithm to maximize the Loglikelihood (fit) the mixture with an i.i.d sample y. The mix input is a mixture that is used to initilize the EM algorithm.

• weights when provided, it will compute a weighted version of the EM. (Useful for fitting mixture of mixtures)
• method determines the algorithm used.
• infos = true returns a Dict with informations on the algorithm (converged, iteration number, loglikelihood).
• robust = true will prevent the (log)likelihood to overflow to -∞ or ∞.
• atol criteria determining the convergence of the algorithm. If the Loglikelihood difference between two iteration i and i+1 is smaller than atol i.e. |ℓ⁽ⁱ⁺¹⁾ - ℓ⁽ⁱ⁾|<atol, the algorithm stops.
• rtol relative tolerance for convergence, |ℓ⁽ⁱ⁺¹⁾ - ℓ⁽ⁱ⁾|<rtol*(|ℓ⁽ⁱ⁺¹⁾| + |ℓ⁽ⁱ⁾|)/2 (does not check if rtol is nothing)
• display value can be :none, :iter, :final to display Loglikelihood evolution at each iterations :iter or just the final one :final

fit_mle(mix::AbstractArray{<:MixtureModel}, y::AbstractVecOrMat, weights...; method = ClassicEM(), display=:none, maxiter=1000, atol=1e-3, rtol=nothing, robust=false, infos=false)

Do the same as fit_mle for each (initial) mixtures in the mix array. Then it selects the one with the largest loglikelihood. Warning: It uses try and catch to avoid errors messages in case EM converges toward a singular solution (probably using robust should be enough in most case to avoid errors).

predict(mix::MixtureModel, y::AbstractVector; robust=false)

Evaluate the most likely category for each observations given a MixtureModel.

• robust = true will prevent the (log)likelihood to overflow to -∞ or ∞.

predict_proba(mix::MixtureModel, y::AbstractVecOrMat; robust=false)

Evaluate the probability for each observations to belong to a category given a MixtureModel..

• robust = true will prevent the (log)likelihood to under(overflow)flow to -∞ (or ∞).

In EM.jl the "instance" version of fit_mle is supported (in addition of the current "type" version). Note this is not supported in Distributions.jl. Example: fit_mle(Bernoulli(0.2), x) is accepted in addition of fit_mle(Bernoulli, x) this allows compatibility with how fit_mle(g::Product) and fit_mle(g::MixtureModel) are written.

fit_mle(g::Product, x::AbstractMatrix)
fit_mle(g::Product, x::AbstractMatrix, γ::AbstractVector)

The fit_mle for multivariate Product distributions g is the product_distribution of fit_mle of each components of g. Product is meant to be depreacated in next versions of Distribution.jl. Use the analog VectorOfUnivariateDistribution type instead.

fit_mle(dists::ArrayOfUnivariateDistribution, x::AbstractArray)
fit_mle(dists::ArrayOfUnivariateDistribution, x::AbstractArray, γ::AbstractVector)

The fit_mle for a ArrayOfUnivariateDistribution distributions dists is the product_distribution of fit_mle of each components of dists. VectorOfUnivariateDistribution should act like old Product while ArrayOfUnivariateDistribution are not really tested yet.

fit_mle(::Type{<:Dirac}, x::AbstractArray{<:Real}[, w::AbstractArray{<:Real}])

fit_mle for Dirac distribution (weighted or not) data sets.

fit_mle(::Type{<:Laplace}, x::AbstractArray{<:Real}, w::AbstractArray{<:Real})

fit_mle for Laplace distribution weighted data sets.

fit_mle(::Type{<:Uniform}, x::AbstractArray{<:Real}, w::AbstractArray{<:Real})

fit_mle for Uniform distribution weighted data sets. It is just the same as unweigted (removing zero weighted data).

fit_mle!(α::AbstractVector, dists::AbstractVector{F} where {F<:Distribution}, y::AbstractVecOrMat, method::ClassicEM; display=:none, maxiter=1000, atol=1e-3, rtol=nothing, robust=false)

Use the EM algorithm to update the Distribution dists and weights α composing a mixture distribution.

• robust = true will prevent the (log)likelihood to overflow to -∞ or ∞.
• atol criteria determining the convergence of the algorithm. If the Loglikelihood difference between two iteration i and i+1 is smaller than atol i.e. |ℓ⁽ⁱ⁺¹⁾ - ℓ⁽ⁱ⁾|<atol, the algorithm stops.
• rtol relative tolerance for convergence, |ℓ⁽ⁱ⁺¹⁾ - ℓ⁽ⁱ⁾|<rtol*(|ℓ⁽ⁱ⁺¹⁾| + |ℓ⁽ⁱ⁾|)/2 (does not check if rtol is nothing)
• display value can be :none, :iter, :final to display Loglikelihood evolution at each iterations :iter or just the final one :final

fit_mle!(α::AbstractVector, dists::AbstractVector{F} where {F<:Distribution}, y::AbstractVecOrMat, method::StochasticEM; display=:none, maxiter=1000, atol=1e-3, robust=false)

Use the stochastic EM algorithm to update the Distribution dists and weights α composing a mixture distribution.

• robust = true will prevent the (log)likelihood to overflow to -∞ or ∞.
• atol criteria determining the convergence of the algorithm. If the Loglikelihood difference between two iteration i and i+1 is smaller than atol i.e. |ℓ⁽ⁱ⁺¹⁾ - ℓ⁽ⁱ⁾|<atol, the algorithm stops.
• rtol relative tolerance for convergence, |ℓ⁽ⁱ⁺¹⁾ - ℓ⁽ⁱ⁾|<rtol*(|ℓ⁽ⁱ⁺¹⁾| + |ℓ⁽ⁱ⁾|)/2 (does not check if rtol is nothing)
• display value can be :none, :iter, :final to display Loglikelihood evolution at each iterations :iter or just the final one :final

M_step!(α, dists, y, γ, method::ClassicEM)

For the ClassicEM the weigths γ computed at E-step for each observation in y are used to update α and dists.

M_step!(α, dists, y, cat, method::StochasticEM)

For the StochasticEM the cat drawn at S-step for each observation in y is used to update α and dists.