Hello!

I’m new to Julia but very excited by the possibilities. I’m particularly interested in using Turing to fit mixture models. However, in building up some non-trivial mixture examples, I’ve encountered an issue I just can’t figure. I wrote up the same model in R and Julia itself and it works fine in both cases so this is something about my understanding about how Turing works.

I have stripped down version of a binomial mixture model with multiple observations (i.e. each sample contains Q distinct observations, each its own binomial). To further reduce the problem, I feed the model the correct parameters except for the indicator variables. However, I cannot get Turing to sample the indicator parameters correctly. What is going wrong?! Any help would be greatly appreciated.

Here is the code for the model itself.

```
@model function binomial_replicate_mixture_model_p(t_counts,W,K,p,lambda)
N,Q = size(t_counts)
# Priors for the mixture weights
#lambda ~ Dirichlet(K, 1.0)
# Priors for the multinomial parameters
z ~ filldist(Categorical(lambda),N)
for n in 1:N
for q in 1:Q
t_counts[n,q,:] ~ Multinomial(W[n,q],vec([1-p[z[n],q],p[z[n],q]]))
end
end
end
```

## Simulation and control script

Here is the code for the data simulation:

```
using Distributions,LinearAlgebra
using CSV,DataFrames
using Random,Turing,MCMCChains
function sim_bmm(N,Q,K)
v = rand(1:K,N)
y_counts = zeros(Int64,N,Q)
n_counts = zeros(Int64,N,Q)
alpha = Matrix{Float64}(undef, K,Q)
beta = Matrix{Float64}(undef, K,Q)
p = Matrix{Float64}(undef, K,Q)
for k in 1:K
for q in 1:Q
alpha[k,q] = rand([3.0,10.0,30.0])
beta[k,q] = rand([3.0,10.0,30.0])
p[k,q] = alpha[k,q]./(alpha[k,q]+beta[k,q])
end
end
#p_mat = alpha ./ (alpha + beta)
for i in 1:N
for q in 1:Q
tmp = rand(Multinomial(10000,vec([p[v[i],q],1-p[v[i],q]])),1)
y_counts[i,q] = tmp[1]
n_counts[i,q] = tmp[2]
end
end
return y_counts,n_counts,p,v
end
```

And here’s the control script:

```
## simulation parameters
N = 100 ## num of samples
Q = 5 ## num of questions
K = 7 ## num of components
y_c,n_c,p,v = sim_bmm(N,Q,K)
t_c = cat(y_c, n_c, dims=3)
W = y_c+n_c
lambda = fill(1/K,K)
## sampler parameters
nsamples=5000
nchains=1
sampler_p = Gibbs(MH(:z))
## model definition
bbmm_p = binomial_replicate_mixture_model_p(t_c,W,K,p,lambda)
chain_p = sample(bbmm_p,sampler_p,MCMCThreads(),nsamples,nchains);
```