Each bacteria splits each second into some number of bacteria that is Poisson distributed with the same parameter for all bacteria. All the splits are not correlated. What would be the distribution of bacteria after $n$ seconds if we start with only one?

My first observation was that $$X_t = \sum_{i=1}^{X_{t-1}} \text{Po}(\lambda) = \text{Po}(X_{t-1}\lambda)$$ And then I attempted to find $X_2$ by conditioning over the poission distribution: $$P(X_2=k) = \sum_{i \geq 0} P(X_2=k | X_1 = i)P(X_1 = i)$$ $$= \sum_{i \geq 0} \frac{e^{-i\lambda}(i\lambda)^k}{k!}\frac{e^{-\lambda}\lambda^i}{i!}= \frac{e^{-\lambda}\lambda^k}{k!} \sum_{i \geq 0} \frac{e^{-i\lambda}i^k\lambda^i}{i!}$$ But then I'm not sure how to evaluate this sum. I was thinking this could maybe be done with MGFs or something similar but I'm not sure how to do that.


1 Answer 1


Interesting, I'll offer a potential solution for the case of $X_2$. I will write as follows so the solution aligns with the notation on the Compound Poisson Distribution on Wikipedia.

$$ Y = \sum_{i=1}^N X_i, N\sim Pois(\lambda),X_i\sim Pois(\lambda) $$

The sums of independent Poisson r.v.s is Poisson as you have used above. $$ Y|N\sim Pois(N\lambda) $$

Consider now the probability generating function of $Y$ unconditionally, and simply marginalize over $N$. The pgf of the Poisson is $\pi(z)=e^{-\lambda(1-z)}$

$$ \pi_Y(z)=\mathbb{E}[z^Y] = \mathbb{E}[\mathbb{E}[z^Y|N]] = \mathbb{E}[e^{-N\lambda(1-z)}] = \sum_{x=0}^\infty e^{-x\lambda(1-z)}\frac{\lambda^x e^{-\lambda}}{x!} $$

Some massaging via Wolfram Alpha, yields the following:

$$ \pi_Y(z) = e^{\lambda(e^{(z-1)\lambda})-1)} = \sum_{x=0}^\infty z^x p(x) $$

Taking derivatives up to order $k$ and evaluating at 0 will yield the desired probability.

$$ \mathbb{P}(Y=k) = \frac{\partial^k}{\partial z^k} \frac{\pi_Y(z)}{k!}\big|_{z=0} $$

Based on the R code below, the derivation above should be correct for this probability mass function where I have computed the mass function and checked it against a simulation for $0,1,2,3$. The higher order derivatives make the subsequent closed form expression prohibitive.

lambda = 5

y = unlist(lapply(c(1:50000),function(x){
  N = rpois(1,lambda)

y_edf = table(y)/length(y)



To derive a feasible and usable pmf, we may consider a Saddlepoint Approximation. Use the fact that $M_Y(t) = \pi_Y(e^t)$, and set $s$ so that $K_Y'(s)=x$. I've overlaid the saddle point approximation with the results from from the simulation, and it matches well. I think it does more poorly with low values of $\lambda$. $$ M_Y(t) = e^{\lambda(e^{(e^t-1)\lambda})-1)}\rightarrow K_Y(t) = \lambda(e^{(e^t-1)\lambda})-1) $$ $$ K'_Y(t) = \lambda^2e^{\lambda (e^{t}-1)+t},K''_Y(t)=\lambda^2e^{\lambda (e^{t}-1)+t}(\lambda e^t+1) $$ $$ \hat{f}(x)=\frac{1}{\sqrt{2\pi K''(s)}}exp(K(s)-sx) $$

cgf <-function(t){

cgf1 <- function(t){

cgf2 <- function(t){

sfind <- function(t,x){

dspois<- function(x){
  s = uniroot(sfind,c(-5,5),x)$root

y_saddle_density = mapply(dspois,c(1:(length(y_edf)-1)))

plot(names(y_edf),y_edf,xlab='Y',ylab='Density'); lines(c(1:length(y_saddle_density)),y_saddle_density)

enter image description here

Edit: Excuse the poor notation for writing density instead of mass.

  • 1
    $\begingroup$ Great job! Just a quick note on the R simulation. The unlisting of lapply's result can be avoided by simply using sapply() instead of lapply(). $\endgroup$
    – Brandmaier
    Jan 24 at 10:02
  • $\begingroup$ @Brandmaier Thanks for the tip! $\endgroup$ Jan 24 at 19:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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