5
$\begingroup$

I'm trying to implement the following model in Stan: $$\begin{align} \text{Pr}(y|n,p) & \sim \text{Binomial}(n,p)\\ \text{Pr}(n|\lambda) & \sim \text{Poisson}(\lambda) \end{align}$$

In this model, $y$ and $n$ are non-negative integers, $0<p<1$, $\lambda > 0$, and $y \le n$. Because Stan does not allow for discrete parameters, I need to marginalize over $n$ (i.e., find the probability mass function of $y$ given $p$ and $\lambda$. The PMFs of $y$ and $n$ are as follows:

$$\begin{align} \text{Pr}(y|n,p) = & \frac{n!}{y!(n-y)!}p^{y}(1-p)^{n-y}\\ \text{Pr}(n|\lambda) =& \frac{\lambda^{n}e^{-\lambda}}{n!}\end{align}$$

My attempt to marginalize over $n$ is as follows:

$$\begin{alignat}{4} \text{Pr}(y|p,\lambda) & = \sum_{n=y}^{\infty}\text{Pr}(y|n,p)\text{Pr}(n|\lambda)\\\\ & = \sum_{n=y}^{\infty}\frac{\lambda^{n}e^{-\lambda}}{y!(n-y)!}p^{y}(1-p)^{n-y}\\\\ & = \frac{e^{-\lambda p}p^{y}\lambda^{y}}{\Gamma(y+1)} \end{alignat} $$

The final step was done by Wolfram Alpha, as my own mathematical skills are a bit rusty. Does this seem correct?

To implement this, I first took the log PMF:

$$\text{log-Pr}(y|p,\lambda)=-\lambda p+y(\text{log}p+\text{log}\lambda)-\text{log}\Gamma(y+1)$$

The Stan code for the function is:

  real binomial_poisson_lpmf(int[] y, vector p, vector lambda) { // vectorized;  assumes that y, p, and lambda are all of the same length.
    vector[size(y)] out;
    for(i in 1:size(y))
      out[i] = -lambda[i] * p[i] + y[i] * (log(p[i]) + log(lambda[i])) - lgamma(y[i]+1));
    return(sum(out));
  }

Anyway, I would appreciate any feedback on whether or not this seems valid; I would also like to know if anyone knows of any citations on the subject.

Thank you for your time.

Improve this question
$\endgroup$
3
  • $\begingroup$ If the above is accurate, then in theory this should also be applicable when the Poisson is replaced with a Negative binomial (Mean parameterization PMF): $$\text{Pr}(y|p,\mu,\phi)=\sum_{n=y}^{\infty}\frac{p^{y}(1-p)^{n-y}(n+\phi-1)!}{y!(n-y)!(\phi-1)!}(\frac{\mu}{\mu+\phi})^{n}(\frac{\phi}{\mu+\phi})^{\phi}$$ The resulting compounded PMF would be $$\text{Pr}(y|p,\mu,\phi)=\frac{p^{y}(\frac{\mu}{\mu+\phi})^{y}(\frac{\phi}{\mu+\phi})^{\phi}(\frac{p\mu+\phi}{\mu+\phi})^{-(y+\phi)}\Gamma(y+\phi)}{\Gamma(y+1)\Gamma(\phi)}$$. $\endgroup$
    – C.R. Peterson
    Aug 20, 2016 at 2:49
  • $\begingroup$ You have equalities between things that are not equal. Consider using $\propto$ when necessary. $\endgroup$
    – Glen_b
    Aug 20, 2016 at 5:26
  • $\begingroup$ I suppose $\theta$ and $p$ refer to the same parameter (you may want to edit to fix the notation). $\endgroup$
    – Juho Kokkala
    Aug 20, 2016 at 20:28

1 Answer 1

Reset to default
5
$\begingroup$

$y \mid p,\lambda$ is Poisson!

Your marginalization, or at least the end result, is correct. The form you have obtained for the distribution is the probability mass function of a Poisson distribution -- just write $p^y\lambda^y$ as $(\lambda\,p)^y$ and behold. That is,

\begin{equation} y \mid p, \lambda \sim \mathrm{Poisson}(\lambda\,p). \end{equation}

You have essentially rediscovered the fact that a Poisson process thinned randomly (so that every point is selected with probability $p$ independent of others) is Poisson. This is a well-known result, a quick Google turned up [1] but I suppose this is in many textbooks. Your situation is analogous to this, since the Poisson distributed random variable $y$ can be interpreted as the number of arrivals in a Poisson process with intensity $\lambda$ in a unit interval. Conditional on $y$ the thinning operation selects each of the $y$ points independently with probability $p$, which is a binomial trial.

So, the marginalization could have been 'derived' without any algebraic manipulations by knowing the thinning-of-a-Poisson-process result and realizing how it applies here.

Note about Stan implementation

Note that this also means you do not have to write your own function in Stan since this is simply

y ~ poisson(lambda .* p).

Negative-binomial case

This extends to the NegBin-case (mentioned in comments), too, since a negative binomial can be represented as a mixture of Poissons where the parameter has a Gamma distribution. Conditional on the gamma-distributed parameter, $y$ is Poisson, too. And if $\lambda$ is gamma-distributed, $p\,\lambda$ is too (fixed $p$), so when marginalizing over the gamma-distributed parameter, $y$ is negative-binomial.

The negative-binomial has multiple parametrizations -- depending on the parametrization one has to work out how multiplying the Poisson rates by $p$ translates into the parameters of the negative-binomial. Left as an exercise to the reader.

Reference

[1] http://www.math.uah.edu/stat/poisson/Splitting.html -- Random (formerly Virtual Laboratories in Probability and Statistics), Section 13.5

Improve this answer
$\endgroup$
1
  • $\begingroup$ Ah, I didn't realize that. Your answer was extraordinarily helpful, thank you. $\endgroup$
    – C.R. Peterson
    Aug 22, 2016 at 19:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.