# Vectorizing Mixture Poisson model

**URL:** <https://discourse.mc-stan.org/t/vectorizing-mixture-poisson-model/8542>\
**Category:** Modeling\
**Created:** [April 19, 2019, 1:54am UTC](https://discourse.mc-stan.org/t/vectorizing-mixture-poisson-model/8542 "2019-04-19T01:54:22Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![yab](https://avatars.discourse-cdn.com/v4/letter/y/ea666f/32.png) [@yab](https://discourse.mc-stan.org/u/yab)\
**Post date:** [April 19, 2019, 1:54am UTC](https://discourse.mc-stan.org/t/vectorizing-mixture-poisson-model/8542/1 "2019-04-19T01:54:22Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![bbbales2](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bbbales2/32/77_2.png) [@bbbales2](https://discourse.mc-stan.org/u/bbbales2)\
**Post date:** [April 20, 2019, 10:15pm UTC](https://discourse.mc-stan.org/t/vectorizing-mixture-poisson-model/8542/2 "2019-04-20T22:15:34Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![yab](https://avatars.discourse-cdn.com/v4/letter/y/ea666f/32.png) [@yab](https://discourse.mc-stan.org/u/yab)\
**Post date:** [April 22, 2019, 2:29am UTC](https://discourse.mc-stan.org/t/vectorizing-mixture-poisson-model/8542/4 "2019-04-22T02:29:02Z")

</div>

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

```

---

<div class="post-metadata">

**Author:** ![bbbales2](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bbbales2/32/77_2.png) [@bbbales2](https://discourse.mc-stan.org/u/bbbales2)\
**Post date:** [April 22, 2019, 4:39pm UTC](https://discourse.mc-stan.org/t/vectorizing-mixture-poisson-model/8542/5 "2019-04-22T16:39:14Z")

</div>

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

```