PLAYGROUND

Skellam distribution

Portrait

The Skellam distribution is the difference between two independent Poisson counts. Its integer support extends in both directions: positive values favour the first process and negative values the second.

Historical trail

J. G. Skellam published the difference distribution for Poisson variables from different populations in 1946. The problem appears in ecology, sports, and flow comparisons.

Two useful connections

  • It is a difference of two Poisson variables and approaches Normal when the sum of rates is large. If one rate tends to zero it approaches Poisson or its reflection; Discrete Laplace is a geometric-tail parallel.
  • Even with positive mean, negative differences remain possible unless one rate is nearly zero; variance satisfies Var(X) ≥ |E(X)|.

A use case

Skellam answers “how many more events occurred in A than B?”. It is used for goal differences in sports and for changes between arrival processes in commerce and epidemiology.

Diagnostic advice

Independence is strong: shared exposure or common shocks call for a dependent bivariate model.

A final score hides two attacks

A score of 2 to 1 and another of 5 to 4 leave the same difference, but not the same uncertainty. Under a Skellam model, the mean depends on the difference between the two rates while the variance depends on their sum. A balance between two intense processes can therefore produce much more volatile differences than a balance between sparse processes.

The distribution is used for goal differences, net inventory changes, and subtraction of images with counting noise. If the event total is known, the first count conditional on that total is Binomial; if the total remains free, their difference is Skellam. This dual view helps decide whether to model the net score or its two generating mechanisms separately.

Decision guide

A good candidate when: the variable is the difference between two independent Poisson counts and may take negative values.

Compare it with: Discrete Laplace for a symmetric shape without a Poisson mechanism and a bivariate Poisson model under shared shocks. Dependence changes the variance of the difference.

References

  • SciPy reference: scipy.stats.skellam — definition and parameterization
  • Johnson, N. L., Kemp, A. W. & Kotz, S. (2005). Univariate Discrete Distributions, 3rd ed. Wiley.
  • Skellam, J. G. (1946). The frequency distribution of the difference between two Poisson variates belonging to different populations. Journal of the Royal Statistical Society, Series A, 109(3), 296.
  • Karlis, D. & Ntzoufras, I. (2009). Bayesian modelling of football outcomes: using the Skellam distribution for the goal difference. IMA Journal of Management Mathematics, 20(2), 133–145.

Skellam Distribution: equations and calculator

Distribution defintion

X∼Skellam(λ1,λ2)X\sim\mathrm{Skellam}\left(\lambda_{1},\lambda_{2}\right)

Distribution domain

x∈Z≡{…,−2,−1,0,1,2,… }x\in\mathbb{Z}\equiv\left\{\dots,-2,-1,0,1,2,\dots\right\}

Parameters domain and parameters constraints

λ1∈R+,λ2∈R+\lambda_{1}\in\mathbb{R}^{+},\lambda_{2}\in\mathbb{R}^{+}

Cumulative distribution function

FX(x)=∑k=−∞xfX(k)F_{X}\left(x\right)=\sum_{k=-\infty}^{x}f_{X}\left(k\right)

Probability mass function

fX(x)=e−(λ1+λ2)(λ1λ2)x/2I∣x∣ ⁣(2λ1λ2)f_{X}\left(x\right)=e^{-\left(\lambda_{1}+\lambda_{2}\right)}\left(\frac{\lambda_{1}}{\lambda_{2}}\right)^{x/2}I_{\left|x\right|}\!\left(2\sqrt{\lambda_{1}\lambda_{2}}\right)

Percent point function/Sample

FX−1(u)=arg⁡min⁡x∣FX(x)−u∣F^{-1}_{X}\left(u\right)=\arg\min_{x}\left| F_{X}\left(x\right)-u \right|

Parametric centered moments

E[Xk]=μk′=∑x=−∞∞xkfX(x)E[X^k]=\mu'_{k}=\sum_{x=-\infty}^{\infty}x^{k}f_{X}\left(x\right)

Parametric mean

Mean(X)=μ1′=λ1−λ2\mathrm{Mean}(X)=\mu'_{1}=\lambda_{1}-\lambda_{2}

Parametric variance

Variance(X)=(μ2′−μ1′2)=λ1+λ2\mathrm{Variance}(X)=(\mu'_{2}-\mu'^{2}_{1})=\lambda_{1}+\lambda_{2}

Parametric skewness

Skewness(X)=μ3′−3μ2′μ1′+2μ1′3(μ2′−μ1′2)1.5=λ1−λ2(λ1+λ2)3/2\mathrm{Skewness}(X)=\frac{\mu'_{3}-3\mu'_{2}\mu'_{1}+2\mu'^{3}_{1}}{(\mu'_{2}-\mu'^{2}_{1})^{1.5}}=\frac{\lambda_{1}-\lambda_{2}}{\left(\lambda_{1}+\lambda_{2}\right)^{3/2}}

Parametric kurtosis

Kurtosis(X)=μ4′−4μ1′μ3′+6μ1′2μ2′−3μ1′4(μ2′−μ1′2)2=3+1λ1+λ2\mathrm{Kurtosis}(X)=\frac{\mu'_{4}-4\mu'_{1}\mu'_{3}+6\mu'^{2}_{1}\mu'_{2}-3\mu'^{4}_{1}}{(\mu'_{2}-\mu'^{2}_{1})^{2}}=3+\frac{1}{\lambda_{1}+\lambda_{2}}

Parametric median

Median(X)=FX−1(0.5)\mathrm{Median}(X)=F^{-1}_{X}\left(0.5\right)

Parametric mode

Mode(X)=arg⁡max⁡k∈{⌊λ1−λ2⌋,⌈λ1−λ2⌉}fX(k)\mathrm{Mode}(X)=\arg\max_{k\in\{\lfloor\lambda_{1}-\lambda_{2}\rfloor,\lceil\lambda_{1}-\lambda_{2}\rceil\}}f_{X}(k)

Additional information and definitions

Computing an analytic expression for the inverse of the cumulative distribution function is not feasible. However, it is possible to calculate the Percentile Point Function by approximating it to the nearest integer.
X=N1−N2,  N1∼Poisson(λ1),  N2∼Poisson(λ2),  N1⊥N2X=N_{1}-N_{2},\;N_{1}\sim\mathrm{Poisson}(\lambda_{1}),\;N_{2}\sim\mathrm{Poisson}(\lambda_{2}),\;N_{1}\perp N_{2}
λ1:Rate parameter of N1\lambda_{1}:\text{Rate parameter of }N_{1}
λ2:Rate parameter of N2\lambda_{2}:\text{Rate parameter of }N_{2}
u:Uniform[0,1] random varibleu:\text{Uniform[0,1] random varible}
Iν(x):Modified Bessel function of the first kind of order νI_{\nu}(x):\text{Modified Bessel function of the first kind of order }\nu
⌊x⌋:Floor function\lfloor{x}\rfloor: \text{Floor function}
⌈x⌉:Ceiling Function\lceil{x}\rceil: \text{Ceiling Function}