Hello,

I have tried this implementation.

```
data {
int<lower=1> G;
int y[G]; // data. gene mupression (integer counts)
}
parameters {
positive_ordered[2] mu; // prior location
real mu_zero[2];
real<lower=0> sigma[2];
real<lower=0, upper=1> lambda;
real<lower=0, upper=1> lambda_zero;
}
model {
for(g in 1:G){
if(y[g]==0) {
target += log_mix(lambda_zero,bernoulli_logit_lpmf(0 | mu_zero[1]),
bernoulli_logit_lpmf(0 | mu_zero[2]));
} else {
target += log_mix(lambda_zero,bernoulli_logit_lpmf(1 | mu_zero[1]),
bernoulli_logit_lpmf(1 | mu_zero[2]));
target += log_mix(lambda,
neg_binomial_2_lpmf(y[g] | mu[1], sigma[1]), //NON
neg_binomial_2_lpmf(y[g] | mu[2], sigma[2]) //EXP
);
}
}
sigma ~ gamma(1.1, 2);
mu[1] ~ gamma(2, 0.1);
mu[2] ~ gamma(2, 0.001);
lambda ~ beta(2,2);
mu_zero ~ cauchy(0,2);
}
generated quantities{
int y_hat_1;
int y_hat_2;
y_hat_1 = neg_binomial_2_rng(mu[1], sigma[1]);
y_hat_2 = neg_binomial_2_rng(mu[2], sigma[2]);
}
```

But although the lambda proportion is plausible, the model have n_eff of 1000. What am I doing wrong?

```
$summary
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff
mu[1] 2.013666e+00 3.115304e-05 9.851458e-04 2.011707e+00 2.013001e+00 2.013643e+00 2.014335e+00 2.015687e+00 1000.00000
mu[2] 2.013817e+00 3.108038e-05 9.828478e-04 2.011894e+00 2.013176e+00 2.013792e+00 2.014474e+00 2.015836e+00 1000.00000
mu_zero[1] 4.950094e-01 3.523073e-01 3.287522e+00 -5.046529e+00 -5.889891e-01 -2.383201e-02 9.439835e-01 1.063341e+01 87.07513
mu_zero[2] 4.976698e-01 3.210963e-01 3.631214e+00 -5.006008e+00 -6.284321e-01 -2.724642e-03 8.411252e-01 1.155646e+01 127.88900
sigma[1] 2.231683e-01 1.176291e-04 3.719758e-03 2.157514e-01 2.206539e-01 2.231813e-01 2.256023e-01 2.306998e-01 1000.00000
sigma[2] 3.706016e-04 9.779321e-08 3.092493e-06 3.643050e-04 3.684997e-04 3.707494e-04 3.727237e-04 3.763318e-04 1000.00000
lambda 5.222136e-01 9.710504e-05 3.070731e-03 5.162808e-01 5.201050e-01 5.222999e-01 5.241679e-01 5.284114e-01 1000.00000
lambda_zero 5.042542e-01 2.797127e-02 2.630855e-01 2.747268e-02 3.109872e-01 5.054078e-01 6.916956e-01 9.690703e-01 88.46455
y_hat_1 1.933000e+00 1.307425e-01 4.134441e+00 0.000000e+00 0.000000e+00 0.000000e+00 2.000000e+00 1.500000e+01 1000.00000
y_hat_2 2.088000e+00 2.049091e+00 6.479796e+01 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 1000.00000
lp__ -3.559931e+05 3.327286e-01 2.442393e+00 -3.559989e+05 -3.559945e+05 -3.559928e+05 -3.559912e+05 -3.559894e+05 53.88289
Rhat
mu[1] 1.0027428
mu[2] 1.0031876
mu_zero[1] 1.0672225
mu_zero[2] 1.0570077
sigma[1] 0.9972462
sigma[2] 0.9980175
lambda 0.9984677
lambda_zero 1.0083174
y_hat_1 0.9967248
y_hat_2 1.0000279
lp__ 1.0763487
```

Can the problem be that I have high count values to start with ?

```
> summary(mix[,1])
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0 11.5 73.3 1394.0 277.0 2559000.0
```