The Skellam distribution is the discrete probability distribution of the difference $N_{1}-N_{2}$ of two statistically independent random variables $N_{1}$ and $N_{2}$, each Poisson-distributed with respective expected values $\mu _{1}$ and $\mu _{2}$.

The probability mass function (pmf) of this distribution is

$$p(k;\mu _{1},\mu _{2})=\Pr\{K=k\}=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. Since k is an integer we have that $I_k(z)=I_{|k|}(z)$.

It is known that the mean is $\mu_1-\mu_2$, namely,


In our problem, we have $N_1$ ride requests and $N_2$ vacant vehicles. If $N_{1}\ge N_{2}$, $N_2$ requests will be satisfied, leaving $N_1-N_2$ requests not satisfied, otherwise all requests will be satisfied, leaving 0 unsatisfied request. That is, the number of unsatisfied requests is,


We'd like to get the expectation of $u$, or equivalently,


where $k=N_1-N_2$,

so as to estimate the average number of excess requests.

Since the form of the pmf is kind of complicated, this sum is believed to be hard to calculate, and was wondering if anyone can help out? Thanks.


2 Answers 2


It appears that the direct approach to the calculation, based on the joint pmf of $N_1,N_2$ --
$$E(N_1-N_2)_+=\sum_{x=1}^\infty\sum_{y=0}^{x-1}(x-y)P(N_1=x)P(N_2=y)$$ -- is about 20 times as fast as the approach based on the Bessel formula.

One can reduce the execution time another $5$ times or so using an approach based on

Lemma. Let $X\sim\text{Poisson}(\mu)$ for some real $\mu>0$. Let $h\colon\mathbb R\to\mathbb R$ such that $Eh(X)$ and $Eh(X+1)$ exist in $\mathbb R$ for all real $\mu>0$. Let $g(\mu):=Eh(X)$. Then \begin{equation} g'(\mu)=E[h(X+1)-h(X)]\quad\text{and hence}\quad g''(\mu)=E([h(X+2)-h(X+1)]-[h(X+1)-h(X)]). \end{equation}

Proof. We have $g(\mu)=Eh(X)=\sum_{x=0}^\infty h(x)\frac{\mu^x}{x!}\,e^{-\mu}$ and hence \begin{equation} g'(\mu)=\sum_{x=1}^\infty h(x)\frac{\mu^{x-1}}{(x-1)!}\,e^{-\mu}-\sum_{x=0}^\infty h(x)\frac{\mu^x}{x!}\,e^{-\mu}=Eh(X+1)-Eh(X). \end{equation} $\Box$

Letting now $h_a(x):=(x-a)_+$ for $a=0,1,\dots$, for $g_a(\mu):=Eh_a(X)$ and $X\sim\text{Poisson}(\mu)$ we have \begin{multline} [h_a(X+2)-h_a(X+1)]-[h_a(X+1)-h_a(X)] \\ =I\{X+1\ge a\}-I\{X\ge a\}=I\{X=a-1\}, \end{multline} where $I$ denotes the indicator, and whence $g_a''(\mu)=P(X=a-1)$. Letting now \begin{equation} G(\mu):=G(\mu,\nu):=E(X-Y)_+ =\sum_{a=0}^\infty E(X-a)_+\, P(Y=a) =\sum_{a=0}^\infty g_a(\mu)P(Y=a), \end{equation} where $X\sim\text{Poisson}(\mu)$ and $Y\sim\text{Poisson}(\nu)$ are independent, we have \begin{multline} G''(\mu)=\sum_{a=0}^\infty g''_a(\mu)P(Y=a)= \sum_{a=0}^\infty P(X=a-1)P(Y=a)= P(X=Y-1) \\ =\sum_{x=0}^\infty\frac{\mu^x\nu^{x+1}}{x!(x+1)!}\,e^{-\mu-\nu} =\nu e^{-\mu-\nu}\frac{I_1(2\sqrt{\mu\nu})}{\sqrt{\mu\nu}}, \end{multline} according to the definition of $I_\alpha(x)$. Also, $G(0+)=E(0-Y)_+=0$ and $G'(0+)=P(0\ge Y)=P(Y=0)=e^{-\nu}$. Using the Taylor theorem and the substitution $\sqrt t=:s$, we see that \begin{equation} G(\mu)=G(0+)+G'(0+)\mu+\int_0^\mu(\mu-t)G''(t)\,dt \end{equation} and thus, for $N_1,N_2,\mu_1,\mu_2$ instead of $X,Y,\mu,\nu$,
\begin{equation} E(N_1-N_2)_+=\mu_2 e^{-\mu_1}+2\sqrt{\mu_2}e^{-\mu_2} \int_0^{\sqrt{\mu_1}}(\mu_1-s^2)e^{-s^2}I_1(2s\sqrt{\mu_2})\,ds. \end{equation} The latter expression for $E(N_1-N_2)_+$ may be referred to as the integral Bessel expression, as opposed to the sum Bessel expression based on the Bessel expression for $p(k;\mu_1,\mu_2)$ in the OP.

To see the comparisons of the three approaches to the calculation of $E(N_1-N_2)_+$, click on the image below of a Mathematica notebook once or twice to magnify it:

enter image description here


Based on your description of the problem, I'm going to assume that you're after the first moment of a Skellam distribution truncated at zero and not the partial first moment of a Skellam (which is the expression you give in the title of your question). If my assumption is wrong and you do want a partial moment, the answer is similar with a slight modification, see the end of this answer before the code block.

Then, changing notation slightly, the pmf is: $$ p(u; \lambda_1, \lambda_2) = (\lambda_1/\lambda_2)^{u/2} I_u(2\sqrt{\lambda_1 \lambda_2}) Z_{\lambda_1, \lambda_2}^{-1} $$

Where $Z_{\lambda_1, \lambda_2} = \sum_{k=0}^\infty (\lambda_1/\lambda_2)^{k/2} I_k (2\sqrt{\lambda_1\lambda_2})$ is the normalizing constant, after cancelling out $e^{-\lambda_1-\lambda_2}$.

I follow the same approach as in my answer to this question on Cross-Validated. Briefly,

The moment generating function (mgf) is, $$ \begin{aligned} \mathcal{M}(t; \lambda_1, \lambda_2) &= Z_{\lambda_1, \lambda_2}^{-1} \sum_{u=0}^\infty e^{tu} (\lambda_1/\lambda_2)^{u/2} I_u (2\sqrt{\lambda_1\lambda_2}) \\ &= Z_{\lambda_1, \lambda_2}^{-1} \sum_{u=0}^\infty (\frac{\lambda_1 e^t}{\lambda_2 e^{-t}})^{u/2} I_u (2\sqrt{\lambda_1 e^t \lambda_2 e^{-t}}) \\ &= Z_{\lambda_1, \lambda_2}^{-1} Z_{\lambda_1 e^t, \lambda_2 e^{-t}} \end{aligned} $$

In practice the normalizing constant can be calculated from the easy-to-compute special function Marcum's Q (or, equivalently, from the noncentral-$\chi^2$ cdf, see code block below), $$ Z_{a,b} = Q(\sqrt{2b},\sqrt{2a}) e^{a+b} $$

Differentiating the mgf around $t=0$ and simplifying gives the expectation of the truncated distribution: $$ \begin{aligned} \mathcal{M}'(t; \lambda_1, \lambda_2)|_{t=0} &= \mathbb{E}[u; \lambda_1, \lambda_2] \\ &= (\lambda_1 - \lambda_2) + Z_{\lambda_1,\lambda_2}^{-1} (\lambda_2 I_0 (2\sqrt{\lambda_1 \lambda_2}) + \sqrt{\lambda_1 \lambda_2} I_1(2\sqrt{\lambda_1\lambda_2})) \end{aligned} $$

If you are actually after the partial first moment, just multiply the expectation above through with the normalizing constant (and cancelled terms): $$ \mathbb{E}_{u \geq 0}[u; \lambda_1, \lambda_2] = \mathbb{E}[u; \lambda_1, \lambda_2] \cdot Z_{\lambda_1,\lambda_2} e^{-\lambda_1 - \lambda_2} $$

Example in R:

Q <- function(b,a) 1-pchisq(b^2, 2, a^2)
Z <- function(a,b) Q(sqrt(2*b),sqrt(2*a)) * exp(a+b)

mean_u <- function(l1,l2) 
  (l1 - l2) + (l2*besselI(2*sqrt(l1*l2),0) +

l1 <- 4; l2 <- 12

# [1] 0.9668963

# compare to naively computing series,
u_grid <- 0:100
pmf_u <- sapply(u_grid, function(u) 
                exp(-l1-l2)*(l1/l2)^(u/2) * besselI(2*sqrt(l1*l2),u))
sum(u_grid * pmf_u)/sum(pmf_u)
# [1] 0.9668963

# or (inefficiently) simulating via rejection sampling
k <- rpois(10^7,l1) - rpois(10^7,l2)
# [1] 0.9637044

## alternatively, partial first moment
mean_u(l1,l2) * Z(l1,l2) * exp(-l1-l2)
# [1] 0.02553668

sum(u_grid * pmf_u) #naive summation
# [1] 0.02553668
  • $\begingroup$ It is nice for me to now find out about the existence of Marcum's $Q$. (You apparently used here certain known expressions for the partial derivatives of $Q$.) At mathworld.wolfram.com/MarcumQ-Function.html, two expressions for $Q$ are given: one equivalent to your series expression for $Z$ and the other one an integral expression not unlike the one derived in my answer. Mathematica knows the expression of the noncentral $\chi^2$ cdf in terms of $Q$, but somehow it is unable to express the mentioned series and integral in terms of $Q$. $\endgroup$ Dec 18, 2018 at 16:58
  • $\begingroup$ Once a function is recognized as someone's special and is built into software packages, specialized methods of calculation are discovered and used. It seems the most advanced source at present on the calculation of $Q$ is at arxiv.org/abs/1311.0681 . $\endgroup$ Dec 18, 2018 at 16:59
  • $\begingroup$ Yes, my approach was only to recognize and use known results about $Q$ and derivatives, which I originally encountered via the noncentral $\chi^2$ distribution function. Thanks for the reference wrt to numerical methods. $\endgroup$
    – Nate Pope
    Dec 18, 2018 at 17:50

Your Answer

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

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