skellam {VGAM}R Documentation

Skellam Distribution Family Function

Description

Estimates the two parameters of a Skellam distribution by maximum likelihood estimation.

Usage

skellam(lmu1="loge", lmu2="loge",
        emu1=list(), emu2=list(),
        imu1=NULL, imu2=NULL,
        nsimEIM=100, parallel=FALSE, zero=NULL)

Arguments

lmu1, emu1 Link function and extra argument for the mu1 parameter. See Links for more choices and for general information.
lmu2, emu2 Link function and extra argument for the mu1 parameter. See Links for more choices and for general information.
imu1, imu2 Optional initial values for the parameters. See CommonVGAMffArguments for more information. If convergence failure occurs (this VGAM family function seems to require good initial values) try using these arguments.
nsimEIM, parallel, zero See CommonVGAMffArguments for more information. In particular, setting parallel=TRUE will constrain the two means to be equal.

Details

The Skellam distribution models the difference between two independent Poisson distributions. It has density function

f(y;mu1,mu2) = ( μ1 / mu_2 )^(y/2) * exp(-mu1-mu2 ) * I_y( 2 * sqrt(mu1*mu2))

where y is an integer, mu1 > 0, mu2 > 0. Here, I_v is the modified Bessel function of the first kind with order v.

The mean is mu1 - mu2 (returned as the fitted values) and the variance is mu1 + mu2. Simulated Fisher scoring is implemented.

Value

An object of class "vglmff" (see vglmff-class). The object is used by modelling functions such as vglm and vgam.

Note

Numerical problems may occur for data if mu1 and/or mu2 are large.

Author(s)

T. W. Yee

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

See Also

dskellam, poissonff.

Examples

x = runif(n <- 1000)
mu1 = exp(1+x); mu2 = exp(1+x);
y = rskellam(n, mu1, mu2)
fit1 = vglm(y ~ x, skellam, trace=TRUE, crit="l")
fit2 = vglm(y ~ x, skellam(parallel=TRUE), trace=TRUE, crit="c")
coef(fit1, matrix=TRUE)
coef(fit2, matrix=TRUE)
summary(fit1)
# Likelihood ratio test for equal means:
1-pchisq(logLik(fit1)-logLik(fit2), df=fit2@df.residual-fit1@df.residual)

[Package VGAM version 0.7-7 Index]