VGAM (version 1.1-14)

skellam: Skellam Distribution Family Function

Description

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

Usage

skellam(lmu1 = "loglink", lmu2 = "loglink", imu1 = NULL, imu2 = NULL,
        nsimEIM = 100, parallel = FALSE, zero = NULL)

Value

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

and vgam.

Details

The Skellam distribution models the difference between two independent Poisson distributions (with means \(\mu_{j}\), say). It has density function $$f(y;\mu_1,\mu_2) = \left( \frac{ \mu_1 }{\mu_2} \right)^{y/2} \, \exp(-\mu_1-\mu_2 ) \, I_{|y|}( 2 \sqrt{ \mu_1 \mu_2}) $$ where \(y\) is an integer, \(\mu_1 > 0\), \(\mu_2 > 0\). Here, \(I_v\) is the modified Bessel function of the first kind with order \(v\).

The mean is \(\mu_1 - \mu_2\) (returned as the fitted values), and the variance is \(\mu_1 + \mu_2\). Simulated Fisher scoring is implemented.

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, dpois, poissonff.

Examples

if (FALSE) {
sdata <- data.frame(x2 = runif(nn <- 1000))
sdata <- transform(sdata, mu1 = exp(1 + x2), mu2 = exp(1 + x2))
sdata <- transform(sdata, y = rskellam(nn, mu1, mu2))
fit1 <- vglm(y ~ x2, skellam, data = sdata, trace = TRUE, crit = "coef")
fit2 <- vglm(y ~ x2, skellam(parallel = TRUE), data = sdata, trace = TRUE)
coef(fit1, matrix = TRUE)
coef(fit2, matrix = TRUE)
summary(fit1)
# Likelihood ratio test for equal means:
pchisq(2 * (logLik(fit1) - logLik(fit2)),
       df = df.residual(fit2) - df.residual(fit1), lower.tail = FALSE)
lrtest(fit1, fit2)  # Alternative
}