What is the expectation of the absolute value of the Skellam distribution?

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

• 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 Commented May 13, 2017 at 9:29
• 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. Commented May 13, 2017 at 9:31
• Very much related (in spite of the title): Is the absolute value of the difference between two Poisson distributions a Poisson distribution? Commented May 14, 2020 at 11:14
• Does this answer your question? Is the absolute value of the difference between two Poisson distributions a Poisson distribution? Commented Jul 16, 2020 at 20:24

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:

set.seed(4)
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