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`.

### 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 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`, `dpois`, `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(2*(logLik(fit1)-logLik(fit2)),
             df=fit2@df.residual-fit1@df.residual)

[Package VGAM version 0.7-9 Index]