From 9caf4a9299a8a7e1687d08c88b952d54fbd3ba15 Mon Sep 17 00:00:00 2001 From: Kai Xu Date: Sat, 27 May 2017 15:40:06 +0100 Subject: [PATCH 01/11] Make uni gmm support generic type --- src/mixtures/unigmm.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/mixtures/unigmm.jl b/src/mixtures/unigmm.jl index 241a232aff..c2394e2583 100644 --- a/src/mixtures/unigmm.jl +++ b/src/mixtures/unigmm.jl @@ -1,12 +1,12 @@ # Univariate Gaussian Mixture Models -immutable UnivariateGMM <: UnivariateMixture{Continuous,Normal} +immutable UnivariateGMM{T<:Real} <: UnivariateMixture{Continuous,Normal} K::Int - means::Vector{Float64} - stds::Vector{Float64} + means::Vector{T} + stds::Vector{T} prior::Categorical - function UnivariateGMM(ms::Vector{Float64}, ss::Vector{Float64}, pri::Categorical) + function UnivariateGMM(ms::Vector{T}, ss::Vector{T}, pri::Categorical) K = length(ms) length(ss) == K || throw(DimensionMismatch()) ncategories(pri) == K || @@ -29,9 +29,9 @@ rand(d::UnivariateGMM) = (k = rand(d.prior); d.means[k] + randn() * d.stds[k]) params(d::UnivariateGMM) = (d.means, d.stds, d.prior) -immutable UnivariateGMMSampler <: Sampleable{Univariate,Continuous} - means::Vector{Float64} - stds::Vector{Float64} +immutable UnivariateGMMSampler{T<:Real} <: Sampleable{Univariate,Continuous} + means::Vector{T} + stds::Vector{T} psampler::AliasTable end From ddf88e67f00e2812692868ec050c6b115ad29bd8 Mon Sep 17 00:00:00 2001 From: Kai Xu Date: Sat, 27 May 2017 15:40:57 +0100 Subject: [PATCH 02/11] Add generic test --- test/mixture.jl | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/test/mixture.jl b/test/mixture.jl index 9386ef9ef1..d6eb107c8d 100644 --- a/test/mixture.jl +++ b/test/mixture.jl @@ -1,5 +1,6 @@ using Distributions using Base.Test +using ForwardDiff: Dual # Core testing procedure @@ -151,13 +152,16 @@ g_u = MixtureModel([TriangularDist(-1,2,0),TriangularDist(-.5,3,1),TriangularDis @test minimum(g_u) == -2.0 @test maximum(g_u) == 3.0 -g_u = UnivariateGMM([0.0, 2.0, -4.0], [1.0, 1.2, 1.5], Categorical([0.2, 0.5, 0.3])) -@test isa(g_u, UnivariateGMM) -@test ncomponents(g_u) == 3 -test_mixture(g_u, 1000, 10^6) -test_params(g_u) -@test minimum(g_u) == -Inf -@test maximum(g_u) == Inf +μ = [0.0, 2.0, -4.0]; σ = [1.0, 1.2, 1.5]; p = [0.2, 0.5, 0.3] +for T = [Float64, Dual] + g_u = UnivariateGMM(T[μ...], T[σ...], Categorical(T[p...])) + @test isa(g_u, UnivariateGMM) + @test ncomponents(g_u) == 3 + test_mixture(g_u, 1000, 10^6) + test_params(g_u) + @test minimum(g_u) == -Inf + @test maximum(g_u) == Inf +end println(" testing MultivariateMixture") g_m = MixtureModel( From 85a41bf4e6913f202475aa29ef19f5505a70de65 Mon Sep 17 00:00:00 2001 From: Kai Xu Date: Thu, 27 Jul 2017 21:39:22 +0100 Subject: [PATCH 03/11] Disable generic test --- test/mixture.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/mixture.jl b/test/mixture.jl index d6eb107c8d..799bc3a85b 100644 --- a/test/mixture.jl +++ b/test/mixture.jl @@ -153,7 +153,7 @@ g_u = MixtureModel([TriangularDist(-1,2,0),TriangularDist(-.5,3,1),TriangularDis @test maximum(g_u) == 3.0 μ = [0.0, 2.0, -4.0]; σ = [1.0, 1.2, 1.5]; p = [0.2, 0.5, 0.3] -for T = [Float64, Dual] +for T = [Float64]#, Dual] g_u = UnivariateGMM(T[μ...], T[σ...], Categorical(T[p...])) @test isa(g_u, UnivariateGMM) @test ncomponents(g_u) == 3 From 227efb8c0ad12072feddbaeba9464907d779b202 Mon Sep 17 00:00:00 2001 From: Kai Xu Date: Thu, 27 Jul 2017 21:39:40 +0100 Subject: [PATCH 04/11] Fix constructor --- src/mixtures/unigmm.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/mixtures/unigmm.jl b/src/mixtures/unigmm.jl index c2394e2583..aad680795f 100644 --- a/src/mixtures/unigmm.jl +++ b/src/mixtures/unigmm.jl @@ -6,15 +6,17 @@ immutable UnivariateGMM{T<:Real} <: UnivariateMixture{Continuous,Normal} stds::Vector{T} prior::Categorical - function UnivariateGMM(ms::Vector{T}, ss::Vector{T}, pri::Categorical) + (::Type{UnivariateGMM{T}}){T}(ms::Vector{T}, ss::Vector{T}, pri::Categorical) = begin K = length(ms) length(ss) == K || throw(DimensionMismatch()) ncategories(pri) == K || error("The number of categories in pri should be equal to the number of components.") - new(K, ms, ss, pri) + new{T}(K, ms, ss, pri) end end +UnivariateGMM{T<:Real}(ms::Vector{T}, ss::Vector{T}, pri::Categorical) = UnivariateGMM{T}(K, ms, ss, pri) + @distr_support UnivariateGMM -Inf Inf ncomponents(d::UnivariateGMM) = d.K From a82592629c0b49a21a93fccc826e021514fe6079 Mon Sep 17 00:00:00 2001 From: Kai Xu Date: Thu, 27 Jul 2017 21:52:22 +0100 Subject: [PATCH 05/11] Add generic test back --- test/mixture.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/mixture.jl b/test/mixture.jl index 799bc3a85b..d6eb107c8d 100644 --- a/test/mixture.jl +++ b/test/mixture.jl @@ -153,7 +153,7 @@ g_u = MixtureModel([TriangularDist(-1,2,0),TriangularDist(-.5,3,1),TriangularDis @test maximum(g_u) == 3.0 μ = [0.0, 2.0, -4.0]; σ = [1.0, 1.2, 1.5]; p = [0.2, 0.5, 0.3] -for T = [Float64]#, Dual] +for T = [Float64, Dual] g_u = UnivariateGMM(T[μ...], T[σ...], Categorical(T[p...])) @test isa(g_u, UnivariateGMM) @test ncomponents(g_u) == 3 From 0e4f5af0c26862d2394825d8ac5f9ea18abbfe05 Mon Sep 17 00:00:00 2001 From: Kai Xu Date: Fri, 28 Jul 2017 16:20:45 +0100 Subject: [PATCH 06/11] Fix K and try to make mixture.jl generic --- src/mixtures/mixturemodel.jl | 18 +++++----- src/mixtures/unigmm.jl | 2 +- test/mixture.jl | 66 ++++++++++++++++++++++++++++++++++++ 3 files changed, 76 insertions(+), 10 deletions(-) diff --git a/src/mixtures/mixturemodel.jl b/src/mixtures/mixturemodel.jl index 2c2ca75e2b..a50dc1503c 100644 --- a/src/mixtures/mixturemodel.jl +++ b/src/mixtures/mixturemodel.jl @@ -226,7 +226,7 @@ function _mixpdf!(r::AbstractArray, d::AbstractMixtureModel, x) return r end -function _mixlogpdf1(d::AbstractMixtureModel, x) +function _mixlogpdf1{T<:Real}(d::AbstractMixtureModel, x::T) # using the formula below for numerical stability # # logpdf(d, x) = log(sum_i pri[i] * pdf(cs[i], x)) @@ -242,11 +242,11 @@ function _mixlogpdf1(d::AbstractMixtureModel, x) p = probs(d) @assert length(p) == K - lp = Vector{Float64}(K) + lp = Vector{T}(K) m = -Inf # m <- the maximum of log(p(cs[i], x)) + log(pri[i]) for i = 1:K @inbounds pi = p[i] - if pi > 0.0 + if pi > one(T) # lp[i] <- log(p(cs[i], x)) + log(pri[i]) lp_i = logpdf(component(d, i), x) + log(pi) @inbounds lp[i] = lp_i @@ -255,25 +255,25 @@ function _mixlogpdf1(d::AbstractMixtureModel, x) end end end - v = 0.0 + v = one(T) @inbounds for i = 1:K - if p[i] > 0.0 + if p[i] > one(T) v += exp(lp[i] - m) end end return m + log(v) end -function _mixlogpdf!(r::AbstractArray, d::AbstractMixtureModel, x) +function _mixlogpdf!{T<:Real}(r::AbstractArray, d::AbstractMixtureModel, x::AbstractArray{T}) K = ncomponents(d) p = probs(d) @assert length(p) == K n = length(r) - Lp = Matrix{Float64}(n, K) + Lp = Matrix{T}(n, K) m = fill(-Inf, n) for i = 1:K @inbounds pi = p[i] - if pi > 0.0 + if pi > one(T) lpri = log(pi) lp_i = view(Lp, :, i) # compute logpdf in batch and store @@ -290,7 +290,7 @@ function _mixlogpdf!(r::AbstractArray, d::AbstractMixtureModel, x) end end - fill!(r, 0.0) + fill!(r, one(T)) @inbounds for i = 1:K if p[i] > 0.0 lp_i = view(Lp, :, i) diff --git a/src/mixtures/unigmm.jl b/src/mixtures/unigmm.jl index aad680795f..17b40b96af 100644 --- a/src/mixtures/unigmm.jl +++ b/src/mixtures/unigmm.jl @@ -15,7 +15,7 @@ immutable UnivariateGMM{T<:Real} <: UnivariateMixture{Continuous,Normal} end end -UnivariateGMM{T<:Real}(ms::Vector{T}, ss::Vector{T}, pri::Categorical) = UnivariateGMM{T}(K, ms, ss, pri) +UnivariateGMM{T<:Real}(ms::Vector{T}, ss::Vector{T}, pri::Categorical) = UnivariateGMM{T}(ms, ss, pri) @distr_support UnivariateGMM -Inf Inf diff --git a/test/mixture.jl b/test/mixture.jl index d6eb107c8d..151198fa32 100644 --- a/test/mixture.jl +++ b/test/mixture.jl @@ -1,5 +1,6 @@ using Distributions using Base.Test +using JSON, ForwardDiff, Calculus, PDMats, Compat # test dependencies using ForwardDiff: Dual # Core testing procedure @@ -69,6 +70,71 @@ function test_mixture(g::UnivariateMixture, n::Int, ns::Int) @test isapprox(mean(Xs), mean(g), atol=0.01) end +function test_mixture(g::UnivariateGMM{Dual}, n::Int, ns::Int) + X = zeros(Dual, n) + for i = 1:n + X[i] = rand(g) + end + + K = ncomponents(g) + pr = probs(g) + @assert length(pr) == K + + # mean + mu = 0.0 + for k = 1:K + mu += pr[k] * mean(component(g, k)) + end + @test mean(g) ≈ mu + + # # evaluation of cdf + # cf = zeros(Dual, n) + # for k = 1:K + # c_k = component(g, k) + # for i = 1:n + # cf[i] += pr[k] * cdf(c_k, X[i]) + # end + # end + # + # for i = 1:n + # @test cdf(g, X[i]) ≈ cf[i] + # end + # @test cdf(g, X) ≈ cf + + # evaluation + P0 = zeros(Dual, n, K) + LP0 = zeros(Dual, n, K) + for k = 1:K + c_k = component(g, k) + for i = 1:n + x_i = X[i] + P0[i,k] = pdf(c_k, x_i) + LP0[i,k] = logpdf(c_k, x_i) + end + end + + mix_p0 = P0 * pr + mix_lp0 = @compat(log.(mix_p0)) + + for i = 1:n + @test pdf(g, X[i]) ≈ mix_p0[i] + @test logpdf(g, X[i]) ≈ mix_lp0[i] + @test componentwise_pdf(g, X[i]) ≈ vec(P0[i,:]) + @test componentwise_logpdf(g, X[i]) ≈ vec(LP0[i,:]) + end + + @test pdf(g, X) ≈ mix_p0 + @test logpdf(g, X) ≈ mix_lp0 + @test componentwise_pdf(g, X) ≈ P0 + @test componentwise_logpdf(g, X) ≈ LP0 + + # sampling + Xs = rand(g, ns) + @test isa(Xs, Vector{Float64}) + @test length(Xs) == ns + @test isapprox(mean(Xs), mean(g), atol=0.01) +end + function test_mixture(g::MultivariateMixture, n::Int, ns::Int) X = zeros(length(g), n) for i = 1:n From ebd7164b1eacbf4fa7481de69ccd93d7b463e0b5 Mon Sep 17 00:00:00 2001 From: Kai Xu Date: Thu, 28 Feb 2019 19:29:51 +0000 Subject: [PATCH 07/11] Support Real --- src/mixtures/mixturemodel.jl | 38 +++++----- src/mixtures/unigmm.jl | 14 ++-- test/mixture.jl | 136 ++++++++++++++++++----------------- 3 files changed, 98 insertions(+), 90 deletions(-) diff --git a/src/mixtures/mixturemodel.jl b/src/mixtures/mixturemodel.jl index 936930d34e..94a311310b 100644 --- a/src/mixtures/mixturemodel.jl +++ b/src/mixtures/mixturemodel.jl @@ -135,13 +135,13 @@ function MixtureModel(components::Vector{C}, prior::Categorical) where C<:Distri MixtureModel{VF,VS,C}(components, prior) end -MixtureModel(components::Vector{C}, p::Vector{Float64}) where {C<:Distribution} = +MixtureModel(components::Vector{C}, p::Vector{T}) where {C<:Distribution,T<:Real} = MixtureModel(components, Categorical(p)) _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{Float64}) where C<:Distribution +function MixtureModel(::Type{C}, params::AbstractArray, p::Vector{T}) where {C<:Distribution,T<:Real} components = C[_construct_component(C, a) for a in params] MixtureModel(components, p) end @@ -343,7 +343,7 @@ function _mixpdf!(r::AbstractArray, d::AbstractMixtureModel, x) return r end -function _mixlogpdf1{T<:Real}(d::AbstractMixtureModel, x::T) +function _mixlogpdf1(d::AbstractMixtureModel, x) # using the formula below for numerical stability # # logpdf(d, x) = log(sum_i pri[i] * pdf(cs[i], x)) @@ -359,11 +359,11 @@ function _mixlogpdf1{T<:Real}(d::AbstractMixtureModel, x::T) p = probs(d) @assert length(p) == K - lp = Vector{Float64}(undef, K) + lp = Vector{eltype(x)}(undef, K) m = -Inf # m <- the maximum of log(p(cs[i], x)) + log(pri[i]) for i = 1:K @inbounds pi = p[i] - if pi > one(T) + if pi > 0.0 # lp[i] <- log(p(cs[i], x)) + log(pri[i]) lp_i = logpdf(component(d, i), x) + log(pi) @inbounds lp[i] = lp_i @@ -372,25 +372,25 @@ function _mixlogpdf1{T<:Real}(d::AbstractMixtureModel, x::T) end end end - v = one(T) + v = 0.0 @inbounds for i = 1:K - if p[i] > one(T) + if p[i] > 0.0 v += exp(lp[i] - m) end end return m + log(v) end -function _mixlogpdf!{T<:Real}(r::AbstractArray, d::AbstractMixtureModel, x::AbstractArray{T}) +function _mixlogpdf!(r::AbstractArray, d::AbstractMixtureModel, x) K = ncomponents(d) p = probs(d) @assert length(p) == K n = length(r) - Lp = Matrix{Float64}(undef, n, K) + Lp = Matrix{eltype(x)}(undef, n, K) m = fill(-Inf, n) for i = 1:K @inbounds pi = p[i] - if pi > one(T) + if pi > 0.0 lpri = log(pi) lp_i = view(Lp, :, i) # compute logpdf in batch and store @@ -412,7 +412,7 @@ function _mixlogpdf!{T<:Real}(r::AbstractArray, d::AbstractMixtureModel, x::Abst end end - fill!(r, one(T)) + fill!(r, 0.0) @inbounds for i = 1:K if p[i] > 0.0 lp_i = view(Lp, :, i) @@ -501,15 +501,15 @@ componentwise_logpdf!(r::AbstractVector, d::MultivariateMixture, x::AbstractVect componentwise_logpdf!(r::AbstractMatrix, d::UnivariateMixture, x::AbstractVector) = _cwise_logpdf!(r, d, x) componentwise_logpdf!(r::AbstractMatrix, d::MultivariateMixture, x::AbstractMatrix) = _cwise_logpdf!(r, d, x) -componentwise_pdf(d::UnivariateMixture, x::Real) = componentwise_pdf!(Vector{Float64}(undef, ncomponents(d)), d, x) -componentwise_pdf(d::UnivariateMixture, x::AbstractVector) = componentwise_pdf!(Matrix{Float64}(undef, length(x), ncomponents(d)), d, x) -componentwise_pdf(d::MultivariateMixture, x::AbstractVector) = componentwise_pdf!(Vector{Float64}(undef, ncomponents(d)), d, x) -componentwise_pdf(d::MultivariateMixture, x::AbstractMatrix) = componentwise_pdf!(Matrix{Float64}(undef, size(x,2), ncomponents(d)), d, x) +componentwise_pdf(d::UnivariateMixture, x::Real) = componentwise_pdf!(Vector{eltype(x)}(undef, ncomponents(d)), d, x) +componentwise_pdf(d::UnivariateMixture, x::AbstractVector) = componentwise_pdf!(Matrix{eltype(x)}(undef, length(x), ncomponents(d)), d, x) +componentwise_pdf(d::MultivariateMixture, x::AbstractVector) = componentwise_pdf!(Vector{eltype(x)}(undef, ncomponents(d)), d, x) +componentwise_pdf(d::MultivariateMixture, x::AbstractMatrix) = componentwise_pdf!(Matrix{eltype(x)}(undef, size(x,2), ncomponents(d)), d, x) -componentwise_logpdf(d::UnivariateMixture, x::Real) = componentwise_logpdf!(Vector{Float64}(undef, ncomponents(d)), d, x) -componentwise_logpdf(d::UnivariateMixture, x::AbstractVector) = componentwise_logpdf!(Matrix{Float64}(undef, length(x), ncomponents(d)), d, x) -componentwise_logpdf(d::MultivariateMixture, x::AbstractVector) = componentwise_logpdf!(Vector{Float64}(undef, ncomponents(d)), d, x) -componentwise_logpdf(d::MultivariateMixture, x::AbstractMatrix) = componentwise_logpdf!(Matrix{Float64}(undef, size(x,2), ncomponents(d)), d, x) +componentwise_logpdf(d::UnivariateMixture, x::Real) = componentwise_logpdf!(Vector{eltype(x)}(undef, ncomponents(d)), d, x) +componentwise_logpdf(d::UnivariateMixture, x::AbstractVector) = componentwise_logpdf!(Matrix{eltype(x)}(undef, length(x), ncomponents(d)), d, x) +componentwise_logpdf(d::MultivariateMixture, x::AbstractVector) = componentwise_logpdf!(Vector{eltype(x)}(undef, ncomponents(d)), d, x) +componentwise_logpdf(d::MultivariateMixture, x::AbstractMatrix) = componentwise_logpdf!(Matrix{eltype(x)}(undef, size(x,2), ncomponents(d)), d, x) ## Sampling diff --git a/src/mixtures/unigmm.jl b/src/mixtures/unigmm.jl index 2eaee69540..71fdd1a6ca 100644 --- a/src/mixtures/unigmm.jl +++ b/src/mixtures/unigmm.jl @@ -1,12 +1,12 @@ # Univariate Gaussian Mixture Models -struct UnivariateGMM <: UnivariateMixture{Continuous,Normal} +struct UnivariateGMM{T<:Real} <: UnivariateMixture{Continuous,Normal} K::Int means::Vector{T} stds::Vector{T} prior::Categorical - (::Type{UnivariateGMM{T}}){T}(ms::Vector{T}, ss::Vector{T}, pri::Categorical) = begin + function UnivariateGMM{T}(ms::Vector{T}, ss::Vector{T}, pri::Categorical) where {T<:Real} K = length(ms) length(ss) == K || throw(DimensionMismatch()) ncategories(pri) == K || @@ -15,7 +15,7 @@ struct UnivariateGMM <: UnivariateMixture{Continuous,Normal} end end -UnivariateGMM{T<:Real}(ms::Vector{T}, ss::Vector{T}, pri::Categorical) = UnivariateGMM{T}(ms, ss, pri) +UnivariateGMM(ms::Vector{T}, ss::Vector{T}, pri::Categorical) where {T<:Real} = UnivariateGMM{T}(ms, ss, pri) @distr_support UnivariateGMM -Inf Inf @@ -27,14 +27,16 @@ probs(d::UnivariateGMM) = probs(d.prior) mean(d::UnivariateGMM) = dot(d.means, probs(d)) +rand(d::UnivariateGMM) = (k = rand(d.prior); d.means[k] + randn() * d.stds[k]) + rand(rng::AbstractRNG, d::UnivariateGMM) = (k = rand(rng, d.prior); d.means[k] + randn(rng) * d.stds[k]) params(d::UnivariateGMM) = (d.means, d.stds, d.prior) -struct UnivariateGMMSampler <: Sampleable{Univariate,Continuous} - means::Vector{Float64} - stds::Vector{Float64} +struct UnivariateGMMSampler{T<:Real} <: Sampleable{Univariate,Continuous} + means::Vector{T} + stds::Vector{T} psampler::AliasTable end diff --git a/test/mixture.jl b/test/mixture.jl index 003cd13ae0..d79f34403c 100644 --- a/test/mixture.jl +++ b/test/mixture.jl @@ -79,8 +79,9 @@ function test_mixture(g::UnivariateMixture, n::Int, ns::Int, @test isapprox(mean(Xs), mean(g), atol=0.01) end -function test_mixture(g::UnivariateGMM{Dual}, n::Int, ns::Int) - X = zeros(Dual, n) +function test_mixture(g::UnivariateGMM{T}, n::Int, ns::Int, + rng::Union{AbstractRNG, Missing} = missing) where {T<:Real} + X = zeros(T, n) for i = 1:n X[i] = rand(g) end @@ -97,7 +98,7 @@ function test_mixture(g::UnivariateGMM{Dual}, n::Int, ns::Int) @test mean(g) ≈ mu # # evaluation of cdf - # cf = zeros(Dual, n) + # cf = zeros(T, n) # for k = 1:K # c_k = component(g, k) # for i = 1:n @@ -111,8 +112,8 @@ function test_mixture(g::UnivariateGMM{Dual}, n::Int, ns::Int) # @test cdf(g, X) ≈ cf # evaluation - P0 = zeros(Dual, n, K) - LP0 = zeros(Dual, n, K) + P0 = zeros(T, n, K) + LP0 = zeros(T, n, K) for k = 1:K c_k = component(g, k) for i = 1:n @@ -123,7 +124,7 @@ function test_mixture(g::UnivariateGMM{Dual}, n::Int, ns::Int) end mix_p0 = P0 * pr - mix_lp0 = @compat(log.(mix_p0)) + mix_lp0 = log.(mix_p0) for i = 1:n @test pdf(g, X[i]) ≈ mix_p0[i] @@ -132,16 +133,20 @@ function test_mixture(g::UnivariateGMM{Dual}, n::Int, ns::Int) @test componentwise_logpdf(g, X[i]) ≈ vec(LP0[i,:]) end - @test pdf(g, X) ≈ mix_p0 - @test logpdf(g, X) ≈ mix_lp0 + @test pdf.(g, X) ≈ mix_p0 + @test logpdf.(g, X) ≈ mix_lp0 @test componentwise_pdf(g, X) ≈ P0 @test componentwise_logpdf(g, X) ≈ LP0 # sampling - Xs = rand(g, ns) - @test isa(Xs, Vector{Float64}) - @test length(Xs) == ns - @test isapprox(mean(Xs), mean(g), atol=0.01) + # if ismissing(rng) + # Xs = rand(g, ns) + # else + # Xs = rand(rng, g, ns) + # end + # @test isa(Xs, Vector{T}) + # @test length(Xs) == ns + # @test isapprox(mean(Xs), mean(g), atol=0.01) end function test_mixture(g::MultivariateMixture, n::Int, ns::Int, @@ -230,58 +235,59 @@ end Dict("rand(...)" => missing, "rand(rng, ...)" => MersenneTwister(123)) -@testset "Testing UnivariateMixture" begin -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 ncomponents(g_u) == 3 -test_mixture(g_u, 1000, 10^6, rng) -test_params(g_u) -@test minimum(g_u) == -Inf -@test maximum(g_u) == Inf -@test extrema(g_u) == (-Inf, Inf) - -g_u = MixtureModel([TriangularDist(-1,2,0),TriangularDist(-.5,3,1),TriangularDist(-2,0,-1)]) -@test minimum(g_u) ≈ -2.0 -@test maximum(g_u) ≈ 3.0 -@test extrema(g_u) == (minimum(g_u), maximum(g_u)) -@test insupport(g_u, 2.5) == true -@test insupport(g_u, 3.5) == false - -μ = [0.0, 2.0, -4.0]; σ = [1.0, 1.2, 1.5]; p = [0.2, 0.5, 0.3] -for T = [Float64, Dual] - g_u = UnivariateGMM(T[μ...], T[σ...], Categorical(T[p...])) - @test isa(g_u, UnivariateGMM) - @test ncomponents(g_u) == 3 - test_mixture(g_u, 1000, 10^6, rng) - test_params(g_u) - @test minimum(g_u) == -Inf - @test maximum(g_u) == Inf - @test extrema(g_u) == (-Inf, Inf) -end + @testset "Testing UnivariateMixture" begin + 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 ncomponents(g_u) == 3 + test_mixture(g_u, 1000, 10^6, rng) + test_params(g_u) + @test minimum(g_u) == -Inf + @test maximum(g_u) == Inf + @test extrema(g_u) == (-Inf, Inf) + + g_u = MixtureModel([TriangularDist(-1,2,0),TriangularDist(-.5,3,1),TriangularDist(-2,0,-1)]) + @test minimum(g_u) ≈ -2.0 + @test maximum(g_u) ≈ 3.0 + @test extrema(g_u) == (minimum(g_u), maximum(g_u)) + @test insupport(g_u, 2.5) == true + @test insupport(g_u, 3.5) == false + + μ = [0.0, 2.0, -4.0]; σ = [1.0, 1.2, 1.5]; p = [0.2, 0.5, 0.3] + for T = [Float64, Dual] + g_u = UnivariateGMM(map(Dual, μ), map(Dual, σ), Categorical(p)) + @test isa(g_u, UnivariateGMM) + @test ncomponents(g_u) == 3 + test_mixture(g_u, 1000, 10^6, rng) + test_params(g_u) + @test minimum(g_u) == -Inf + @test maximum(g_u) == Inf + @test extrema(g_u) == (-Inf, Inf) + end + end -@testset "Testing MultivariatevariateMixture" begin -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(components(g_m)) == 3 -@test length(g_m) == 2 -@test insupport(g_m, [0.0, 0.0]) == true -test_mixture(g_m, 1000, 10^6, rng) -test_params(g_m) - -u1 = Uniform() -u2 = Uniform(1.0, 2.0) -utot =Uniform(0.0, 2.0) - -# mixture supposed to be a uniform on [0.0,2.0] -unif_mixt = MixtureModel([u1,u2]) -@test var(utot) ≈ var(unif_mixt) -@test mean(utot) ≈ mean(unif_mixt) -for x in -1.0:0.5:2.5 - @test cdf(utot,x) ≈ cdf(utot,x) -end -end + @testset "Testing MultivariatevariateMixture" begin + 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(components(g_m)) == 3 + @test length(g_m) == 2 + @test insupport(g_m, [0.0, 0.0]) == true + test_mixture(g_m, 1000, 10^6, rng) + test_params(g_m) + + u1 = Uniform() + u2 = Uniform(1.0, 2.0) + utot =Uniform(0.0, 2.0) + + # mixture supposed to be a uniform on [0.0,2.0] + unif_mixt = MixtureModel([u1,u2]) + @test var(utot) ≈ var(unif_mixt) + @test mean(utot) ≈ mean(unif_mixt) + for x in -1.0:0.5:2.5 + @test cdf(utot,x) ≈ cdf(utot,x) + end + end end From baf67d56a327424bb93a9c1b9e25a6c20a83b297 Mon Sep 17 00:00:00 2001 From: Kai Xu Date: Thu, 28 Feb 2019 19:36:45 +0000 Subject: [PATCH 08/11] enable cdf tests --- src/mixtures/mixturemodel.jl | 2 +- test/mixture.jl | 26 +++++++++++++------------- 2 files changed, 14 insertions(+), 14 deletions(-) diff --git a/src/mixtures/mixturemodel.jl b/src/mixtures/mixturemodel.jl index 94a311310b..5ef0cf89aa 100644 --- a/src/mixtures/mixturemodel.jl +++ b/src/mixtures/mixturemodel.jl @@ -304,7 +304,7 @@ function _cdf(d::UnivariateMixture, x::Real) return r end -cdf(d::UnivariateMixture{Continuous}, x::Float64) = _cdf(d, x) +cdf(d::UnivariateMixture{Continuous}, x::Real) = _cdf(d, x) cdf(d::UnivariateMixture{Discrete}, x::Int) = _cdf(d, x) diff --git a/test/mixture.jl b/test/mixture.jl index d79f34403c..689a533895 100644 --- a/test/mixture.jl +++ b/test/mixture.jl @@ -97,19 +97,19 @@ function test_mixture(g::UnivariateGMM{T}, n::Int, ns::Int, end @test mean(g) ≈ mu - # # evaluation of cdf - # cf = zeros(T, n) - # for k = 1:K - # c_k = component(g, k) - # for i = 1:n - # cf[i] += pr[k] * cdf(c_k, X[i]) - # end - # end - # - # for i = 1:n - # @test cdf(g, X[i]) ≈ cf[i] - # end - # @test cdf(g, X) ≈ cf + # evaluation of cdf + cf = zeros(T, n) + for k = 1:K + c_k = component(g, k) + for i = 1:n + cf[i] += pr[k] * cdf(c_k, X[i]) + end + end + + for i = 1:n + @test cdf(g, X[i]) ≈ cf[i] + end + @test cdf.(g, X) ≈ cf # evaluation P0 = zeros(T, n, K) From 90d090581c6b363e651a8399f8998ff467fd7f47 Mon Sep 17 00:00:00 2001 From: Kai Xu Date: Thu, 28 Feb 2019 19:44:45 +0000 Subject: [PATCH 09/11] merge tests --- test/mixture.jl | 93 +++++++------------------------------------------ 1 file changed, 13 insertions(+), 80 deletions(-) diff --git a/test/mixture.jl b/test/mixture.jl index 689a533895..dbc8c46a96 100644 --- a/test/mixture.jl +++ b/test/mixture.jl @@ -7,80 +7,11 @@ using ForwardDiff: Dual function test_mixture(g::UnivariateMixture, n::Int, ns::Int, rng::Union{AbstractRNG, Missing} = missing) - X = zeros(n) - for i = 1:n - if ismissing(rng) - X[i] = rand(g) - else - X[i] = rand(rng, g) - end - end - - K = ncomponents(g) - pr = probs(g) - @assert length(pr) == K - - # mean - mu = 0.0 - for k = 1:K - mu += pr[k] * mean(component(g, k)) - end - @test mean(g) ≈ mu - - # evaluation of cdf - cf = zeros(n) - for k = 1:K - c_k = component(g, k) - for i = 1:n - cf[i] += pr[k] * cdf(c_k, X[i]) - end - end - - for i = 1:n - @test cdf(g, X[i]) ≈ cf[i] - end - @test cdf.(g, X) ≈ cf - - # evaluation - P0 = zeros(n, K) - LP0 = zeros(n, K) - for k = 1:K - c_k = component(g, k) - for i = 1:n - x_i = X[i] - P0[i,k] = pdf(c_k, x_i) - LP0[i,k] = logpdf(c_k, x_i) - end - end - - mix_p0 = P0 * pr - mix_lp0 = log.(mix_p0) - - for i = 1:n - @test pdf(g, X[i]) ≈ mix_p0[i] - @test logpdf(g, X[i]) ≈ mix_lp0[i] - @test componentwise_pdf(g, X[i]) ≈ vec(P0[i,:]) - @test componentwise_logpdf(g, X[i]) ≈ vec(LP0[i,:]) - end - - @test pdf.(g, X) ≈ mix_p0 - @test logpdf.(g, X) ≈ mix_lp0 - @test componentwise_pdf(g, X) ≈ P0 - @test componentwise_logpdf(g, X) ≈ LP0 - - # sampling - if ismissing(rng) - Xs = rand(g, ns) + if g isa UnivariateGMM + T = eltype(g.means) else - Xs = rand(rng, g, ns) + T = eltype(g) end - @test isa(Xs, Vector{Float64}) - @test length(Xs) == ns - @test isapprox(mean(Xs), mean(g), atol=0.01) -end - -function test_mixture(g::UnivariateGMM{T}, n::Int, ns::Int, - rng::Union{AbstractRNG, Missing} = missing) where {T<:Real} X = zeros(T, n) for i = 1:n X[i] = rand(g) @@ -139,14 +70,16 @@ function test_mixture(g::UnivariateGMM{T}, n::Int, ns::Int, @test componentwise_logpdf(g, X) ≈ LP0 # sampling - # if ismissing(rng) - # Xs = rand(g, ns) - # else - # Xs = rand(rng, g, ns) - # end - # @test isa(Xs, Vector{T}) - # @test length(Xs) == ns - # @test isapprox(mean(Xs), mean(g), atol=0.01) + if (T isa AbstractFloat) + if ismissing(rng) + Xs = rand(g, ns) + else + Xs = rand(rng, g, ns) + end + @test isa(Xs, Vector{T}) + @test length(Xs) == ns + @test isapprox(mean(Xs), mean(g), atol=0.01) + end end function test_mixture(g::MultivariateMixture, n::Int, ns::Int, From 7f41b0f5ec7233f383f1c2a36c1231845747116f Mon Sep 17 00:00:00 2001 From: Kai Xu Date: Thu, 14 Mar 2019 21:47:11 +0000 Subject: [PATCH 10/11] allow mean and std to have differetn types --- src/mixtures/unigmm.jl | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/src/mixtures/unigmm.jl b/src/mixtures/unigmm.jl index 71fdd1a6ca..f869031940 100644 --- a/src/mixtures/unigmm.jl +++ b/src/mixtures/unigmm.jl @@ -1,22 +1,20 @@ # Univariate Gaussian Mixture Models -struct UnivariateGMM{T<:Real} <: UnivariateMixture{Continuous,Normal} +struct UnivariateGMM{VT1<:AbstractVector{<:Real},VT2<:AbstractVector{<:Real}} <: UnivariateMixture{Continuous,Normal} K::Int - means::Vector{T} - stds::Vector{T} + means::VT1 + stds::VT2 prior::Categorical - function UnivariateGMM{T}(ms::Vector{T}, ss::Vector{T}, pri::Categorical) where {T<:Real} + function UnivariateGMM(ms::VT1, ss::VT2, pri::Categorical) where {VT1<:AbstractVector{<:Real},VT2<:AbstractVector{<:Real}} K = length(ms) length(ss) == K || throw(DimensionMismatch()) ncategories(pri) == K || error("The number of categories in pri should be equal to the number of components.") - new{T}(K, ms, ss, pri) + new{VT1,VT2}(K, ms, ss, pri) end end -UnivariateGMM(ms::Vector{T}, ss::Vector{T}, pri::Categorical) where {T<:Real} = UnivariateGMM{T}(ms, ss, pri) - @distr_support UnivariateGMM -Inf Inf ncomponents(d::UnivariateGMM) = d.K @@ -34,9 +32,9 @@ rand(rng::AbstractRNG, d::UnivariateGMM) = params(d::UnivariateGMM) = (d.means, d.stds, d.prior) -struct UnivariateGMMSampler{T<:Real} <: Sampleable{Univariate,Continuous} - means::Vector{T} - stds::Vector{T} +struct UnivariateGMMSampler{VT1<:AbstractVector{<:Real},VT2<:AbstractVector{<:Real}} <: Sampleable{Univariate,Continuous} + means::VT1 + stds::VT2 psampler::AliasTable end From 875b4e517832475427b3c841365f3b23edf659e9 Mon Sep 17 00:00:00 2001 From: Kai Xu Date: Thu, 14 Mar 2019 21:51:56 +0000 Subject: [PATCH 11/11] allow p to be abstract vector --- src/mixtures/mixturemodel.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mixtures/mixturemodel.jl b/src/mixtures/mixturemodel.jl index 5ef0cf89aa..d2eaba9f0c 100644 --- a/src/mixtures/mixturemodel.jl +++ b/src/mixtures/mixturemodel.jl @@ -135,7 +135,7 @@ function MixtureModel(components::Vector{C}, prior::Categorical) where C<:Distri MixtureModel{VF,VS,C}(components, prior) end -MixtureModel(components::Vector{C}, p::Vector{T}) where {C<:Distribution,T<:Real} = +MixtureModel(components::Vector{C}, p::VT) where {C<:Distribution,VT<:AbstractVector{<:Real}} = MixtureModel(components, Categorical(p)) _construct_component(::Type{C}, arg) where {C<:Distribution} = C(arg)