Ce serveur Gitlab sera éteint le 30 juin 2020, pensez à migrer vos projets vers les serveurs gitlab-research.centralesupelec.fr et gitlab-student.centralesupelec.fr !

Commit 273b0eaa authored by Dahua Lin's avatar Dahua Lin

Merge pull request #303 from JuliaStats/dh/mixture2

Reimplement Mixture models
parents 95cee4da b0b951d6
......@@ -39,6 +39,8 @@ export
ContinuousMatrixDistribution,
SufficientStats,
AbstractMvNormal,
UnivariateMixture,
MultivariateMixture,
# distribution types
Arcsine,
......@@ -141,6 +143,7 @@ export
cquantile, # complementary quantile (i.e. using prob in right hand tail)
cumulant, # cumulants of distribution
complete, # turn an incomplete formulation into a complete distribution
components, # get components from a mixture model
concentration, # the concentration parameter
dim, # sample dimension of multivariate distribution
entropy, # entropy of distribution in nats
......@@ -195,6 +198,7 @@ export
ncategories, # the number of categories in a Categorical distribution
pdf, # probability density function (ContinuousDistribution)
pmf, # probability mass function (DiscreteDistribution)
priorprobs, # prior probabilities
quantile, # inverse of cdf (defined for p in (0,1))
qqbuild, # build a paired quantiles data structure for qqplots
sampler, # create a Sampler object for efficient samples
......
immutable MixtureModel{VF,VS,Component<:Distribution} <: Distribution{VF,VS}
# finite mixture models
immutable MixtureModel{VF<:VariateForm, VS<:ValueSupport, Component<:Distribution} <: Distribution{VF,VS}
components::Vector{Component}
probs::Vector{Float64}
aliastable::AliasTable
function MixtureModel(c::Vector{Component}, p::Vector{Float64})
if !(Component <: Distribution{VF,VS})
throw(TypeError(:MixtureModel,
"Mixture components type mismatch",
Distribution{VF,VS},
Component))
end
if length(c) != length(p)
throw(ArgumentError(string("components and probs must have",
" the same number of elements")))
end
sizes = [size(component)::Tuple for component in c]
if !all(sizes .== sizes[1])
error("MixtureModel: mixture components have different dimensions")
end
sump = 0.0
for i in 1:length(p)
if p[i] < 0.0
error("MixtureModel: probabilities must be non-negative")
end
sump += p[i]
end
table = AliasTable(p ./ sump)
new(c, p ./ sump, table)
end
prior::Categorical
end
function MixtureModel{C<:Distribution}(cs::Vector{C}, ps::Vector{Float64})
typealias UnivariateMixture MixtureModel{Univariate}
typealias MultivariateMixture MixtureModel{Multivariate}
typealias MatrixvariateMixture MixtureModel{Matrixvariate}
#### Constructors
function MixtureModel{C<:Distribution}(components::Vector{C}, prior::Categorical)
length(components) == ncategories(prior) ||
error("Inconsistent sizes of components and prior.")
VF = variate_form(C)
VS = value_support(C)
MixtureModel{VF,VS,C}(cs, ps)
MixtureModel{VF,VS,C}(components, prior)
end
dim(d::MixtureModel{Multivariate}) = dim(d.components[1])
MixtureModel{C<:Distribution}(components::Vector{C}, p::Vector{Float64}) =
MixtureModel(components, Categorical(p))
# all components have the same prior probabilities
MixtureModel{C<:Distribution}(components::Vector{C}) =
MixtureModel(components, Categorical(length(components)))
#### Basic properties
length(d::MultivariateMixture) = length(d.components[1])
size(d::MatrixvariateMixture) = size(d.components[1])
components(d::MixtureModel) = d.components
priorprobs(d::MixtureModel) = d.prior.prob
function mean(d::MixtureModel)
function mean(d::UnivariateMixture)
cs = components(d)
p = priorprobs(d)
m = 0.0
for i in 1:length(d.components)
if d.probs[i] > 0.
m += mean(d.components[i]) * d.probs[i]
for i = 1:length(cs)
pi = p[i]
if pi > 0.0
m += mean(cs[i]) * pi
end
end
return m
end
function _mixpdf(d::MixtureModel, x)
p = 0.0
for i in 1:length(d.components)
p += pdf(d.components[i], x) * d.probs[i]
function mean(d::MultivariateMixture)
cs = components(d)
p = priorprobs(d)
m = zeros(length(d))
for i = 1:length(cs)
pi = p[i]
if pi > 0.0
BLAS.axpy!(pi, mean(cs[i]), m)
end
end
return p
return m
end
function _mixlogpdf(d::MixtureModel, x)
l = [(logpdf(d.components[i],x)+log(d.probs[i]))::Float64
for i in find(d.probs .> 0.)]
m = maximum(l)
log(sum(exp(l.-m))) + m
function var(d::UnivariateMixture)
cs = components(d)
p = priorprobs(d)
K = length(cs)
means = Array(Float64, K)
m = 0.0
v = 0.0
for i = 1:K
pi = p[i]
if pi > 0.0
ci = cs[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
pdf(d::MixtureModel{Univariate}, x::Real) = _mixpdf(d, x)
logpdf(d::MixtureModel{Univariate}, x::Real) = _mixlogpdf(d, x)
_pdf(d::MixtureModel{Multivariate}, x::AbstractVector) = _mixpdf(d, x)
_logpdf(d::MixtureModel{Multivariate}, x::AbstractVector) = _mixpdf(d, x)
#### Evaluation
function rand(d::MixtureModel)
i = rand(d.aliastable)
return rand(d.components[i])
function _mixpdf1(d::MixtureModel, x)
cs = components(d)
p = priorprobs(d)
v = 0.0
for i = 1:length(cs)
pi = p[i]
if pi > 0.0
v += pdf(cs[i], x) * pi
end
end
return v
end
function var(d::MixtureModel)
m = 0.0
squared_mean_mixture = mean(d).^2
for i in 1:length(d.components)
if d.probs[i] > 0.
m += (var(d.components[i]) .- squared_mean_mixture .+ mean(d.components[i]).^2) * d.probs[i]
function _mixpdf!(r::DenseArray, d::MixtureModel, x)
cs = components(d)
p = priorprobs(d)
fill!(r, 0.0)
t = Array(Float64, size(r))
for i = 1:length(cs)
pi = p[i]
if pi > 0.0
# compute pdf in batch
pdf!(t, cs[i], x)
# accumulate pdf in batch
BLAS.axpy!(pi, t, r)
end
end
return m
return r
end
size(d::MixtureModel) = size(d.components[1])
length(d::MixtureModel) = length(d.components[1])
function _mixlogpdf1(d::MixtureModel, x)
# using the formula below for numerical stability
#
# logpdf(d, x) = log(sum_i pri[i] * pdf(cs[i], x))
# = log(sum_i pri[i] * exp(logpdf(cs[i], x)))
# = log(sum_i exp(logpri[i] + logpdf(cs[i], x)))
# = m + log(sum_i exp(logpri[i] + logpdf(cs[i], x) - m))
#
# m is chosen to be the maximum of logpri[i] + logpdf(cs[i], x) - m
# such that the argument of exp is in a reasonable range
#
cs = components(d)
p = priorprobs(d)
K = length(cs)
lp = Array(Float64, K)
m = -Inf # m <- the maximum of log(p(cs[i], x)) + log(pri[i])
for i = 1:K
pi = p[i]
if pi > 0.0
# lp[i] <- log(p(cs[i], x)) + log(pri[i])
lp[i] = lp_i = logpdf(cs[i], x) + log(pi)
if lp_i > m
m = lp_i
end
end
end
v = 0.0
for i = 1:K
if p[i] > 0.0
v += exp(lp[i] - m)
end
end
return m + log(v)
end
function _mixlogpdf!(r::DenseArray, d::MixtureModel, x)
cs = components(d)
p = priorprobs(d)
K = length(cs)
n = length(r)
Lp = Array(Float64, n, K)
m = fill(-Inf, n)
for i = 1:K
pi = p[i]
if pi > 0.0
lpri = log(pi)
lp_i = view(Lp, :, i)
# compute logpdf in batch and store
logpdf!(lp_i, cs[i], x)
# in the mean time, add log(prior) to lp and
# update the maximum for each sample
for j = 1:n
lp_i[j] += lpri
if lp_i[j] > m[j]
m[j] = lp_i[j]
end
end
end
end
fill!(r, 0.0)
for i = 1:K
if p[i] > 0.0
lp_i = view(Lp, :, i)
for j = 1:n
r[j] += exp(lp_i[j] - m[j])
end
end
end
for j = 1:n
r[j] = log(r[j]) + m[j]
end
return r
end
pdf(d::UnivariateMixture, x::Real) = _mixpdf1(d, x)
logpdf(d::UnivariateMixture, x::Real) = _mixlogpdf1(d, x)
_pdf!(r::AbstractArray, d::UnivariateMixture, x::DenseArray) = _mixpdf!(r, d, x)
_logpdf!(r::AbstractArray, d::UnivariateMixture, x::DenseArray) = _mixlogpdf!(r, d, x)
_pdf(d::MultivariateMixture, x::AbstractVector) = _mixpdf1(d, x)
_logpdf(d::MultivariateMixture, x::AbstractVector) = _mixlogpdf1(d, x)
_pdf!(r::AbstractArray, d::MultivariateMixture, x::DenseMatrix) = _mixpdf!(r, d, x)
_lodpdf!(r::AbstractArray, d::MultivariateMixture, x::DenseMatrix) = _mixlogpdf!(r, d, x)
## Sampling
immutable MixtureSampler{VF,VS,Sampler} <: Sampleable{VF,VS}
csamplers::Vector{Sampler}
psampler::AliasTable
end
function MixtureSampler{VF,VS}(d::MixtureModel{VF,VS})
csamplers = map(sampler, d.components)
psampler = sampler(d.prior)
MixtureSampler{VF,VS,eltype(csamplers)}(csamplers, psampler)
end
rand(d::MixtureModel) = rand(d.components[rand(d.prior)])
rand(s::MixtureSampler) = rand(s.csamplers[rand(s.psampler)])
_rand!(s::MixtureSampler{Multivariate}, x::DenseVector) = _rand!(s.csamplers[rand(s.psampler)], x)
sampler(d::MixtureModel) = MixtureSampler(d)
......@@ -126,8 +126,15 @@ gradlogpdf(d::MvNormal, x::Vector{Float64}) = -(d.Σ \ (x - d.μ))
# Sampling (for GenericMvNormal)
_rand!(d::MvNormal, x::DenseVecOrMat{Float64}) = add!(unwhiten!(d.Σ, randn!(x)), d.μ)
_rand!(d::MvNormal, x::VecOrMat{Float64}) = add!(unwhiten!(d.Σ, randn!(x)), d.μ)
# Workaround: randn! only works for Array, but not generally for DenseArray
function _rand!(d::MvNormal, x::DenseVecOrMat{Float64})
for i = 1:length(x)
@inbounds x[i] = randn()
end
add!(unwhiten!(d.Σ, x), d.μ)
end
###########################################################
#
......
using Distributions
using Base.Test
m = MixtureModel([Poisson(2.0), Binomial(10, 0.3)],
[0.4, 0.6])
x = rand(m, (100,))
@test_approx_eq exp(logpdf(m,x)) pdf(m,x)
@test Distributions.variate_form(typeof(m))=== Univariate
@test Distributions.value_support(typeof(m))=== Discrete
# Core testing procedure
function test_mixture(g::UnivariateMixture, n::Int, ns::Int)
X = zeros(n)
for i = 1:10
X[i] = rand(g)
end
cs = components(g)
pr = priorprobs(g)
@assert length(cs) == length(pr)
# mean
mu = 0.0
for k = 1:length(cs)
mu += pr[k] * mean(cs[k])
end
@test_approx_eq mean(g) mu
# ground-truth
p0 = zeros(n)
for i = 1:n
p0_i = 0.0
x_i = X[i]
for k = 1:length(cs)
p0_i += pr[k] * pdf(cs[k], x_i)
end
p0[i] = p0_i
end
lp0 = log(p0)
for i = 1:n
@test_approx_eq pdf(g, X[i]) p0[i]
@test_approx_eq logpdf(g, X[i]) lp0[i]
end
p_e = pdf(g, X)
lp_e = logpdf(g, X)
@assert isa(p_e, Vector{Float64}) && length(p_e) == n
@assert isa(lp_e, Vector{Float64}) && length(lp_e) == n
@test_approx_eq p_e p0
@test_approx_eq lp_e lp0
# sampling
Xs = rand(g, ns)
@test isa(Xs, Vector{Float64})
@test length(Xs) == ns
@test_approx_eq_eps mean(Xs) mean(g) 0.01
end
function test_mixture(g::MultivariateMixture, n::Int, ns::Int)
X = zeros(length(g), n)
for i = 1:10
X[:,i] = rand(g)
end
cs = components(g)
pr = priorprobs(g)
@assert length(cs) == length(pr)
# mean
mu = 0.0
for k = 1:length(cs)
mu += pr[k] * mean(cs[k])
end
@test_approx_eq mean(g) mu
# ground-truth
p0 = zeros(n)
for i = 1:n
p0_i = 0.0
x_i = X[:,i]
for k = 1:length(cs)
p0_i += pr[k] * pdf(cs[k], x_i)
end
p0[i] = p0_i
end
lp0 = log(p0)
for i = 1:n
@test_approx_eq pdf(g, X[:,i]) p0[i]
@test_approx_eq logpdf(g, X[:,i]) lp0[i]
end
p_e = pdf(g, X)
lp_e = logpdf(g, X)
@assert isa(p_e, Vector{Float64}) && length(p_e) == n
@assert isa(lp_e, Vector{Float64}) && length(lp_e) == n
@test_approx_eq p_e p0
@test_approx_eq lp_e lp0
# sampling
Xs = rand(g, ns)
@test isa(Xs, Matrix{Float64})
@test size(Xs) == (length(g), ns)
@test_approx_eq_eps vec(mean(Xs, 2)) mean(g) 0.01
end
# Tests
println(" testing UnivariateMixture")
g_u = MixtureModel(
Normal[Normal(0.0, 1.0), Normal(2.0, 1.0), Normal(-4.0, 1.5)],
[0.2, 0.5, 0.3])
@test isa(g_u, MixtureModel{Univariate, Continuous, Normal})
test_mixture(g_u, 1000, 10^6)
println(" testing MultivariateMixture")
g_m = MixtureModel(
IsoNormal[ MvNormal([0.0, 0.0], 1.0),
MvNormal([0.2, 1.0], 1.0),
MvNormal([-0.5, -3.0], 1.6) ],
[0.2, 0.5, 0.3])
@test isa(g_m, MixtureModel{Multivariate, Continuous, IsoNormal})
@test length(g_m) == 2
test_mixture(g_m, 1000, 10^6)
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment