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

Question

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

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:

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