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