The Skellam Distribution

Description

Density, distribution function, quantile function, and random generation for the Skellam distribution.

Usage

dskellam(x, lambda1, lambda2 = lambda1, log = FALSE)

dskellam.sp(x, lambda1, lambda2 = lambda1, log = FALSE)

pskellam(q, lambda1, lambda2 = lambda1, lower.tail = TRUE, log.p = FALSE)

pskellam.sp(q, lambda1, lambda2 = lambda1, lower.tail = TRUE, log.p = FALSE)

qskellam(p, lambda1, lambda2 = lambda1, lower.tail = TRUE, log.p = FALSE)

rskellam(n, lambda1, lambda2 = lambda1)

Arguments

Details

The Skellam distribution describes the difference between two independent Poisson random variables. This documentation covers:

Density:

dskellam(x, lambda1, lambda2 = lambda1, log = FALSE)

Distribution Function:

pskellam(q, lambda1, lambda2 = lambda1, lower.tail = TRUE, log.p = FALSE)

Quantile Function:

qskellam(p, lambda1, lambda2 = lambda1, lower.tail = TRUE, log.p = FALSE)

Random Generation:

rskellam(n, lambda1, lambda2 = lambda1)

Saddlepoint Approximations:

dskellam.sp(x, lambda1, lambda2 = lambda1, log = FALSE)
pskellam.sp(q, lambda1, lambda2 = lambda1, lower.tail = TRUE, log.p = FALSE)

If Y_1 and Y_2 are Poisson variables with means \mu_1 and \mu_2 and correlation \rho, then X = Y_1 - Y_2 is Skellam with parameters:

\lambda_1 = \mu_1 - \rho \sqrt{\mu_1 \mu_2}

\lambda_2 = \mu_2 - \rho \sqrt{\mu_1 \mu_2}

The density is given by:

I(2 \sqrt{\lambda_1 \lambda_2}, |x|) (\lambda_1/\lambda_2)^{x/2} \exp(-\lambda_1-\lambda_2)

where I(y,\nu) is the modified Bessel function of the first kind.

Value

• dskellam returns the (log) density.

• pskellam returns the (log) cumulative distribution function.

• qskellam returns the quantile function.

• rskellam generates random deviates.

Invalid lambda values will return NaN with a warning.

Note

The VGAM package also provides Skellam functions. This implementation offers a broader working range, correct handling when one rate parameter is zero, enhanced argument checking, and improved accuracy for x < 0 (in R versions prior to 2.9). Use skellam::dskellam or VGAM::dskellam to specify which implementation to use.

References

• Butler, R. (2007) Saddlepoint Approximations with Applications, Cambridge University Press.

• Johnson, N. L. (1959) On an extension of the connection between Poisson and \chi^2 distributions. Biometrika 46, 352-362.

• Johnson, N. L., Kotz, S., & Kemp, A. W. (1993) Univariate Discrete Distributions, 2nd ed., John Wiley and Sons.

• Skellam, J. G. (1946) The frequency distribution of the difference between two Poisson variates. Journal of the Royal Statistical Society, Series A 109(3), 296.

• Strackee, J. & van der Gon, J. J. D. (1962) The frequency distribution of the difference between two Poisson variates. Statistica Neerlandica 16(1), 17-23.

• Wikipedia: https://en.wikipedia.org/wiki/Skellam_distribution

Examples

# Compare with Poisson when one lambda = 0
dskellam(0:10, 5, 0)
dpois(0:10, 5)

# Both lambdas non-zero
dskellam(c(-1,1), c(12,10), c(10,12))
pskellam(c(-1,0), c(12,10), c(10,12))

# Quantile function
qskellam(c(0.05, 0.95), 3, 4)

# Random generation
rskellam(10, 8.5, 10.25)

Maximum Likelihood Estimation for the Skellam Distribution

Description

Estimates the parameters of a Skellam distribution using maximum likelihood.

Usage

skellam.mle(x)

Arguments

Details

Instead of having to maximize the log-likelihood with respect to both parameters (\lambda_1 and \lambda_2), the function maximizes with respect to \lambda_2 while setting \lambda_1 = \lambda_2 + \bar{x}. This approach improves computational efficiency. The optimization is performed using nlm as it proved faster than optimise.

Value

A list with components:

iters

Number of iterations required by nlm.

loglik

Maximized log-likelihood value.

param

Estimated parameters (\hat{\lambda}_1, \hat{\lambda}_2).

Author(s)

Michail Tsagris

References

• Butler, R. (2007) Saddlepoint Approximations with Applications, Cambridge University Press.

• Johnson, N. L. (1959) On an extension of the connection between Poisson and \chi^2 distributions. Biometrika 46, 352-362.

• Johnson, N. L.; Kotz, S.; Kemp, A. W. (1993) Univariate Discrete Distributions, 2nd ed., John Wiley and Sons.

• 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.

• Strackee, J.; van der Gon, J. J. D. (1962) The frequency distribution of the difference between two Poisson variates. Statistica Neerlandica 16(1), 17-23.

• Abdulhamid, A. A.; Maha, A. O. (2010) On The Poisson Difference Distribution Inference and Applications. Bulletin of the Malaysian Mathematical Sciences Society 33(1), 17-45.

• Wikipedia: Skellam distribution https://en.wikipedia.org/wiki/Skellam_distribution

Examples

# Basic example
x1 <- rpois(1000, 10)
x2 <- rpois(1000, 6)
x <- x1 - x2
skellam.mle(x)

# Larger sample size
x1 <- rpois(10000, 10)
x2 <- rpois(10000, 6)
x <- x1 - x2
skellam.mle(x)

Skellam Regression

Description

Fits a regression model assuming a Skellam distribution for the response variable.

Usage

skellam.reg(y, x)

Arguments

Details

The function uses an exponential link function to ensure positive values for both rate parameters (\lambda_1 and \lambda_2). Optimization is performed using nlm.

Value

A list with components:

loglik

Maximized log-likelihood value

param1

Matrix for \lambda_1 parameters:

• Column 1: Estimated coefficients

• Column 2: Standard errors

• Column 3: t-values (coef/se)

• Column 4: p-values (Wald test)

param2

Matrix for \lambda_2 parameters (same structure as param1)

Author(s)

Michail Tsagris

References

• 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.

• Strackee, J.; van der Gon, J. J. D. (1962) The frequency distribution of the difference between two Poisson variates. Statistica Neerlandica 16(1), 17-23.

• Karlis D. and Ntzoufras I. (2009) Analysis of sports data using bivariate Poisson models. IMA Conference Presentation. http://www2.stat-athens.aueb.gr/~jbn/papers/files/20_Karlis_Ntzoufras_2009_IMA_presentation_handouts_v01.pdf

Examples

set.seed(0)
x <- rnorm(100)
y1 <- rpois(100, exp(1 + 1 * x))
y2 <- rpois(100, exp(-1 + 1 * x))
y <- y2 - y1
skellam.reg(y, x)