I have a Poisson process with a fixed (large) $\lambda$. If I run the process twice, what is the probability that the two runs have the same outcome?
That is, how can I approximate $$f(\lambda)=e^{-2\lambda}\sum_{k=0}^\infty\frac{\lambda^{2k}}{k!^2}$$ for $\lambda\gg1$? If there's a simple expression about $+\infty$ that would be best, but I'm open to whatever can be suggested.
Fourier transforms yield a fully rigorous proof.
First recall that, as explained here, for every integer valued random variable $Z$, $$ P(Z=0)=\int_{-1/2}^{1/2}E(\mathrm{e}^{2\mathrm{i}\pi tZ})\mathrm{d}t. $$ Hence, if $X_\lambda$ and $Y_\lambda$ are independent Poisson random variables with parameter $\lambda$, $$ f(\lambda)=P(X_\lambda=Y_\lambda)=\int_{-1/2}^{1/2}E(\mathrm{e}^{2\mathrm{i}\pi tX_\lambda})E(\mathrm{e}^{-2\mathrm{i}\pi tY_\lambda})\mathrm{d}t. $$ For Poisson distributions, one knows that $E(s^{X_\lambda})=\mathrm{e}^{-\lambda(1-s)}$ for every complex number $s$. This yields $$ f(\lambda)=\int_{-1/2}^{1/2}\mathrm{e}^{-2\lambda(1-\cos(2\pi t))}\mathrm{d}t=\int_{-1/2}^{1/2}\mathrm{e}^{-4\lambda\sin^2(\pi t)}\mathrm{d}t. $$ Consider the change of variable $u=2\pi\sqrt{2\lambda}t$. One gets $$ f(\lambda)=\frac1{\sqrt{4\pi\lambda}}\int_\mathbb{R} g_\lambda(u)\mathrm{d}u, $$ with $$ g_\lambda(u)=\frac1{\sqrt{2\pi}}\mathrm{e}^{-4\lambda\sin^2(u/\sqrt{8\lambda})}\,[|u|\le\pi\sqrt{2\lambda}]. $$ When $\lambda\to+\infty$, $g_\lambda(u)\to g(u)$ where $g$ is the standard Gaussian density, defined by $$ g(u)=\frac1{\sqrt{2\pi}}\mathrm{e}^{-u^2/2}. $$ Furthermore, the inequality $$4\lambda\sin^2(u/\sqrt{8\lambda})\ge2u^2/\pi^2, $$ valid for every $|u|\le\pi\sqrt{2\lambda}$, shows that the functions $g_\lambda$ are uniformly dominated by an integrable function. Lebesgue dominated convergence theorem and the fact that $g$ is a probability density yield finally $$ \int_\mathbb{R} g_\lambda(u)\mathrm{d}u\to1,\qquad\text{hence}\ \sqrt{4\pi\lambda}f(\lambda)\to1. $$
So I've come across this problem before, and ran into the Bessel-function exact solution; but I figured that there must be some heuristic explanation. Here's the best I was able to come up with.
Let $X$ and $Y$ be independent Poisson with mean $\lambda$ (and therefore variance $\lambda$). You want $P(X=Y)$.
Now, $X$ and $Y$ are both approximately normal with mean $\lambda$ and variance $\lambda$, if $\lambda$ is large.
Let $U = X + \lambda$ and $V = 3 \lambda - Y$. Then $U$ and $V$ are both approximately normal with mean $2 \lambda$ and variance $\lambda$. But these are the mean and variance of a binomial distribution with parameters $4\lambda$ and $1/2$. So $U$ and $V$ are both approximately binomial with parameters $4\lambda$ and $1/2$; thus $U + V$ is approximately binomial with parameters $8\lambda$ and $1/2$. But by construction $U + V = 4 \lambda + X - Y$. So $X = Y$ if and only if $U + V = 4 \lambda$.
The probability that $U + V$ is $4 \lambda$ is approximately ${8 \lambda \choose 4 \lambda} 2^{-8\lambda}$. And by Stirling's approximation this is of order $(4 \pi \lambda)^{-1/2}$, which is what we wanted.
(The fact that the probability should decay like $\lambda^{-1/2}$ is fairly easy to see -- since $X$ and $Y$ have standard deviation $\lambda^{1/2}$, the number of values they take ``regularly'' is of order $\lambda^{1/2}$, so the chance of collision is of order $\lambda^{-1/2}$. The constant $(4\pi)^{-1/2}$ is in my opinion much harder to guess.)
Another way, similar to Michael, more quick and informal, to get an asymptotic:
We have two independent poisson variables $X$ and $Y$, with parameter $\lambda \gg 1$ and we want to estimate $P(Z=0)$ where $Z=X-Y$. But using the normal approximation, we have $X \approx N(\lambda,\lambda)$, and so $Y$; and hence $Z \approx N(0,2\lambda)$ ($N$ here denotes the Normal distribution). We are approximating here a discrete variable by a continuous one, and so we can further approximate:
$$ P(Z=0) \approx \int_{-1/2}^{1/2} f_z(Z) dZ \approx f_z(0) = N(0,2\lambda)|_{z=0} = \frac{1}{\sqrt{4 \pi \lambda}} $$
According to Wolfram alpha, $$ f(\lambda) = e^{-2\lambda} I_0(2\lambda). $$ It also gives a first-order approximation $$ f(\lambda) = \frac{1}{\sqrt{4\pi\lambda}} + O(\lambda^{-3/2}). $$
