VGAM (version 1.0-4)

skellam: Skellam Distribution Family Function


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


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


lmu1, lmu2

Link functions for the \(\mu_1\) and \(\mu_2\) parameters. 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 information. In particular, setting parallel=TRUE will constrain the two means to be equal.


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


This VGAM family function seems fragile and very sensitive to the initial values. Use very cautiously!!


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.


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.


Run this code
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)
# 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
# }

Run the code above in your browser using DataCamp Workspace