Package 'skellam'

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

x, q	For functions dskellam, dskellam.sp, and pskellam.sp: a numeric vector of quantiles.
lambda1, lambda2	Numeric vectors of (non-negative) means; lambda2 defaults to lambda1 if not provided.
log, log.p	Logical; if TRUE, returns the logarithm of the computed value.
lower.tail	Logical; if TRUE (default), returns $P(X \le x)$; otherwise, returns $P(X > x)$.
p	For qskellam: a numeric vector of probabilities.
n	For rskellam: a non-negative integer specifying the number of observations.

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

x	A vector of integers (positive or negative).

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

y	A vector of integers (positive or negative)
x	A matrix, vector or data.frame of covariates

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)

---