A Poisson–Lognormal 1D integrator using mode finding (GSL Lambert-W), shifted Gauss–Hermite, or Laplace fallback.

Poisson–Lognormal 1D Distribution

Definitions

The Poisson component for an integer count R ≥ 0 and rate λ > 0 is $$ P_{\rm P}(R \mid \lambda) = e^{-\lambda} \frac{\lambda^{R}}{R!}. $$ The lognormal prior for λ is defined by $$ P_{\rm LN}(\lambda \mid \mu, \sigma) = \frac{1}{\lambda \, \sigma \sqrt{2\pi}} \exp!\left[-\frac{(\ln\lambda - \mu)^2}{2\sigma^2}\right], $$ where $\mu$ is the mean of $\log\lambda$ and $\sigma$ is the standard deviation of $\log\lambda$.

The Poisson–Lognormal likelihood is the integral $$ P(R \mid \mu,\sigma) = \int_{0}^{\infty} P_{\rm P}(R\mid\lambda)\, P_{\rm LN}(\lambda\mid\mu,\sigma)\; \mathrm{d}\lambda. $$

Change of Variables

We rewrite $\lambda$ using the substitution $$ \lambda = e^{\sigma x}, $$ where $x \in (-\infty,\infty)$ becomes the integration variable.

The Jacobian is $$ \mathrm{d}\lambda = \sigma\, e^{\sigma x}\, \mathrm{d}x. $$ Applying this change of variables and simplifying yields the integrand $$ I(x; R, \mu, \sigma) = -\frac{R^2\sigma^2}{2} + u R \sigma - \left(x - u\right)^2 - e^{\sigma x}, $$ where $u = \mu/\sigma + R\sigma$. Thus, $$ P(R\mid\mu,\sigma)= \int_{-\infty}^{\infty} \frac{\exp(I(x;R,\mu,\sigma))}{\sqrt{2\pi}\, R!}\, \mathrm{d}x. $$

Integrand Mode Finding

To find the mode of the integrand, we set the derivative of I with respect to x to zero: $$ \frac{\mathrm{d}}{\mathrm{d}x}I(x;R,\mu,\sigma) = 0. $$ Solving for x gives $$ x = u - \frac{W_0(\sigma^2 \exp(\sigma u))}{\sigma}, $$ Here, $W_0$ is the principal branch of the Lambert $W$ function.