MathOverflow is a question and answer site for professional mathematicians. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

The probability mass function for the Skellam distribution for a count difference $k=n_1-n_2$ from two Poisson-distributed variables with means $\mu_1$ and $\mu_2$ is given by:

$$ f(k;\mu_1,\mu_2)= e^{-(\mu_1+\mu_2)} \left({\mu_1\over\mu_2}\right)^{k/2}I_{|k|}(2\sqrt{\mu_1\mu_2}) $$

where $I_k(z)$ is the modified Bessel function of the first kind.

My question: Is it just for convenience to get to grips with the resulting infinite summation terms that the Bessel function appears in this formula or is there a deeper mathematical reason connecting Poisson distributions and Bessel functions or even the Poisson distribution with Bessel differential equations? Is there perhaps even some physical interpretation or intution?

share|cite|improve this question
up vote 8 down vote accepted

A mathematical reason is as follows.

On the one hand, the Laurent series for the modified Bessel functions of the first kind $I_k$ can be deduced from the Laurent series for the Bessel functions of the first kind $J_k$ given here. It reads $$ \sum_{k\in\mathbb{Z}}I_k(x)t^k=\mathrm{e}^{(x/2)(t+1/t)}. $$ On the other hand, the characteristic function of a Poisson random variable $Y$ with mean $\mu$ is $E(z^{Y})=\mathrm{e}^{-\mu(1-z)}$, at least for every complex number $z$ of modulus $1$. Hence the characterization, which you recalled in your post, of Skellam distribution as the distribution of $X=Y_1-Y_2$ for independent Poisson random variables $Y_1$ and $Y_2$ with means $\mu_1$ and $\mu_2$ shows that $$ E(z^X)=\mathrm{e}^{-(\mu_1+\mu_2)}\mathrm{e}^{\mu_1 z+\mu_2/z}. $$ Solving $(x/2)(t+1/t)=\mu_1 z+\mu_2/z$ for $x$ fixed and $t$ depending on $z$ yields $$ x=2\sqrt{\mu_1\mu_2},\qquad t=z\sqrt{\mu_1/\mu_2}. $$ The value of $\mathbb{P}(X=k)$ for every integer $k$ follows, which involves $I_k(2\sqrt{\mu_1\mu_2})$.

Finally, one should not worry too much about the appearance of $I_k$ in this answer versus $I_{|k|}$ in the OP's post because $J_{-k}(-x)=(-1)^kJ_k(x)$ for every integer $k$, hence $I_{-k}(x)=I_k(x)$ (a relation which is also a consequence of the invariance by $t\to1/t$ of the Laurent series for the functions $I_k$).

share|cite|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.