In particular, for a Skellam distribution obtained by substracting two iid Poisson Processes. Thank you!

  • 1
    $\begingroup$ Taking $\lambda$ as the Poisson mean, I don't know how to evaluate $\sum_{k=1}^\infty k I_{k}(2\lambda)$ -- at least I didn't see an immediate thing to try when looking at the double sum... and Wolfram Alpha didn't succeed with the sum even when fixing $\lambda$ (though cutting it off at finite values well above $\lambda$ yielded accurate numerical values). I wonder if you might get somewhere by looking at the joint distribution of the min and max of two independent Poisson($\lambda$) variates $(X_{(1)},X_{(2)})$ and then looking at $X_{(2)}-X_{(1)}$... ctd $\endgroup$
    – Glen_b
    Commented May 13, 2017 at 9:29
  • $\begingroup$ ctd... though you'll probably end up with the exact same infinite sum to deal with I'd guess. For large $\lambda$ there's a normal approximation for the ordinary difference $X_2-X_1$ that would imply that the mean deviation goes to $\sqrt{\lambda/\pi}$. This works quite well, but it's a bit too large, roughly by a small multiple of $1/\sqrt{\lambda}$; it seems to be a useful bound. $\endgroup$
    – Glen_b
    Commented May 13, 2017 at 9:31
  • $\begingroup$ Very much related (in spite of the title): Is the absolute value of the difference between two Poisson distributions a Poisson distribution? $\endgroup$ Commented May 14, 2020 at 11:14
  • $\begingroup$ Does this answer your question? Is the absolute value of the difference between two Poisson distributions a Poisson distribution? $\endgroup$ Commented Jul 16, 2020 at 20:24

1 Answer 1


It's possible to write the expectation in terms of easy-to-compute special functions.

Let $z$ follow a Skellam distribution with rates $\lambda_1$ and $\lambda_2$, and $k = |z|$. The pmf for $k$ is: $$p(k; \lambda_1, \lambda_2) = \begin{cases} e^{-\lambda_1 - \lambda_2} \left( \left(\frac{\lambda_1}{\lambda_2}\right)^{\frac{k}{2}} I_k(2\sqrt{\lambda_1 \lambda_2}) + \left(\frac{\lambda_2}{\lambda_1}\right)^{\frac{k}{2}} I_{-k}(2\sqrt{\lambda_1 \lambda_2}) \right) &\text{if } k > 0\\ e^{-\lambda_1 - \lambda_2}I_0 (2\sqrt{\lambda_1 \lambda_2})& \text{if } k = 0\end{cases}$$

Here $I_k(a)$ is a modified Bessel function of the first kind and has the symmetry property $I_{k}(a) = I_{-k}(a)$, so the moment generating function of $k$ is

$$\begin{aligned} \mathcal{M}(t; \lambda_1, \lambda_2) = e^{-\lambda_1 - \lambda_2} \left(\sum_{k=0}^{\infty} e^{tk} I_k(2\sqrt{\lambda_1 \lambda_2}) \big[\big(\frac{\lambda_1}{\lambda_2}\big)^{\frac{k}{2}} + \big(\frac{\lambda_2}{\lambda_1}\big)^{\frac{k}{2}} \big] - I_0 (2\sqrt{\lambda_1 \lambda_2}) \right) \end{aligned} $$

Written in this form, recognize that the sum can be written in terms of a special function known as Marcum's $Q$ (used, for example, in the cdf of the noncentral $\chi^2$ distribution). A definition of $Q$ is:

$$ Q(\sqrt{2b},\sqrt{2a}) = e^{-a - b} \sum_{k=0}^\infty \left(\frac{a}{b}\right)^{\frac{k}{2}} I_k(2\sqrt{a b}) $$

So that the moment-generating function becomes:

$$\begin{aligned} \mathcal{M}(t;\lambda_1, \lambda_2) = e^{-\lambda_1 - \lambda_2} \big(&Q(\sqrt{2\lambda_2e^{-t}},\sqrt{2\lambda_1e^t}) e^{\lambda_1e^t + \lambda_2e^{-t}} + \\ &Q(\sqrt{2\lambda_1e^{-t}},\sqrt{2\lambda_2e^t}) e^{\lambda_2e^t + \lambda_1e^{-t}} - \\ &I_0 (2\sqrt{\lambda_1 \lambda_2})\big) \end{aligned}$$

The derivative of $Q(\sqrt{2\lambda_1e^{-t}}, \sqrt{2\lambda_2e^t})$ w.r.t. $t$ is:

$$Q'(\sqrt{2\lambda_1e^{-t}}, \sqrt{2\lambda_2e^t}) = e^{ -\lambda_1 e^t - \lambda_2 e^{-t}} (\lambda_2e^{-t} I_0(2\sqrt{\lambda_1 \lambda_2 }) + \sqrt{\lambda_2 \lambda_1} I_1(2\sqrt{\lambda_1 \lambda_2 }) )$$

Differentiating the mgf around $t=0$ and simplifying gives the expectation of $k$:

$$ \begin{aligned} \mathbb{E}[k; \lambda_1, \lambda_2] = 2 &e^{-\lambda_1 - \lambda_2} \big( \lambda_2 I_0(2\sqrt{\lambda_1 \lambda_2 }) + \sqrt{\lambda_1 \lambda_2} I_1(2\sqrt{\lambda_1 \lambda_2 }) \big) + \\ &(\lambda_2 - \lambda_1)\left(1 - 2 Q(\sqrt{2\lambda_1}, \sqrt{2\lambda_2}) \right) \end{aligned} $$

The $Q$ function can be calculated using any statistical package that implements the noncentral $\chi^2$ distribution function (see below for an R example).

In the case where $\lambda_1 = \lambda_2 = \lambda$, the expectation reduces to:

$$ \mathbb{E}[k; \lambda] = 2\lambda e^{-2\lambda} \left( I_0(2\lambda) + I_1(2\lambda) \right) $$

A numerical example in R:

MarcumQ <- function(a, b) 
  1-pchisq( b^2, 2, a^2) 

# case where l1 \neq l2 
exp_k <- function(l1, l2) {
  a <- 2*sqrt(l1*l2)
  2*exp(-l1-l2)*(l2*besselI(a,0) + a/2*besselI(a,1)) +
       (l2-l1)*(1 - 2*MarcumQ(sqrt(2*l1),sqrt(2*l2)))

exp_k(5,20) # analytical 
# [1] 15.00187 

mean(abs(rpois(100000,5) - rpois(100000,20))) # simulated 
# [1] 15.0018 

# case where l1 = l2 
exp_k2 <- function(l) 
  exp(-2*l)*2*l*(besselI(2*l,0) + besselI(2*l,1)) 

exp_k2(20) # analytical 
# [1] 5.03042 

mean(abs(rpois(100000,20) - rpois(100000,20))) # simulated 
# [1] 5.03498 

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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