Benini Distribution Family Function

Estimating the 1-parameter Benini distribution by maximum likelihood estimation.

Usage

benini1(y0 = stop("argument 'y0' must be specified"),
        lshape = "loglink", ishape = NULL, imethod = 1,
        zero = NULL, parallel = FALSE,
        type.fitted = c("percentiles", "Qlink"),
        percentiles = 50)

Arguments

y0

Positive scale parameter.

lshape

Parameter link function and extra argument of the parameter $b$, which is the shape parameter. See Links for more choices. A log link is the default because $b$ is positive.

ishape

Optional initial value for the shape parameter. The default is to compute the value internally.

imethod, zero, parallel

Details at CommonVGAMffArguments.

type.fitted, percentiles

See CommonVGAMffArguments for information. Using "Qlink" is for quantile-links in VGAMextra.

Details

The Benini distribution has a probability density function that can be written $$f(y) = 2 s \exp(-s[(\log(y/y_0))^2]) \log(y/y_0) / y $$ for $0 < y_0 < y$, and shape $s > 0$. The cumulative distribution function for $Y$ is $$F(y) = 1 - \exp(-s[(\log(y/y_0))^2]).$$ Here, Newton-Raphson and Fisher scoring coincide. The median of $Y$ is now returned as the fitted values, by default. This VGAM family function can handle a multiple responses, which is inputted as a matrix.

On fitting, the extra slot has a component called y0 which contains the value of the y0 argument.

Value

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

References

Kleiber, C. and Kotz, S. (2003). Statistical Size Distributions in Economics and Actuarial Sciences, Hoboken, NJ, USA: Wiley-Interscience.

Author

T. W. Yee

Note

Yet to do: the 2-parameter Benini distribution estimates another shape parameter $a$ too. Hence, the code may change in the future.

Examples

y0 <- 1; nn <- 3000
bdata <- data.frame(y  = rbenini(nn, y0 = y0, shape = exp(2)))
fit <- vglm(y ~ 1, benini1(y0 = y0), data = bdata, trace = TRUE)
#> Iteration 1: loglikelihood = 179.036359
#> Iteration 2: loglikelihood = 179.036359
coef(fit, matrix = TRUE)
#>             loglink(shape)
#> (Intercept)       1.976216
Coef(fit)
#>    shape 
#> 7.215386 
fit@extra$y0
#> [1] 1
c(head(fitted(fit), 1), with(bdata, median(y)))  # Should be equal
#> [1] 1.363348 1.370697