I’m a new Stan user but I was not able to resolve my basic question by reading the manual or searching through the forums. I guess I’m missing something obvious.

I’m trying to fit a mixture a negative binomials to the distribution of the count of fail events in a system.

Here is my stan program :

```
data {
int<lower=1> N;
int fail[N]; // data (integer counts)
}
parameters {
positive_ordered[2] mu;
real<lower=0> phi[2];
real<lower=0, upper=1> lambda;
}
model {
phi[1] ~ gamma(2, 4);
mu[1] ~ gamma(5, 0.008);
phi[2] ~ gamma(8, 5);
mu[2] ~ gamma(8, 0.0001);
lambda ~ beta(2,2);
for(g in 1:N){
target += log_mix(lambda,
neg_binomial_2_lpmf(fail[g] | mu[1], phi[1]),
neg_binomial_2_lpmf(fail[g] | mu[2], phi[2])
);
}
}
generated quantities{
real y_rep;
int y_hat_1;
int y_hat_2;
y_hat_1 = neg_binomial_2_rng(mu[1], phi[1]); // first component of the mixture
y_hat_2 = neg_binomial_2_rng(mu[2], phi[2]); // second component of the mixture
//y_rep = log_mix(lambda, y_hat_1, y_hat_2); // My issue is here
}
```

The model seems to converge correctly and inspecting the y_hat distribution confirm this result

```
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
mu[1] 1431.58 1.44 70.21 1303.78 1382.13 1427.65 1478.18 1577.25 2388 1
mu[2] 55837.46 12.21 672.51 54508.23 55378.65 55834.62 56295.59 57143.64 3036 1
phi[1] 1.27 0.00 0.07 1.15 1.22 1.27 1.32 1.42 2515 1
phi[2] 3.15 0.00 0.14 2.90 3.06 3.15 3.25 3.43 2444 1
lambda 0.31 0.00 0.01 0.29 0.30 0.31 0.31 0.32 2327 1
```

but I am unable to draw samples from the posterior distribution in the last line of the generated quantities block.

I tried to adapt the example p190 in the reference manual (*Inference under Label Switching*), in vain.

What am I missing ?