Say you want to adapt a typical zero inflated poisson model such as:

```
model {
for (n in 1:N) {
if (y[n] == 0) {
target += log_sum_exp(bernoulli_lpmf(1 | theta),
bernoulli_lpmf(0 | theta)
+ poisson_lpmf(y[n] | lambda));
} else {
target += bernoulli_lpmf(0 | theta)
+ poisson_lpmf(y[n] | lambda);
}
}
}
```

to instead have two inflated values, say the 0 and and the 7.

Is the adaptation as simple as the below code?:

```
model {
for (n in 1:N) {
if (y[n] == 0) {
target += log_sum_exp(bernoulli_lpmf(1 | theta0),
bernoulli_lpmf(0 | theta0) +
bernoulli_lpmf(0 | theta7) +
poisson_lpmf(y[n] | lambda));
} else if (y[n] == 7) {
target += log_sum_exp(bernoulli_lpmf(1 | theta7),
bernoulli_lpmf(0 | theta7) +
bernoulli_lpmf(0 | theta0) +
poisson_lpmf(y[n] | lambda));
} else {
target += bernoulli_lpmf(0 | theta0) +
bernoulli_lpmf(0 | theta7) +
poisson_lpmf(y[n] | lambda);
}
}
```

I’ve tried running the model and it seems to converge ok but when I recreate the model in R and run, I seem to be getting an even higher frequency of the inflated values than occurs in the model dataset.

I’m also not quite sure how to write the generated quantities code to create the predictions in Stan. Something like:

```
generated quantities {
int y_test[N];
int <lower=0, upper=1> zero;
int <lower=0, upper=1> seven;
for(i in 1:N) {
y_test[i] = poisson_rng(lambda[i])
zero = bernoulli_rng(theta0[i]);
seven = bernoulli_rng(theta7[i]);
if (zero == 1) {
y_test[i] = 0
} else if (seven == 1) {
y_test[i] = 7
}
}
}
```

It’s the last part I’m really not sure about, especially when both `zero`

and `seven`

equal 1.

Any input would be much appreciated.