Skellam Laws for Differences of Poisson Counts

Differences of Independent Poisson Counts

Many counts are naturally compared by subtraction: events occurring in one region minus events occurring in another, for example. If the two counts are independent Poisson variables, their difference is integer-valued but can be negative as well as positive. Its distribution is the Skellam distribution. The two Poisson rates determine both the direction of any imbalance and the amount of variability.

Let $X\sim\operatorname{Poisson}(\lambda_1)$ and $Y\sim\operatorname{Poisson}(\lambda_2)$ be independent, where $\lambda_1,\lambda_2>0$, and define $D=X-Y$. The support of $D$ is all the integers: each integer can occur because suitable nonnegative integer values of $X$ and $Y$ can have that difference. We derive its probability mass function by summing over all pairs of counts that produce the same difference.

Definition: For $\lambda_1,\lambda_2>0$, the difference $D=X-Y$ of independent variables $X\sim\operatorname{Poisson}(\lambda_1)$ and $Y\sim\operatorname{Poisson}(\lambda_2)$ has a Skellam distribution with parameters $(\lambda_1,\lambda_2)$. Its probability mass function is $$ P(D=k)=e^{-(\lambda_1+\lambda_2)} \left(\frac{\lambda_1}{\lambda_2}\right)^{k/2} I_{|k|}\!\left(2\sqrt{\lambda_1\lambda_2}\right), \qquad k\in\mathbb Z, $$ where, for a nonnegative integer $m$, the modified Bessel function of the first kind is defined here by $$ I_m(z)=\sum_{j=0}^{\infty} \frac{(z/2)^{2j+m}}{j!(j+m)!}. $$

Deriving the Probability Mass Function

For $k\geq 0$, the event $D=k$ consists of the disjoint possibilities $X=j+k$, $Y=j$, for $j=0,1,2,\ldots$. Independence and the Poisson mass functions therefore give

\begin{aligned} P(D=k) &=\sum_{j=0}^{\infty}P(X=j+k)P(Y=j)\\ &=e^{-(\lambda_1+\lambda_2)} \sum_{j=0}^{\infty} \frac{\lambda_1^{j+k}\lambda_2^j}{(j+k)!j!}\\ &=e^{-(\lambda_1+\lambda_2)} \left(\frac{\lambda_1}{\lambda_2}\right)^{k/2} \sum_{j=0}^{\infty} \frac{(\sqrt{\lambda_1\lambda_2})^{2j+k}}{(j+k)!j!}\\ &=e^{-(\lambda_1+\lambda_2)} \left(\frac{\lambda_1}{\lambda_2}\right)^{k/2} I_k\!\left(2\sqrt{\lambda_1\lambda_2}\right). \end{aligned}

The factorization in the third line follows because $\lambda_1^{j+k}\lambda_2^j =(\lambda_1/\lambda_2)^{k/2}(\sqrt{\lambda_1\lambda_2})^{2j+k}$. The last line is exactly the defining series for $I_k$.

For a negative integer $k=-m$, where $m>0$, the event $D=-m$ consists of $X=j$, $Y=j+m$. Summing those possibilities yields

\begin{aligned} P(D=-m) &=e^{-(\lambda_1+\lambda_2)} \sum_{j=0}^{\infty} \frac{\lambda_1^j\lambda_2^{j+m}}{j!(j+m)!}\\ &=e^{-(\lambda_1+\lambda_2)} \left(\frac{\lambda_1}{\lambda_2}\right)^{-m/2} I_m\!\left(2\sqrt{\lambda_1\lambda_2}\right). \end{aligned}

Combining the nonnegative and negative cases gives the formula with $I_{|k|}$ for every integer $k$. In particular, the absolute value on the Bessel order is essential: its order is nonnegative even when the difference $k$ is negative.

Worked Example: Calculating a Negative Difference Probability

Let $X\sim\operatorname{Poisson}(4)$ and $Y\sim\operatorname{Poisson}(1)$ independently. Find $P(X-Y=-2)$.

Here $k=-2$, so the mass function gives

P(X-Y=-2)=e^{-5}\left(\frac{4}{1}\right)^{-1}I_2(4) =\frac{e^{-5}}{4}I_2(4).

From the defining series,

