# What is the distribution of a Poisson variable, where the Poisson rate is Normal (or Binomial)?

What is the distribution of $X$ if $$X \sim \text{Poisson}(\lambda), \quad \text{where }\lambda \sim N(\mu,\sigma^2)$$

or

$$X \sim \text{Poisson} (\lambda), \quad \text{where }\lambda \sim Bin(n,p)$$ where $n$ is very large.

# Detail and motivation

In biology, next-generation sequencing is often used to quantify the level of RNA in a cell. The experiment goes like this:

Imagine that a cell (or a extract made from many cells), is a mixed bag of 200,000 RNAs molecules from 20,000 classes. We wish to estimate the proportions of the various classes.

For each replicate of an experiment we sample, with replacement, RNA "reads" from the bag.

For a class of RNA we get a count of "reads".

Each replicate of an experiment is performed with RNA from a different cell (or an extract from a different group of many cells), drawn from the same underlying population of cells.

If X is the count of "reads" from a particular replicate of the experiment, in practice we often model X as being negatively binomial as when we look at it, it looks like an overdispersed poisson.

There has been some discussion recently if their is any theoretical basis behind using the negative bionomial other than "it seems to fit". One obvious answer is that the negative binomial is the result of a poisson distribution where the poisson rate follows a gamma distribution. One interpretation is that our measurement error is poisson (counts from a sample), and the underlying variance between biological samples is gamma. But why gamma: it does not seem like an obvious answer as to the distribution to use.

A model that we are more used to might be that the amount of RNA in a sample is normally distributed (or equivalently binomially distributed with large n).

# Previous attempts

I'm sure I've seen this written down somewhere before, using PGFs to figure out the result.

$$\Pi_x(s) = \sum_x{P(X=x)s^x}$$

$$P(X=x) = \int_r{P(X=x|\lambda=r)P(\lambda=r)}$$

but writing this down has only resulted in an intergral I don't know how to solve!

• 1. The Poission $\lambda$ parameter is required to be positive, so some care is needed if you want to assume that $\lambda$ is Normal, to ensure the Normail tails do not practically become negative. _______ 2. One can solve the parameter-mix distribution for your desired Poisson-Normal as a closed form, but it is a bit messy with Hypergeometric1F1 functions. ___3. What is the meaning of your title: TL:DR? – wolfies Nov 21 '17 at 13:42
• @wolfies it stands for "too long; didn't read." – Taylor Nov 21 '17 at 14:14
• You write "A model that we are more used to might be [...] binomially distributed with large n)": And that is Poisson. Citing from Wikipedia for simplicity: "The Poisson distribution can be derived as a limiting case to the binomial distribution as the number of trials goes to infinity and the expected number of successes remains fixed" or en.wikipedia.org/wiki/Poisson_distribution#law_of_rare_events . Binomial and Poisson are always positive. – Bernhard Nov 21 '17 at 14:21
• A Poisson parameter $\lambda$ cannot be negative hence cannot be Normal. – Xi'an Nov 21 '17 at 14:44
• Regarding the negbin: On the one hand, I personally am good with "because it works". On the other hand, Hilbe's Negative Binomial Regression lists many, many motivations and constructions that all lead to the negbin. (But that doesn't have much to do with your original question, which I find interesting and wait for an answer on.) – Stephan Kolassa Nov 21 '17 at 18:17

## 2 Answers

So you assume a continuous mixture of Poisson distributions, say $$X \mid \Lambda=\lambda \sim \mathcal{Pois}(\lambda),$$ where $$\Lambda$$ then has some distribution on the positive line.

If we assume that $$\Lambda$$ has a gamma distribution, results that $$X$$ has a marginal negative binomial distribution. You asks what happens if $$\Lambda$$ has a normal distribution. That doesn't seem like a natural model since a normal distribution can take negative values, which does not make sense here. Also, the proposal of a binomial mixture distribution seems strange, since why restrict $$\lambda$$ to integers? Here is a review paper on continuous Poisson mixtures you could have a look at it and try some of the more prominent examples. Table 1 page 50 in the cited paper is a list of examples, many of them difficult to use due to numerical issues.

An interesting idea would be to estimate the mixing density $$g$$ on $$(0,\infty)$$ nonparametrically. There are many possible methods and many papers.

### Normal approximation

One approach is to use a normal approximation by matching the moments.

The compound Poisson distribution is a mixture distribution of Poisson distributions with weights distributed according to a binomial distribution $$w_i = P_{binom}(i\vert n,p)$$

$$P(x) = \sum_{i=0}^n w_iP_{Poisson}(x\vert \lambda = i )$$

• The mean of the distribution will be equal to $$E[X] = \sum_{i=0}^n w_i E[X\vert i] = \sum w_i i = np$$ the mean of the binomial distribution.
• The variance will be $$\begin{array}{} E[(X-\mu_X)^2] &=& \sum_{i=0}^n w_i (\sigma_i^2 + \mu_i^2 - \mu^2) \\& = &\sum_{i=0}^n w_i (i + i^2 - (np)^2) \\&=& \sum_{i=0}^n w_i (i + i^2) - \sum_{i=0}^n (np)^2\\ &=& np(1-p) + (np) \end{array}$$ the variance plus the mean of the binomial distribution

### Negative binomial approximation

We can do the same sort of approximation (matching the moments) with a negative binomial distribution.

In the examples below, we have plotted this approximation as well. The curve resembles more closely the histogram than the normal approximation, but it is still not exactly the same.

### Computational

The example below is computed with $$n=40$$ and $$p=0.5$$ giving $$\mu_X = 20$$ and $$\sigma_X^2 = 30$$ ## parameters
set.seed(1)
ns <- 10^5
nbinom <- 40
pbinom <- 0.5

### generate sample
Y <- rbinom(ns, nbinom, pbinom)
X <- rpois(ns, Y)

### plot histogram
h <- hist(X, breaks = seq(0,max(X)+1)-0.5, freq = 0)

### add curve
xs <- seq(0,max(X))
#lines(xs,dpois(xs,nbinom*pbinom)*ns, col = 2)
lines(xs,dnorm(xs+0.5, mean = nbinom*pbinom,
sd = sqrt(nbinom*pbinom*(2-pbinom))), col = 2)

### binomial approximation
mu = nbinom*pbinom
vr = (nbinom*pbinom*(2-pbinom))
p = 1-mu/vr
r = mu*(1-p)/p
lines(xs,dnbinom(xs, size = r, p= 1-p ), col = 3)

legend(max(X),max(h\$density), c("normal approximation", "negative binomial approximation"),
cex= 0.7, col = c(2,3),lty = 1, xjust = 1)


The problem with the normal approximation is that it does not work so well when $$\sigma_X$$ is not a lot smaller than $$\mu_X$$. This occurs when $$\mu_X = np$$ is small.

See for instance the plot below with $$n=4000$$ and $$p=0.0005$$ giving $$\mu_X = 2$$ and $$\sigma_X^2 = 1.99975$$ 