skellam {VGAM} | R Documentation |
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)
Arguments
lmu1 , lmu2 |
Link functions for the |
imu1 , imu2 |
Optional initial values for the parameters.
See |
nsimEIM , parallel , zero |
See |
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.
Value
An object of class "vglmff"
(see vglmff-class
).
The object is used by modelling functions such as vglm
and vgam
.
Warning
This VGAM family function seems fragile and very sensitive to the initial values. Use very cautiously!!
Note
Numerical problems may occur for data if \mu_1
and/or
\mu_2
are large.
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
Examples
## 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)