skellam (version 0.2.3)

Skellam: The Skellam Distribution

Description

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

Usage

dskellam(x, lambda1, lambda2 = lambda1, log = FALSE)
pskellam(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)
dskellam.sp(x, lambda1, lambda2 = lambda1, log = FALSE)
pskellam.sp(q, lambda1, lambda2 = lambda1, 
    lower.tail = TRUE, log.p = FALSE)

Value

dskellam gives the (log) density, pskellam gives the (log) distribution function,

qskellam gives the quantile function, and rskellam generates random deviates. Invalid lambdas will result in return value NaN, with a warning.

Arguments

x, q

vector of quantiles.

p

vector of probabilities.

n

number of observations. If length(n) > 1, the length is taken to be the number required.

lambda1, lambda2

vectors of (non-negative) means.

log, log.p

logical; if TRUE, probabilities p are given as log(p).

lower.tail

logical; if TRUE (default), probabilities are \(P[X \le x]\), otherwise, \(P[X > x]\).

Author

Jerry W. Lewis, Patrick E. Brown

Details

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}\) and \(\lambda_2 = \mu_2 - \rho \sqrt{\mu_1 \mu_2}\).

dskellam returns a value equivalent to $$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. The \(|x|\) differs from most Skellam expressions in the literature, but is correct since \(x\) is an integer, resulting in improved portability and (in R versions prior to 2.9) improved accuracy for \(x<0\). Exponential scaling is used in besselI to postpone numerical problems. When numerical problems do occur, a saddlepoint approximation is substituted, which typically gives at least 4-figure accuracy. An alternative representation is \(dchisq(2 \lambda_1,2(x+1),2\lambda_2) 2\) for \(x \ge 0\), and \(dchisq(2 \lambda_2,2(1-x),2\lambda_1) 2\) for \(x \le 0\); but in R besselI appears to be more accurately implemented (for very small probabilities) than dchisq.

pskellam(x,lambda1,lambda2) returns pchisq(2*lambda2, -2*x, 2*lambda1) for \(x<=0\) and 1 - pchisq(2*lambda1, 2*(x+1), 2*lambda2) for \(x>=0\). When pchisq incorrectly returns 0, a saddlepoint approximation is substituted, which typically gives at least 2-figure accuracy.

The quantile is defined as the smallest value \(x\) such that \(F(x) \ge p\), where \(F\) is the distribution function. For lower.tail=FALSE, the quantile is defined as the largest value \(x\) such that \(F(x,\)lower.tail=FALSE\() \le p\).

rskellam is calculated as rpois(n,lambda1)-rpois(n,lambda2)

dskellam.sp and pskellam.sp return saddlepoint approximations to the pmf and cdf. They are called by dskellam and pskellam when results from primary methods are in doubt.

References

Butler, R. (2007) Saddlepoint Approximations with Applications, Cambridge University Press, Cambridge & New York, p.17.

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, New York, pp.190-192.

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, 26.

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. Skellam distribution https://en.wikipedia.org/wiki/Skellam_distribution

Examples

require('skellam')

# one lambda = 0 ~ Poisson
c(dskellam(0:10,5,0), dpois(0:10,5))
c(dskellam(-(0:10),0,5), dpois(0:10,5))
c(pskellam(0:10,5,0,lower.tail=TRUE), 
    ppois(0:10,5,lower.tail=TRUE))
c(pskellam(0:10,5,0,lower.tail=FALSE),
    ppois(0:10,5,lower.tail=FALSE))
c(pskellam(-(0:10),0,5,lower.tail=FALSE),
    ppois(0:10-1,5,lower.tail=TRUE))
c(pskellam(-(0:10),0,5,lower.tail=TRUE),
    ppois(0:10-1,5,lower.tail=FALSE))

# both lambdas != 0 ~ convolution of Poissons
dskellam(1,0.5,0.75)  # sum(dpois(1+0:10,0.5)*dpois(0:10,0.75))
pskellam(1,0.5,0.75)  # sum(dskellam(-10:1,0.5,0.75))
dskellam(c(-1,1),c(12,10),c(10,12))  # c(0.0697968,0.0697968)
dskellam(c(-1,1),c(12,10),c(10,12),log=TRUE) 
 # log(dskellam(c(-1,1),c(12,10),c(10,12)))
dskellam(256,257,1)  
# 0.024829348733183769  
# exact result for comparison with saddlepoint
dskellam(-3724,2000,3000)  
# 3.1058145363400105e-308  
# exact result for comparison with saddlepoint 
# (still accurate in extreme tail)
pskellam(c(-1,0),c(12,10),c(10,12))  # c(0.2965079,0.7034921)
pskellam(c(-1,0),c(12,10),c(10,12),lower.tail=FALSE) 
 # 1-pskellam(c(-1,0),c(12,10),c(10,12))
pskellam(-2:2,8.5,10.25,log.p=TRUE)  # log(pskellam(-2:2,8.5,10.25))
qskellam(c(0.05,0.95),3,4)  # c(-5,3); pskellam(cbind(-6:-5,2:3),3,4)
qskellam(c(0.05,0.95),3,0)  # c(1,6); qpois(c(0.05,0.95),3)
rskellam(35,8.5,10.25)