What is $\sum_{k=0}^{+\infty}{k⋅p(k;\mu_1,\mu_2)}$, where $p$ is the pmf of Skellam distribution?

Question

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,

$$\sum_{k=-\infty}^{+\infty}{k⋅p(k;\mu_1,\mu_2)}=\mu_1-\mu_2$$

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,

$$u=\max{(0,N_1-N_2)}$$

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

$$\sum_{k=0}^{+\infty}{k⋅p(k;\mu_1,\mu_2)}$$

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.

Answer 1

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:

Answer 2

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) +
  sqrt(l1*l2)*besselI(2*sqrt(l1*l2),1))/Z(l1,l2)

l1 <- 4; l2 <- 12

mean_u(l1,l2) 
# [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
set.seed(1)
k <- rpois(10^7,l1) - rpois(10^7,l2)
mean(k[k>=0])
# [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