# Vectorizing Mixture Poisson model

which one is a right way of vectorizing the mixture of Poisson model (truncated at (y > 1)).

    target += log_mix(pi[i],
poisson_lpmf(y[i]|mu[1,i])-log1m_exp(poisson_lccdf(1|mu[1,i])),
poisson_lpmf(y[i]|mu[2,i])-log1m_exp(poisson_lccdf(1|mu[2,i])));


equivalently

   target += log_mix(pi[i],
poisson_lpmf(y[i]|mu[1,i])-log1m_exp(-mu[1,i] + log(1+mu[1,i])),
poisson_lpmf(y[i]|mu[2,i])-log1m_exp(-mu[2,i] + log(1+mu[2,i])));


OR

target += log_mix(pi[i],
poisson_lpmf(y[i]|mu[1,i])-log1m_exp(poisson_lpmf(0|mu[1,i]))
-log1m_exp(poisson_lpmf(1|mu[1,i])),
poisson_lpmf(y[i]|mu[2,i])-log1m_exp(poisson_lpmf(0|mu[2,i]))
-log1m_exp(poisson_lpmf(1|mu[2,i])));


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]);
}
}

If you’re trying to write:

p(y_i | \mu_1, \mu_2) = \pi p_\text{truncated}(y_i | \mu_1) + (1 - \pi) p_\text{truncated}(y | \mu_2)

And:

p_\text{truncated}(y_i | \mu) = \frac{p(y_i | \mu)}{P(y_i > 1 | \mu)}

I think it should be:

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]);
}

1 Like