4

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.

mrepic1123
  • 489
  • 6

1 Answers1

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)
  sum(rpois(N,lambda))
}))


y_edf = table(y)/length(y)


head(y_edf,4)


exp((exp(-lambda)-1)lambda)
exp((exp(-lambda)-2)lambda)lambda^2
exp((exp(-lambda)-2)lambda)lambda^3(exp(-lambda)lambda+1)/factorial(2)
exp(lambda(exp(-lambda)-2))lambda^4(exp(-2lambda)lambda^2+3exp(-lambda)lambda+1)/factorial(3)

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){
  lambda*(exp((exp(t)-1)*lambda)-1)
}

cgf1 <- function(t){
  lambda^2exp(lambda(exp(t)-1)+t)
}


cgf2 <- function(t){
  lambda^2exp(lambda(exp(t)-1)+t)(lambdaexp(t)+1)
}


sfind <- function(t,x){
  cgf1(t)-x
}


dspois<- function(x){
  s = uniroot(sfind,c(-5,5),x)$root
  exp(cgf(s)-sx)/sqrt(2pi*cgf2(s))
}


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.

user1848065
  • 345