and after we vectorize it in the right way, is it possible to define the generated quantities in this way

generated quantities{
real y_hat_1[N];
real y_hat_2[N];
for (i in 1:N){
y_hat_1[i] = poisson_rng(mu[1,i]); // should T[2,] included here or not?
y_hat_1[i] = poisson_rng(mu[2,i]);
}
}

for(i in 1:size(y)) {
target += log_mix(pi[i],
poisson_lpmf(y[i]|mu[1,i])-poisson_lccdf(1|mu[1,i]),
poisson_lpmf(y[i]|mu[2,i])-poisson_lccdf(1|mu[2,i]));
}

I checked log_mix in the docs and it doesn’t list a vectorized version yet. Do you intend to have a different mixing proportion for every y[i]? Or should that just be a single mu in the log_mix?

Yes, I intend to have a different mixing proportion for every y[i]. If we state the truncation like this, what will be the generated quantities? Is it make sense…?

generated quantities{
real y_hat_1[N];
real y_hat_2[N];
for (i in 1:N){
y_hat_1[i] = poisson_rng(mu[1,i]); // should T[2,] included here or not?
y_hat_1[i] = poisson_rng(mu[2,i]);
}
}

OR

generated quantities{
real y_hat_1[N];
for (i in 1:N){
y_hat_1[i] = log_mix(pi[i],
poisson_rng(mu[1,i]),
poisson_rng(mu[2,i]));
}
}

It depends on what you’re looking for. There are a bunch of different things you could generate and look at.

As I understand it, your generative process is to first randomly pick one of two distributions and then sample from the distribution you picked. This looks like:

int y_hat[N];
for(i in 1:N) {
int which = bernoulli_rng(pi[i]);
y[i] = poisson_rng(mu[which + 1, i]);
}