# 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[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 ?

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

It was indeed obvious.
Thank you.