\begin{aligned} I_2(4) &=\sum_{j=0}^{\infty}\frac{2^{2j+2}}{j!(j+2)!}\\ &=2+\frac{8}{3}+\frac{4}{3}+\frac{16}{45} +\frac{8}{135}+\frac{32}{4725}+\cdots\\ &\approx 6.4222. \end{aligned}

Thus $P(X-Y=-2)\approx (0.0067379)(0.25)(6.4222)=0.01082$. As a check, the direct sum for this event starts with

e^{-5}\left( \frac{1}{2}+\frac{4}{6}+\frac{16}{48} +\frac{64}{720}+\cdots \right) \approx 0.01082.

The two calculations agree; the direct sum enumerates $X=j$, $Y=j+2$, while the Bessel formula packages that infinite sum into one function.

Characteristic Function

A characteristic function provides a compact description of the distribution. For real $t$, the characteristic function of $D$ is $\phi_D(t)=E[e^{itD}]$. Independence lets us factor the expectation, and the Poisson characteristic function follows directly by summing its probability mass function.

Theorem: If $X\sim\operatorname{Poisson}(\lambda_1)$ and $Y\sim\operatorname{Poisson}(\lambda_2)$ are independent, then the characteristic function of $D=X-Y$ is $$ \phi_D(t)= \exp\!\left\{ \lambda_1(e^{it}-1)+\lambda_2(e^{-it}-1) \right\},\qquad t\in\mathbb R. $$

To verify the formula, first use the Poisson mass function to calculate, for any real $t$,

\begin{aligned} E[e^{itX}] &=\sum_{r=0}^{\infty}e^{itr}e^{-\lambda_1}\frac{\lambda_1^r}{r!}\\ &=e^{-\lambda_1}\exp(\lambda_1e^{it}) =\exp\{\lambda_1(e^{it}-1)\}. \end{aligned}

The same calculation for $Y$, with $-t$ in place of $t$, gives $E[e^{-itY}]=\exp\{\lambda_2(e^{-it}-1)\}$. Since $X$ and $Y$ are independent and $e^{itD}=e^{itX}e^{-itY}$,

\begin{aligned} \phi_D(t) &=E[e^{itX}e^{-itY}]\\ &=E[e^{itX}]E[e^{-itY}]\\ &=\exp\!\left\{\lambda_1(e^{it}-1)+\lambda_2(e^{-it}-1)\right\}. \end{aligned}

This calculation establishes the characteristic function for every real $t$, without requiring a moment generating function. Its real and imaginary parts also show the roles of the two rates:

\phi_D(t)= \exp\!\left\{(\lambda_1+\lambda_2)(\cos t-1) +i(\lambda_1-\lambda_2)\sin t\right\}.

The rate difference appears in the phase, while the rate sum determines the real exponential factor. This is consistent with the mean and variance derived next.

Worked Example: Evaluating the Characteristic Function

Suppose $X$ and $Y$ are independent with $X\sim\operatorname{Poisson}(3)$ and $Y\sim\operatorname{Poisson}(2)$. Evaluate the characteristic function at $t=\pi/2$.

Substitute $e^{i\pi/2}=i$ and $e^{-i\pi/2}=-i$:

\begin{aligned} \phi_D(\pi/2) &=\exp\{3(i-1)+2(-i-1)\}\\ &=\exp(-5+i)\\ &=e^{-5}(\cos 1+i\sin 1). \end{aligned}

The modulus is $e^{-5}$, and the argument is $1$ radian. In particular, $\phi_D(-t)$ is the complex conjugate of $\phi_D(t)$, as required for the characteristic function of a real-valued random variable.

Mean and Variance

The difference inherits a simple mean and variance from its two components. For a Poisson variable with rate $\lambda$, the mean and variance are both $\lambda$. Independence makes the covariance between $X$ and $Y$ zero, so subtraction changes the mean by a minus sign but does not subtract the variances.

Theorem: For independent $X\sim\operatorname{Poisson}(\lambda_1)$ and $Y\sim\operatorname{Poisson}(\lambda_2)$, the difference $D=X-Y$ has $$ E[D]=\lambda_1-\lambda_2,\qquad \operatorname{Var}(D)=\lambda_1+\lambda_2. $$

Both Poisson variables have finite second moments, so $D$ also has finite second moment. Linearity of expectation gives

E[D]=E[X-Y]=E[X]-E[Y]=\lambda_1-\lambda_2.

