The Skellam Distribution

Overview

This vignette is intended to validate the outputs of the skellam package. Specifically, we verify that the following functions:

• rskellam: Random variate generation from the Skellam distribution.
• dskellam: Evaluation of the Skellam probability mass function.
• qskellam: Calculation of quantiles for the Skellam distribution.

yield results consistent with the distribution of the difference between two independent Poisson random variables. Two sets of parameters are examined:

• Low Lambda Scenario: \(\lambda_1 = 1.5\) and \(\lambda_2 = 0.5\).
• Large Lambda Scenario: \(\lambda_1 = 12\) and \(\lambda_2 = 8\).

Theoretical Background

Consider two independent random variables:

\[ X \sim \mathrm{Poisson}(\lambda_1) \]

and

\[ Y \sim \mathrm{Poisson}(\lambda_2) \]

Then the difference between these two variables is given by:

\[ Z = X - Y \]

It can be shown that \(Z\) follows a Skellam distribution with parameters \(\lambda_1\) and \(\lambda_2\):

\[ Z \sim \mathrm{Skellam}(\lambda_1, \lambda_2) \]

Small Poisson Parameters

In this section, we use a small value for the Poisson parameters. We simulate two Poisson processes, compute their difference, and then compare that result with data generated directly from the Skellam distribution.

N <- 5000
lambda1 <- 1.5
lambda2 <- 0.5

X <- rpois(N, lambda1)
Y <- rpois(N, lambda2)
XminusY <- X - Y
Z <- rskellam(N, lambda1, lambda2)

# Your plot code here
# Density comparison
plot(table(XminusY), xlab = "X - Y", ylab = "", type = "p", pch = 1)
points(table(Z), col = "red", type = "p", pch = 3, cex = 2)
xseq <- seq(floor(par("usr")[1]), ceiling(par("usr")[2]))
points(xseq, N * dskellam(xseq, lambda1, lambda2), col = "blue",
       pch = 4, cex = 3)
legend("topright", pch = c(1, 3, 4), col = c("black", "red", "blue"),
       legend = c("rpois-rpois", "rskellam", "dskellam"))

# Quantile comparison
Sprob <- seq(0, 1, by = 1/100)
qZ <- quantile(Z, prob = Sprob)
plot(qZ, qskellam(Sprob, lambda1, lambda2))
abline(0, 1, col = "#FF000040")

Large Poisson Parameters

Now we repeat the above checks using a scenario with larger Poisson parameters. This serves to ensure that the functions perform correctly across a wider range of parameter values.

lambda1 <- 12
lambda2 <- 8

X <- rpois(N, lambda1)
Y <- rpois(N, lambda2)
XminusY <- X - Y
Z <- rskellam(N, lambda1, lambda2)

# Density comparison
plot(table(XminusY), xlab = "X - Y", ylab = "", type = "p", pch = 1)
points(table(Z), col = "red", type = "p", pch = 3, cex = 2)
xseq <- seq(floor(par("usr")[1]), ceiling(par("usr")[2]))
points(xseq, N * dskellam(xseq, lambda1, lambda2), col = "blue",
       pch = 4, cex = 3)
legend("topright", pch = c(1, 3, 4), col = c("black", "red", "blue"),
       legend = c("rpois-rpois", "rskellam", "dskellam"))

# Quantile comparison
Sprob <- seq(0, 1, by = 1/100)
qZ <- quantile(Z, prob = Sprob)
plot(qZ, qskellam(Sprob, lambda1, lambda2))
abline(0, 1, col = "#FF000040")

Conclusion

This vignette has demonstrated that the functions provided in the skellam package produce outputs consistent with the theoretical behavior of the difference between two independent Poisson random variables. Both for low and large values of the Poisson parameters, the following are confirmed:

• The density of the difference of two Poisson distributions matches that of the Skellam distribution.
• The empirical quantiles of the simulated Skellam data are in line with the theoretical quantiles computed by qskellam.

These checks validate that the rskellam, dskellam, and qskellam functions perform as expected, providing users with a reliable toolset for working with the Skellam distribution.

Getting help

If you encounter a clear bug, please file an issue with a minimal reproducible example on GitHub.