Comparing Poisson variables

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)$$

-
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

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$$

or

$$\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.

-
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.

-
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