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