Help with Poisson intermediate parameters

Modeling
, , ,
1

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:

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!

1 Like
2

Hi,

There are few approaches. First, you might want to do something similar to what @bgoodri does here, 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.

1 Like
3

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. It seems this is now implemented in Stan, so you might be able to write your marginalisation in closed-form in Stan after all.

EDIT: see @nhuurre’s post.

2 Likes
4

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.

1 Like
5

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

6

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!

1 Like
7

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)
1 Like