# How to do posterior draws from a mixture of negative binomials

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 mu;
real<lower=0> phi;
real<lower=0, upper=1> lambda;
}
model {
phi ~ gamma(2, 4);
mu ~ gamma(5, 0.008);
phi ~ gamma(8, 5);
mu ~ gamma(8, 0.0001);
lambda ~ beta(2,2);

for(g in 1:N){
target += log_mix(lambda,
neg_binomial_2_lpmf(fail[g] | mu, phi),
neg_binomial_2_lpmf(fail[g] | mu, phi)
);
}
}
generated quantities{
real y_rep;
int y_hat_1;
int y_hat_2;

y_hat_1 = neg_binomial_2_rng(mu, phi); // first component of the mixture
y_hat_2 = neg_binomial_2_rng(mu, phi); // 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     1431.58    1.44    70.21   1303.78   1382.13   1427.65   1478.18   1577.25  2388    1
mu    55837.46   12.21   672.51  54508.23  55378.65  55834.62  56295.59  57143.64  3036    1
phi       1.27    0.00     0.07      1.15      1.22      1.27      1.32      1.42  2515    1
phi       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 ?

``````if(bernoulli_rng(lambda))
y_rep = neg_binomial_2_rng(mu, phi);
else
y_rep = neg_binomial_2_rng(mu, phi);
``````
1 Like

It was indeed obvious.
Thank you.