Package {skellam}

The Skellam Distribution

Description

Density, distribution function, quantile function and random generation for the Skellam distribution, the distribution of the difference of two independent Poisson variables, together with saddlepoint approximations to the density and distribution function.

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

If Y_1 \sim \mathrm{Poisson}(\lambda_1) and Y_2 \sim \mathrm{Poisson}(\lambda_2) are independent, then X = Y_1 - Y_2 has the Skellam distribution with probability mass function

P(X = x) = e^{-(\lambda_1+\lambda_2)} (\lambda_1/\lambda_2)^{x/2} I_{|x|}(2\sqrt{\lambda_1 \lambda_2}), \qquad x = \ldots, -1, 0, 1, \ldots

where I_\nu is the modified Bessel function of the first kind. The mean is \lambda_1 - \lambda_2 and the variance \lambda_1 + \lambda_2. (If Y_1 and Y_2 are correlated Poisson variables with means \mu_i and correlation \rho, the difference is still Skellam, with \lambda_i = \mu_i - \rho\sqrt{\mu_1 \mu_2}.)

All functions follow the conventions of the stats distribution functions: arguments are recycled to a common length (with a warning when the lengths are not multiples), missing values propagate, invalid rates and probabilities give NaN with a warning, a non-integer x has density 0 (with a warning), and qskellam returns the smallest integer x with P(X \le x) \ge p (for lower.tail = TRUE), so that qskellam(pskellam(x, ...), ...) recovers x.

Numerical methods. dskellam evaluates the log of the Bessel function directly, so that log = TRUE is finite and accurate for any input: through besselI where that is reliable, through the uniform asymptotic expansion for large orders, a stable downward recurrence for large arguments, and the power series where besselI underflows. The relative error is below 10^{-12} for rates up to at least 10^7, and below 10^{-9} for extremely unequal large rates such as (10^6, 10^3). pskellam uses the non-central chi-square identity of Johnson (1959) through pchisq when \lambda_1 + \lambda_2 \le 10^4, and otherwise, as well as in tails where pchisq underflows, sums the log-pmf directly toward the nearer tail, so that far tails are exact on the log scale. qskellam starts from a Cornish-Fisher approximation and searches the integers with pskellam. rskellam returns the difference of two rpois draws.

dskellam.sp and pskellam.sp are the second-order saddlepoint approximation to the pmf and the Lugannani-Rice approximation to the cdf with Butler's second continuity correction (Butler 2007). Their relative error is of order 1/(\lambda_1+\lambda_2)^2 and 1/(\lambda_1+\lambda_2) respectively; they are provided for users who want the approximations themselves and are not used by dskellam or pskellam.

Value

dskellam gives the (log) probability mass function, pskellam the (log) distribution function, qskellam the quantile function and rskellam generates random deviates (an integer vector). dskellam.sp and pskellam.sp give the saddlepoint approximations. Names of the first argument are kept.

Invalid lambda values return NaN with a warning; a zero rate is allowed and reduces the distribution to a (reflected) Poisson.

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.

• Olver, F. W. J. et al. (eds.) NIST Digital Library of Mathematical Functions, Section 10.41 (uniform asymptotic expansions of Bessel functions). https://dlmf.nist.gov/10.41

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

See Also

skellam.mle and skellam.reg for estimation; dpois for the Poisson 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 inverts the distribution function
qskellam(c(0.05, 0.95), 3, 4)
qskellam(pskellam(-3:3, 3, 4), 3, 4)

# Log-scale values stay finite far in the tails and at large rates
dskellam(500, 3, 2, log = TRUE)
pskellam(0, 1e6, 1e6)

# Random generation
set.seed(1)
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

At the maximum likelihood estimate the fitted mean equals the sample mean, \hat\lambda_1 - \hat\lambda_2 = \bar{x}, so the function maximizes the profile log-likelihood over \lambda_2 alone, with \lambda_1 = \lambda_2 + \bar{x}, using optimize on the admissible range \lambda_2 \ge \max(0, -\bar{x}), which keeps both rates non-negative.

The estimate can lie on the boundary, with one rate equal to 0 (the data are then best described by a Poisson or a reflected Poisson). This happens when the sample variance does not exceed the absolute sample mean, which a Skellam with two positive rates cannot produce (its variance is \lambda_1 + \lambda_2 \ge |\lambda_1 - \lambda_2|); for example, a sample of zeros gives (0, 0) and a constant sample c > 0 gives (c, 0).

Value

A list with components:

iters

Number of log-likelihood evaluations used by the optimizer.

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
set.seed(1)
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 response is modelled as y_i = Y_{1i} - Y_{2i} with Y_{1i} \sim \mathrm{Poisson}(\lambda_{1i}) and Y_{2i} \sim \mathrm{Poisson}(\lambda_{2i}), so that param1 describes the rate of the counts that are added and param2 the rate of the counts that are subtracted, as in dskellam(y, lambda1, lambda2). Both rates use a log link, \log\lambda_{ki} = x_i^\top\beta_k, with the same covariates (an intercept is added).

The log-likelihood is computed with dskellam(log = TRUE) and maximized by optim (BFGS) with the analytic gradient, starting deterministically from the moment estimates of an intercept-only model, so repeated calls give identical results. Standard errors come from the inverse of the observed information (optimHess). Observations with a missing response or covariate are removed with a warning. A warning is also given if the optimizer does not converge or the information matrix is not invertible (standard errors are then NA).

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: Wald z statistics (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) Bayesian modelling of football outcomes: using the Skellam's distribution for the goal difference. IMA Journal of Management Mathematics 20(2), 133-145.

Examples

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