For the variance, use the identity $\operatorname{Var}(X-Y)=\operatorname{Var}(X)+\operatorname{Var}(Y)-2\operatorname{Cov}(X,Y)$. Independence implies $\operatorname{Cov}(X,Y)=0$, and the Poisson variances are their respective rates. Therefore

\operatorname{Var}(D) =\lambda_1+\lambda_2-2(0) =\lambda_1+\lambda_2.

The standard deviation is consequently $\sqrt{\lambda_1+\lambda_2}$. The mean can be positive, negative, or zero, but the variance is positive whenever at least one rate is positive.

Worked Example: Interpreting the Mean and Variance

Let $X\sim\operatorname{Poisson}(8)$ and $Y\sim\operatorname{Poisson}(5)$ independently, where $X-Y$ is the difference in counts over a fixed time interval. The expected difference is

$$ E[X-Y]=8-5=3. $$

The variance and standard deviation are

\operatorname{Var}(X-Y)=8+5=13, \qquad \operatorname{SD}(X-Y)=\sqrt{13}\approx 3.606.

Thus the first count exceeds the second by three on average, but the spread is governed by the sum of the rates, not their difference. A positive expected difference does not mean every observed difference is positive.

Symmetry, Boundary Cases, and Common Pitfalls

When $\lambda_1=\lambda_2=\lambda$, exchanging $X$ and $Y$ does not change their joint distribution, but it changes $D$ into $-D$. Hence the distribution is symmetric about zero. This also follows from the mass formula: the rate-ratio factor is $1$, leaving the same Bessel term for $k$ and $-k$. If the rates differ, the law is generally asymmetric; a larger $\lambda_1$ favors positive differences, and a larger $\lambda_2$ favors negative ones.

Worked Example: Checking Symmetry and the Zero Probability

Let $X$ and $Y$ be independent $\operatorname{Poisson}(2)$ variables. For each integer $k$, the mass function simplifies to

P(X-Y=k)=e^{-4}I_{|k|}(4).

Since $|k|=|-k|$, this gives $P(X-Y=k)=P(X-Y=-k)$. In particular,

P(X-Y=0)=e^{-4}I_0(4)\approx 0.2070.

For a nonzero example, the series gives $I_1(4)=2+4+\frac{8}{3}+\frac{8}{9}+\cdots\approx9.7595$. Therefore $P(X-Y=1)=P(X-Y=-1)=e^{-4}I_1(4)\approx0.1788$, illustrating the symmetry of the probabilities around zero.

If one rate is zero, the ratio in the Bessel-form mass function is undefined, so use the original difference instead. If $\lambda_1=0$, then $X=0$ with probability one and $D=-Y$; if $\lambda_2=0$, then $Y=0$ with probability one and $D=X$. When both rates are zero, $D=0$ with probability one. The characteristic-function and mean-variance formulas remain valid in all these boundary cases.

Common mistake: Subtracting the variances because the random variables are subtracted gives the incorrect expression $\lambda_1-\lambda_2$. Variance measures squared deviation, and the correct identity is $\operatorname{Var}(X-Y)=\operatorname{Var}(X)+\operatorname{Var}(Y)-2\operatorname{Cov}(X,Y)$. For independent counts the covariance term is zero, so the variances add. A related error is using the rate ratio with exponent $|k|/2$ for every $k$; for negative $k$, the exponent must retain the sign of $k$, as shown by the direct sum.

Key takeaway: The difference of independent Poisson counts has a Skellam distribution. Its characteristic function is the product of the two Poisson transforms, with $t$ and $-t$, and its mean and variance are the rate difference and rate sum, respectively.

Check Your Understanding

Use the definitions and derivations above to answer these questions.

If $X\sim\operatorname{Poisson}(6)$ and $Y\sim\operatorname{Poisson}(4)$ independently, what are the mean and variance of $X-Y$?
In the convolution for $P(X-Y=k)$ when $k\geq 0$, which values of $X$ and $Y$ must be summed over?
What changes in the characteristic-function formula when the difference is defined as $Y-X$ instead?
Why is the Skellam mass formula as written for positive rates not directly usable when $\lambda_2=0$? What is the distribution of the difference in that case?
For equal positive rates, explain from the mass function why the probabilities at $k$ and $-k$ agree.