Comparing Poisson variables

Question

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

Answer 1

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.

Answer 2

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.