Multivariate dirichlet mixture with Turing

Specific Domains Probabilistic programming
1

Hi – I’m trying to fit a mixture of 4 bivariate gaussians with Turing and getting a strange error? Can someone help with the syntax here:

@model GaussianMixtureModel(x, K) = begin
# K is expected number of clusters
# x is bivariate data, typically 2 x 1000

D, N = size(x)

# Draw the paramters
μ ~ [MvNormal([0.,0.], 5.) for i in 1:K]
τ ~ [Gamma() for i in 1:K]

# Dirichlet prior
α = 1.0
w ~ Dirichlet(K, α)

# Draw assignments for each datum and generate it from a multivariate normal.
k = Vector{Int}(undef, N)
for i in 1:N
    k[i] ~ Categorical(w)
    x[:, i] ~ MvNormal(μ[k[i]], τ[k[i]])
end
    
return k

end

I’m getting error (when sampling): MethodError: no method matching assume(::DynamicPPL.SampleFromPrior,

Thanks

2

This should work.

@model GaussianMixtureModel(x, K) = begin
    # K is expected number of clusters
    # x is bivariate data, typically 2 x 1000

    D, N = size(x)

    # Draw the paramters
    μ ~ filldist(MvNormal([0.,0.], 5.), K)
    τ ~ filldist(Gamma(), K)

    # Dirichlet prior
    α = 1.0
    w ~ Dirichlet(K, α)

    # Draw assignments for each datum and generate it from a multivariate normal.
    k = Vector{Int}(undef, N)
    for i in 1:N
        k[i] ~ Categorical(w)
        x[:, i] ~ MvNormal(μ[:, k[i]], τ[k[i]])
    end

    return k
end
3

Awesome - this works! Thanks

1 Like
4

Are there any issues sampling from MvNormal() as used in the model above?

gmm_model = GaussianMixtureModel(x, 8);
num_iter = 1000
samp_mh = sample(gmm_model, MH(), num_iter);
samp_hmc = sample(gmm_model, HMC(0.1, 5), num_iter);
samp_nuts = sample(gmm_model, NUTS(100, 0.65), num_iter);

If I try these three samplers MH() runs ok and completes. However both NUTS() and HMC() complain with the same error message: “TypeError: in typeassert, expected Int64, got ForwardDiff.Dual{Nothing,Int64,10}”

Any ideas? The synthetic dataset I’m trying to run this on is generated as follows:

y1 = rand(MvNormal([0., 0.], [[1, 0.] [0., 0.02]]), 250);
y2 = rand(MvNormal([0., 0.], [[0.02, 0.] [0., 1]]), 400);
y3 = rand(MvNormal([2., 2.], [[1., -0.9] [-0.9, 1]]), 500);
y4 = rand(MvNormal([-2., -2.], [[0.1, 0.] [0., 0.1]]), 300);
x = reduce(hcat, [y1, y2, y3, y4])

Thanks

5

You can’t use pure HMC when you have integer random variables, k above. You need Gibbs for that.

6

Does filldist() work as well for InverseWishart()? I’m having some issues setting up covariance prior as:

τ ~ filldist(InverseWishart(2+1, Matrix{Float64}([[3,0.5] [0.5,3]])), K)

and getting error message: “MethodError: no method matching filldist(::InverseWishart{Float64,PDMats.PDMat{Float64,Array{Float64,2}}}, ::Int64)”

Here’s the full model:

@model GaussianMixtureModel(x, K) = begin
    # K is expected number of clusters
    # x is bivariate data, typically 2 x 1000
    
    D, N = size(x)
    
     # Draw the paramters
    μ ~ filldist(MvNormal([0.,0.], 5.), K)
    τ ~ filldist(InverseWishart(2+1, Matrix{Float64}([[3,0.5] [0.5,3]])), K)
    
    # Dirichlet prior
    α = 1.0
    w ~ Dirichlet(K, α)
    
    # Draw assignments for each datum and generate it from a multivariate normal.
    k = Vector{Int}(undef, N)
    for i in 1:N
        k[i] ~ Categorical(w)
        x[:, i] ~ MvNormal(μ[:, k[i]], τ[:, :, k[i]])
    end
        
    return k
end

Thanks

7

Please open an issue in DistributionsAD.jl.