Skellam {skellam} | R Documentation |
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)
Arguments
x , q |
vector of quantiles. |
p |
vector of probabilities. |
n |
number of observations. If |
lambda1 , lambda2 |
vectors of (non-negative) means. |
log , log.p |
logical; if |
lower.tail |
logical; if |
Details
If and
are Poisson variables with means
and
and correlation
, then
is Skellam with parameters
and
.
dskellam
returns a value equivalent to
where is the modified Bessel function of the first kind.
The
differs from most Skellam expressions in the literature, but is correct since
is an integer,
resulting in improved portability and (in R versions prior to 2.9) improved accuracy for
.
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
for
, and
for
; 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 and
1 - pchisq(2*lambda1, 2*(x+1), 2*lambda2)
for .
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 such that
, where
is the distribution function.
For
lower.tail=FALSE
, the quantile is defined as the largest value such that
lower.tail=FALSE
.
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.
Value
dskellam
gives the (log) density, pskellam
gives the (log) distribution function,
qskellam
gives the quantile function, and rskellam
generates random deviates.
Invalid lambda
s will result in return value NaN
, with a warning.
Note
The VGAM package
also contains Skellam functions, which are syntactically similar;
independently developed versions are included here for completeness.
Moreover, this dskellam
function offers a broader working range,
correct handling of cases where at least one rate parameter is zero,
enhanced argument checking,
and (in R versions prior to 2.9) improved accuracy for .
If both packages are loaded,
skellam::dskellam
or
VGAM::dskellam
can unambiguously specify which implementation to use.
Author(s)
Jerry W. Lewis, Patrick E. Brown
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 lambda
s is zero,
in order to verify that this search algorithm has been implemented properly.
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 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)