Skellam Distribution Family Function
Estimates the two parameters of a Skellam distribution by maximum likelihood estimation.
skellam(lmu1 = "loglink", lmu2 = "loglink", imu1 = NULL, imu2 = NULL,
nsimEIM = 100, parallel = FALSE, zero = NULL)lmu1, lmu2 |
Link functions for the mu1 and mu2 parameters.
See |
imu1, imu2 |
Optional initial values for the parameters.
See |
nsimEIM, parallel, zero |
See |
The Skellam distribution models the difference between two independent Poisson distributions (with means mu_j, say). 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.
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!!
Numerical problems may occur for data if mu1 and/or mu2 are large.
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.
## Not run:
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
## End(Not run)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.