How to do posterior draws from a mixture of negative binomials

Modeling
1

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 ?

2 
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
3

It was indeed obvious.
Thank you.