skellam (version 0.1.3)

Skellam: The Skellam Distribution

Description

Density, distribution function, quantile function and random number generation for the Skellam distribution with parameters lambda1 and lambda2. If $Y1$ and $Y2$ are Poisson variables with means $mu1$ and $mu2$ and correlation $rho$, then $X = Y1 - Y2$ is Skellam with parameters $lambda1 = mu1 - rho*sqrt(mu1*mu2)$ and $lambda2 = mu2 - rho*sqrt(mu1*mu2)$.

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)

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]$.

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.

source

The relation of dgamma to the modified Bessel function of the first kind was given by Skellam (1946). The relation of pgamma to the noncentral chi-square was given by Johnson (1959). Tables are given by Strackee and van der Gon (1962), which can be used to verify this implementation (cf. direct calculation in the examples below). qskellam uses the Cornish--Fisher expansion to include skewness and kurtosis corrections to a normal approximation, followed by a search. If getOption("verbose")==TRUE, then qskellam will not use qpois when one of the lambdas is zero, in order to verify that this search algorithm has been implemented properly.

Details

dskellam returns a value equivalent to $$I(2 \sqrt(lambda1 lambda2), abs(x)) (lambda1/lambda2)^(x/2) exp(-lambda1-lambda2)$$ where $I(y,nu)$ is the modified Bessel function of the first kind. The abs(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="" problems="" do="" occur,="" a="" saddlepoint="" approximation="" substituted,="" which="" typically="" gives="" at="" least="" 4-figure="" accuracy.="" an="" alternative="" representation="" $dchisq(2*lambda1,2*(x+1),2*lambda2)*2$="" for="" $x="" \ge="" 0$,="" and="" $dchisq(2*lambda2,2*(1-x),2*lambda1)*2$="" \le="" 0$;="" but="" 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 connexion 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 http://en.wikipedia.org/wiki/Skellam_distribution

Examples

require('skellam')

# one lambda = 0 ~ Poisson
dskellam(0:10,5,0)    # dpois(0:10,5)
dskellam(-(0:10),0,5) # dpois(0:10,5)
pskellam(0:10,5,0,lower.tail=TRUE)      # ppois(0:10,5,lower.tail=TRUE)
pskellam(0:10,5,0,lower.tail=FALSE)     # ppois(0:10,5,lower.tail=FALSE)
pskellam(-(0:10),0,5,lower.tail=FALSE)  # ppois(0:10-1,5,lower.tail=TRUE)
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)