# Using MixtureModels (from Distributions.jl) in Turing? Issues with posterior prediction

**URL:** <https://discourse.julialang.org/t/using-mixturemodels-from-distributions-jl-in-turing-issues-with-posterior-prediction/71617>\
**Category:** Probabilistic Programming\
**Tags:** turing\
**Created:** [November 16, 2021, 10:13pm UTC](https://discourse.julialang.org/t/using-mixturemodels-from-distributions-jl-in-turing-issues-with-posterior-prediction/71617 "2021-11-16T22:13:03Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![opera\_malenky](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/opera_malenky/32/8213_2.png) [@opera\_malenky](https://discourse.julialang.org/u/opera_malenky)\
**Post date:** [November 16, 2021, 10:13pm UTC](https://discourse.julialang.org/t/using-mixturemodels-from-distributions-jl-in-turing-issues-with-posterior-prediction/71617/1 "2021-11-16T22:13:03Z")

</div>

I am trying to fit a kernel mixture model in Turing. It seems to work for estimation, but gives an error when trying to use `predict`. I’d be curious if anyone has a suggestion for a fix.

```julia
using Distributions, Turing

# Unidimensional Kernel Mixture model with K pre-specified components
# that cover the space from min_x to max_x
@model function KMM(x, min_x, max_x, k, σ)

    N = size(x, 1)
    linspan = range(min_x, stop=max_x, length=k)
    kernels = map(u -> Normal(u, σ), linspan)

    ω ~ Dirichlet(k, 1.0)
    mixdist = MixtureModel(kernels, ω)

    x ~ filldist(mixdist, N)

end

# Simulate data from a bimodal distribution
data = vcat(rand(Normal(-1, 0.5), 50), rand(Normal(1, 0.5), 50))

# Define a kernel mixture with 10 gaussian components, with means covering -2:2
model = KMM(data, -2.0, 2.0, 10, 0.5)

# Estimate weights
m1 = sample(model, NUTS(0.65), 1000)

```

That seems to work (although it is quite slow for a very small model). But when I want to get the posterior predictive distribution of the original data, it gives an error, saying that the method `loglikelihood` does not exist for the filldist of mixtures.

```julia
pp_data = predict(KMM(Vector{Union{Missing, Float64}}(missing, length(data)), -2.0, 2.0, 10, 0.5), m1)

MethodError: no method matching loglikelihood(::Product{Continuous, MixtureModel{Univariate, Continuous, Normal{Float64}, Categorical{Float64, Vector{Float64}}}, FillArrays.Fill{MixtureModel{Univariate, Continuous, Normal{Float64}, Categorical{Float64, Vector{Float64}}}, 1, Tuple{Base.OneTo{Int64}}}}, ::Vector{Union{Missing, Float64}})

```

`MixtureModel` from `Distributions.jl` apparently doesn’t have a `loglikelihood` function, only a `logpdf` function, which does the job of both. I thought if I added a `loglikelihood` function for the mixture, it might fix it:

```julia
loglikelihood(d::Union{UnivariateMixture, MultivariateMixture}, x) = logpdf(d, x)

```

But it doesn’t change the error.

Any ideas? Alternatives? I seem to remember that the MixtureModel distribution is a bit of an odd duck that often doesn’t play well with packages like Turing, but if I can avoid writing a custom logpdf function for mixtures, it’d be nice…

Also, the fact that the lack of a `loglikelihood` function for a `MixtureModel` is causing an issue for posterior prediction makes me suspicious of whether the original model estimates are actually correct. (Although I haven’t seen anything that looks obviously wrong.)

---

<div class="post-metadata">

**Author:** ![devmotion](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/devmotion/32/9201_2.png) [@devmotion](https://discourse.julialang.org/u/devmotion)\
**Post date:** [November 18, 2021, 3:14pm UTC](https://discourse.julialang.org/t/using-mixturemodels-from-distributions-jl-in-turing-issues-with-posterior-prediction/71617/2 "2021-11-18T15:14:42Z")

</div>

`loglikelihood` is already defined for `MixtureModel` - it just uses the default `loglikelihood` implementations for `UnivariateDistribution`s and `MultivariateDistribution`s which sums `logpdf` values. The main problem here is that the second argument is of type `Vector{Union{Missing, Float64}}`: Distributions (and also Turing) only supports evaluation of `loglikelihood` with `Real` (for univariate distributions), `AbstractArray{<:Real}`, or e.g. `AbstractArray{<:AbstractVector{<:Real}}` (for multivariate distributions).

More concretely, the problem is that you pass a `Vector{Union{Missing, Float64}}` for `x` in the `predict` call. Generally, Turing samples variables on the LHS of a `~` statement that are not an argument of the model or `missing`. In the expression `x ~ ...` this means `x` would be sampled if `x === missing` (since it is an argument of the model). I.e., you should pass `missing` for `x` in the `predict` call but not a vector of `missing`. Of course, this means that you would have to define `N` independent of `x`. You could e.g. define your model as

```julia
@model function KMM(x, min_x, max_x, k, σ, N=size(x, 1))
    linspan = range(min_x, stop=max_x, length=k)
    kernels = map(u -> Normal(u, σ), linspan)

    ω ~ Dirichlet(k, 1.0)
    mixdist = MixtureModel(kernels, ω)

    x ~ filldist(mixdist, N)
end

```

Then you should be able to perform inference without any changes and could call `predict` e.g. as

```julia
pp_data = predict(KMM(missing, -2.0, 2.0, 10, 0.5, size(data, 1)), m1)

```

---

<div class="post-metadata">

**Author:** ![opera\_malenky](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/opera_malenky/32/8213_2.png) [@opera\_malenky](https://discourse.julialang.org/u/opera_malenky)\
**Post date:** [November 18, 2021, 3:48pm UTC](https://discourse.julialang.org/t/using-mixturemodels-from-distributions-jl-in-turing-issues-with-posterior-prediction/71617/3 "2021-11-18T15:48:23Z")

</div>

Thanks! That worked perfectly.

And yeah, the idea that there was no `loglikelihood` function defined didn’t sound right (otherwise how would the inference work?) But the error – combined with the fact that the docs for the MixtureModel type mention a [`logpdf` function but not `loglikelihood`](https://juliastats.org/Distributions.jl/stable/mixture/#Distributions.logpdf-Tuple%7BAbstractMixtureModel,%20Any%7D) – had me wondering. 🙂