Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

Let $X$ and $Y$ be two Poisson variables with different mean.

Is there a better (as in more concise or numerically faster) way to compute $P(X\leq Y)$ than using

$$ P(X\leq Y) = \sum_{y=0}^\infty P(Y=y) P(X\leq y) $$

share|cite|improve this question
I've been able to reduce it to solving an integral of a product of an exponential and a modified Bessel function evaluated in 2 times the square root of the product of the means. Is that numerically tractable enough? – Raskolnikov Feb 16 '11 at 12:17
Independent Poisson random variables. – Did May 7 '11 at 10:47
up vote 3 down vote accepted

OK, this is how I've done it for now. It's very likely that this is not the best way to do it, but I just decided to post it, for the fun of it. There are probably some mistakes, so go ahead and point them out.

First, for the sake of definiteness, I will posit $X \sim Poiss(\lambda)$ and $Y \sim Poiss(\mu)$. In particular, this means that

$$ \mathbb{P}(X=k) = \frac{e^{-\lambda} \lambda^k}{k!} \; .$$

Then, notice that

$$ \frac{d}{d\lambda} \mathbb{P}(X=k) = \mathbb{P}(X=k-1) - \mathbb{P}(X=k) $$

with special case

$$ \frac{d}{d\lambda} \mathbb{P}(X=0) = - \mathbb{P}(X=0) \; . $$

This implies that

$$ \frac{d}{d\lambda} \mathbb{P}(X \leq y) = - \mathbb{P}(X=y) $$

by summing our first identity over $k$ from $0$ to $y$.

Now, differentiating the sought after probability with respect to $\lambda$ gives us

$$ \frac{d}{d\lambda}\mathbb{P}(X\leq Y) = - \sum_{y=0}^{\infty} \mathbb{P}(X=y)\mathbb{P}(Y=y) $$

The sum on the right hand side can be expressed as

$$ - \sum_{y=0}^{\infty} \mathbb{P}(X=y)\mathbb{P}(Y=y) = - e^{-(\lambda + \mu)} \sum_{y=0}^{\infty} \frac{\lambda^y \mu^y}{(y!)^2} \; . $$

The last part can be rewritten noting that the modified Bessel function $I_0(x)$ has the following Taylor expansion around $0$

$$I_0(x) = \sum_{n=0}^{\infty} \frac{(x/2)^{2n}}{(n!)^2} \; .$$

We then obtain

$$ \frac{d}{d\lambda}\mathbb{P}(X\leq Y) = - e^{-(\lambda + \mu)} I_0(2 \sqrt{\lambda \mu})$$

Now, integrating this we get the following expression

$$ \mathbb{P}(X\leq Y) - \mathbb{P}(0\leq Y)= - \int_0^{\lambda} e^{-(t + \mu)} I_0(2 \sqrt{t \mu}) dt $$


$$ \mathbb{P}(X\leq Y) = 1 - e^{-\mu}\int_0^{\lambda} e^{-t} I_0(2 \sqrt{\mu t}) dt \; . $$

This integral should be in principle more numerically tractable than the previous sum, provided you have modified Bessel functions implemented in whatever program you are using to solve this.

share|cite|improve this answer
Nice stuff! Thanks much, that will indeed be easier to calculate with the tools I have. – Jens Feb 16 '11 at 14:44

Rasholnikov has shown you a method using modified Bessel functions. You could also get that sort of result from a Skellam distribution, which is the difference between two Poisson distributions.

If you want a quick and simple approximation for $X \sim \text{Poiss}(\lambda)$ and $Y \sim \text{Poiss}(\mu)$ then

$$\text{Pr}(X \leq Y) \approx \Phi\left( \frac{1/2 - (\lambda - \mu )}{\sqrt{\lambda + \mu}} \right)$$

where $\Phi()$ is the cumulative distribution function of a standard normal distribution. The $1/2$ is there as a continuity adjustment to deal with the possibility that $X=Y$.

For example with $\lambda =5$ and $\mu =10$, the true result is about 0.9256 while the approximation gives 0.9222.

share|cite|improve this answer
Nice, it confirms my formula which is a special case of the Skellam distribution. – Raskolnikov Feb 16 '11 at 13:17
Thanks! This approximation is indeed.. simpler =) I'll try if its good enough for my case. – Jens Feb 16 '11 at 14:49

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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