Hey everyone,

I’m trying to create a MGMM model from paper http://proceedings.mlr.press/v15/xiong11a/xiong11a.pdf. The idea is that I have a dataset D = [G_1, \dots, G_M], where each group consists of vectors G_M = [x_1, \dots, x_n] and I want to use Turing.jl to try and sample from the model. The algorithm for sampling looks like this (from the paper):

I took inspiration from the tutorial on Gaussian Mixture model and I managed to get the model going for one G_m separately. This model works and returns some estimates. They are not very good, but that is expected since I need to do this for all G_m in D.

```
@model MGMM(Gm) = begin
D, N = size(Gm)
K, T = 3, 2
# draw α and χ
α ~ Dirichlet(T,2)
χ1 ~ Dirichlet(K,2)
χ2 ~ Dirichlet(K,2)
χ = hcat(χ1, χ2)
# choose topic
Y ~ Categorical(α)
θ = χ[:,Y]
# sample mean values and variance
μ1 ~ MvNormal(zeros(D),1)
μ2 ~ MvNormal(zeros(D),1)
μ3 ~ MvNormal(zeros(D),1)
μ = [μ1, μ2, μ3]
σ ~ truncated(Normal(0, 100), 0, Inf)
# draw assignments for each point in x and generate it from a multivariate normal
for i in 1:N
k ~ Categorical(θ[:])
Gm[:,i] ~ MvNormal(μ[k], sqrt(σ))
end
end
model = MGMM(Gm);
gmm_sampler = Gibbs(PG(100, :k, :Y), HMC(0.05, 10, :μ1, :μ2, :μ3, :σ, :α, :χ1, :χ2))
chain = sample(model, gmm_sampler, 10)
```

The problem arrises when trying to use this for my data, which are basically `data = [G1, G2, ..., GM]`

where every G_m is a vector of size 2 \times n, n \sim Poisson(100). I rewrote the model like this:

```
@model MGMM2(data) = begin
K, T = 3, 2
# draw α and χ
α ~ Dirichlet(T,2)
χ1 ~ Dirichlet(K,2)
χ2 ~ Dirichlet(K,2)
χ = hcat(χ1, χ2)
for Gm in data
D, N = size(Gm)
# choose topic
Y ~ Categorical(α)
θ = χ[:,Y]
# sample mean values and variance
μ1 ~ MvNormal(zeros(D),1)
μ2 ~ MvNormal(zeros(D),1)
μ3 ~ MvNormal(zeros(D),1)
μ = [μ1, μ2, μ3]
σ ~ truncated(Normal(0, 10), 0, Inf)
# draw assignments for each point in x and generate it from a multivariate normal
for i in 1:N
k ~ Categorical(θ[:])
Gm[:,i] ~ MvNormal(μ[k], sqrt(σ))
end
end
end
model2 = MGMM2(data);
gmm_sampler = Gibbs(PG(100, :k, :Y), HMC(0.05, 10, :μ1, :μ2, :μ3, :σ, :α, :χ1, :χ2))
chain2 = sample(model2, gmm_sampler, 2)
```

Unfortunatelly, I get `StackOverFlow`

error (showing only the first 2 lines):

```
ERROR: StackOverflowError:
Stacktrace:
[1] MvNormal(::Array{Real,1}, ::PDMats.ScalMat{Float64}) at C:\Users\masen\.julia\packages\Distributions\HjzA0\src\multivariate\mvnormal.jl:202 (repeats 52089 times)
[2] MvNormal(::Array{Real,1}, ::Float64) at C:\Users\masen\.julia\packages\Distributions\HjzA0\src\multivariate\mvnormal.jl:220
etc...
```

I’m not sure where exactly the error comes from and how to possibly try and solve it. Do you have any ideas? I’m a Turing.jl beginner and I’m not even entirely sure that the code is correct.