# Help with Poisson intermediate parameters

**URL:** <https://discourse.mc-stan.org/t/help-with-poisson-intermediate-parameters/20131>\
**Category:** Modeling\
**Tags:** rstan, fitting-issues, specification, techniques\
**Created:** [January 10, 2021, 4:05am UTC](https://discourse.mc-stan.org/t/help-with-poisson-intermediate-parameters/20131 "2021-01-10T04:05:41Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![pengyu](https://avatars.discourse-cdn.com/v4/letter/p/8baadc/32.png) [@pengyu](https://discourse.mc-stan.org/u/pengyu)\
**Post date:** [January 10, 2021, 4:05am UTC](https://discourse.mc-stan.org/t/help-with-poisson-intermediate-parameters/20131/1 "2021-01-10T04:05:41Z")

</div>

Hi everyone,

I was hoping to include a Poisson variable as intermediate parameter in my model, though I know Stan doesn’t allow integer parameters. I was wondering if log\_sum\_exp or similar functions can be applied on Poisson…

The math expression of this part of model is:  
\lambda\_t \sim \Gamma(\alpha, 1);  
N\_t |\lambda\_t \sim Poisson(\lambda\_t\frac{1-\rho}{\rho});  
\zeta\_t|N\_t \sim \Gamma(N\_t, \beta/\rho)

This is a time series model, with X\_t be observations and model as X\_t = \rho X\_{t-1} + \zeta\_t. If the integer parameter is supported, I was hoping some model looks like:

* * *

```stan
data{
  int<lower=1> T;
  int<lower=0> x[T];
  real<lower=0> rho;
  real<lower=0> alpha;
  real<lower=0> beta;
}
parameters{
  real<lower=0> lambda[T-1];
  int<lower=0> N[T-1];
}
transformed data{
  real<lower=0> zeta[T-1];
  for (t in 2:T){
    zeta[t-1] = x[t] - rho * x[t-1]
  }
}
model {
  lambda ~ gamma(alpha, 1);
  N ~ poisson(lambda);
  zeta ~ gamma(N, beta)
} 

```

Will log\_sum\_exp work on Poisson?

Thanks for any reply!

---

<div class="post-metadata">

**Author:** ![maxbiostat](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/maxbiostat/32/1326_2.png) [@maxbiostat](https://discourse.mc-stan.org/u/maxbiostat)\
**Post date:** [January 10, 2021, 9:15pm UTC](https://discourse.mc-stan.org/t/help-with-poisson-intermediate-parameters/20131/2 "2021-01-10T21:15:18Z")

</div>

Hi,

There are few approaches. First, you might want to do something similar to what @bgoodri does [here](https://discourse.mc-stan.org/t/conway-maxwell-poisson-distribution/2370), which is to sum terms until they’re really small.

Another thing you can do is adaptive truncation, which I discussed in my StanCon2020 [contribution](https://github.com/maxbiostat/StanCon2020/blob/master/poster/poster.pdf).

---

<div class="post-metadata">

**Author:** ![maxbiostat](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/maxbiostat/32/1326_2.png) [@maxbiostat](https://discourse.mc-stan.org/u/maxbiostat)\
**Post date:** [January 12, 2021, 3:50am UTC](https://discourse.mc-stan.org/t/help-with-poisson-intermediate-parameters/20131/3 "2021-01-12T03:50:37Z")

</div>

This might also lend itself to closed-form marginalisation. Let’s do a few calculations and see where we end up. Again, let’s consider the model.

N\_t |\lambda\_t \sim \operatorname{Poisson}(\lambda\_t\frac{1-\rho}{\rho})

\zeta\_t|N\_t \sim \operatorname{Gamma}(N\_t, \beta/\rho)

What we want is the marginal for \zeta\_t, i.e.

f\_Z(\zeta\_t) = \sum\_{k=0}^\infty f\_{Z, N}(\zeta\_t, N\_t = k) = \sum\_{k=0}^\infty f\_{Z|N}(\zeta\_t| N\_t = k)f\_N(N\_t = k).

Substituting the relevant functions in, we get

f\_Z(\zeta\_t) = \sum\_{k=0}^\infty \frac{b^k}{(k-1)!} \zeta\_t^{k-1} \exp(-b\zeta\_t)\frac{\theta\_t^k \exp(-\theta\_t)}{k!},

where b := \beta/\rho and \theta\_t := \lambda\_t(1-\rho)/\rho.

Rerranging we get

f\_Z(\zeta\_t) = \exp(-(b\zeta\_t + \theta\_t))\zeta\_t^{-1}\sum\_{k=0}^\infty \frac{[b\zeta\_t\theta\_t]^k}{(k-1)!k!}.

It turns out that the sum in the last expression has a nice closed-form, which leads to

f\_Z(\zeta\_t) = \exp(-(b\zeta\_t + \theta\_t)) \zeta\_t^{-1}\sqrt{b\zeta\_t\theta\_t} I\_1(2\sqrt{b\zeta\_t\theta\_t}),

where I\_1(\cdot) is the [modified Bessel function of the first kind](https://mathworld.wolfram.com/ModifiedBesselFunctionoftheFirstKind.html). It seems this is now [implemented](https://mc-stan.org/docs/2_25/functions-reference/betafun.html) in Stan, so you might be able to write your marginalisation in closed-form in Stan after all.

EDIT: see @nhuurre’s post.

---

<div class="post-metadata">

**Author:** ![nhuurre](https://avatars.discourse-cdn.com/v4/letter/n/ad7895/32.png) [@nhuurre](https://discourse.mc-stan.org/u/nhuurre)\
**Post date:** [January 12, 2021, 9:27am UTC](https://discourse.mc-stan.org/t/help-with-poisson-intermediate-parameters/20131/4 "2021-01-12T09:27:54Z")

</div>

> [@maxbiostat](#):
>
> Substituting the relevant functions in, we get  
> f\_Z(\zeta\_t) = \sum\_{k=0}^\infty \frac{b^k}{(k-1)!} \zeta\_t^{k-1} \exp(-b\zeta\_t)\frac{\theta\_t^k \exp(-\theta\_t)}{k!},  
> Rerranging we get  
> f\_Z(\zeta\_t) = \exp(-(b\zeta\_t + \theta\_t))\sum\_{k=0}^\infty \frac{[b\zeta\_t\theta\_t]^k}{(k-1)!k!}.  
> It turns out that the sum in the last expression has a nice closed-form, which leads to  
> f\_Z(\zeta\_t) = \exp(-(b\zeta\_t + \theta\_t)) \sqrt{b\zeta\_t\theta\_t} I\_1(2\sqrt{b\zeta\_t\theta\_t}),

Neat, but I think \zeta\_t^{-1} went missing during rearranging.  
Here’s what I get

f\_{Z}\left(\zeta\_{t}\right)=\exp\left(-\theta\_{t}-b\zeta\_{t}\right)\sqrt{\frac{b\theta\_{t}}{\zeta\_{t}}}I\_{1}\left(2\sqrt{b\theta\_{t}\zeta\_{t}}\right)

The Bessel function series

\sum\_{N=1}^\infty \frac{\left(\frac{1}{4}x^2\right)^{N-1}}{(N-1)!\, N!} = \sum\_{k=0}^\infty \frac{\left(\frac{1}{4}x^2\right)^k}{k!\, (k+1)!} = \frac{2}{x} I\_1(x)

corresponds to the N\>0 part of the marginalization. Poisson distribution has non-negligible chance of producing a zero; `gamma(0,beta)` is a degenerate always-zero distribution that needs to be considered separately.

---

<div class="post-metadata">

**Author:** ![maxbiostat](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/maxbiostat/32/1326_2.png) [@maxbiostat](https://discourse.mc-stan.org/u/maxbiostat)\
**Post date:** [January 12, 2021, 11:00am UTC](https://discourse.mc-stan.org/t/help-with-poisson-intermediate-parameters/20131/5 "2021-01-12T11:00:16Z")

</div>

> [@nhuurre](#):
>
> Neat, but I think ζ−1t\zeta\_t^{-1} went missing during rearranging.

Oh yeah. By the time I realised I had forgotten to type that in, I was already in bed. Thanks!

---

<div class="post-metadata">

**Author:** ![pengyu](https://avatars.discourse-cdn.com/v4/letter/p/8baadc/32.png) [@pengyu](https://discourse.mc-stan.org/u/pengyu)\
**Post date:** [January 12, 2021, 8:04pm UTC](https://discourse.mc-stan.org/t/help-with-poisson-intermediate-parameters/20131/6 "2021-01-12T20:04:12Z")

</div>

Thanks @maxbiostat! Sorry about the late response…my PC was sent for repair so I can’t check the forum quite easily…  
I never thought about Bessel function, it looks like a perfect closed form solution. I’ll try coding as soon as my computer comes back. Thanks again for the help!

---

<div class="post-metadata">

**Author:** ![maxbiostat](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/maxbiostat/32/1326_2.png) [@maxbiostat](https://discourse.mc-stan.org/u/maxbiostat)\
**Post date:** [January 21, 2021, 1:53pm UTC](https://discourse.mc-stan.org/t/help-with-poisson-intermediate-parameters/20131/7 "2021-01-21T13:53:19Z")

</div>

> [@nhuurre](#):
>
> Poisson distribution has non-negligible chance of producing a zero; `gamma(0,beta)` is a degenerate always-zero distribution that needs to be considered separately.

To anyone reading this post in the future, I think one way to handle this complication cleanly is to truncate the distribution of N\_t at zero and then the expression becomes

f\_{Z}\left(\zeta\_{t}\right)=\frac{\exp\left(-\theta\_{t}-b\zeta\_{t}\right)}{1-\exp\left(-\lambda\_t\frac{(1-\rho)}{\rho}\right)}\sqrt{\frac{b\theta\_{t}}{\zeta\_{t}}}I\_{1}\left(2\sqrt{b\theta\_{t}\zeta\_{t}}\right)