sinmad: Singh-Maddala Distribution Family Function

Description Usage Arguments Details Value Note Author(s) References See Also Examples

View source: R/family.actuary.R

Description

Maximum likelihood estimation of the 3-parameter Singh-Maddala distribution.

Usage

1
2
3
4
5
sinmad(lscale = "loglink", lshape1.a = "loglink", lshape3.q = "loglink",
       iscale = NULL, ishape1.a = NULL, ishape3.q = NULL, imethod = 1,
       lss = TRUE, gscale = exp(-5:5), gshape1.a = exp(-5:5),
       gshape3.q = exp(-5:5), probs.y = c(0.25, 0.5, 0.75),
       zero = "shape")

Arguments

lss

See CommonVGAMffArguments for important information.

lshape1.a, lscale, lshape3.q

Parameter link functions applied to the (positive) parameters a, scale, and q. See Links for more choices.

iscale, ishape1.a, ishape3.q, imethod, zero

See CommonVGAMffArguments for information. For imethod = 2 a good initial value for ishape3.q is needed to obtain good estimates for the other parameters.

gscale, gshape1.a, gshape3.q

See CommonVGAMffArguments for information.

probs.y

See CommonVGAMffArguments for information.

Details

The 3-parameter Singh-Maddala distribution is the 4-parameter generalized beta II distribution with shape parameter p=1. It is known under various other names, such as the Burr XII (or just the Burr distribution), Pareto IV, beta-P, and generalized log-logistic distribution. More details can be found in Kleiber and Kotz (2003).

Some distributions which are special cases of the 3-parameter Singh-Maddala are the Lomax (a=1), Fisk (q=1), and paralogistic (a=q).

The Singh-Maddala distribution has density

f(y) = aq y^(a-1) / [b^a (1 + (y/b)^a)^(1+q)]

for a > 0, b > 0, q > 0, y >= 0. Here, b is the scale parameter scale, and the others are shape parameters. The cumulative distribution function is

F(y) = 1 - [1 + (y/b)^a]^(-q).

The mean is

E(Y) = b gamma(1 + 1/a) gamma(q - 1/a) / gamma(q)

provided -a < 1 < aq; these are returned as the fitted values. This family function handles multiple responses.

Value

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

Note

See the notes in genbetaII.

Author(s)

T. W. Yee

References

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

See Also

Sinmad, genbetaII, betaII, dagum, fisk, inv.lomax, lomax, paralogistic, inv.paralogistic, simulate.vlm.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
sdata <- data.frame(y = rsinmad(n = 1000, shape1 = exp(1),
                    scale = exp(2), shape3 = exp(0)))
fit <- vglm(y ~ 1, sinmad(lss = FALSE), data = sdata, trace = TRUE)
fit <- vglm(y ~ 1, sinmad(lss = FALSE, ishape1.a = exp(1)),
            data = sdata, trace = TRUE)
coef(fit, matrix = TRUE)
Coef(fit)
summary(fit)

# Harder problem (has the shape3.q parameter going to infinity)

set.seed(3)
sdata <- data.frame(y1 = rbeta(1000, 6, 6))
# hist(with(sdata, y1))
if (FALSE) {
# These struggle
  fit1 <- vglm(y1 ~ 1, sinmad(lss = FALSE), data = sdata, trace = TRUE)
  fit1 <- vglm(y1 ~ 1, sinmad(lss = FALSE), data = sdata, trace = TRUE,
               crit = "coef")
  Coef(fit1)
}
# Try this remedy:
fit2 <- vglm(y1 ~ 1, data = sdata, trace = TRUE, stepsize = 0.05, maxit = 99,
             sinmad(lss = FALSE, ishape3.q = 3, lshape3.q = "logloglink"))
             
coef(fit2, matrix = TRUE)
Coef(fit2)

Example output

Loading required package: stats4
Loading required package: splines
VGLM    linear loop  1 :  loglikelihood = -2944.9284
VGLM    linear loop  2 :  loglikelihood = -2944.9212
VGLM    linear loop  3 :  loglikelihood = -2944.9212
VGLM    linear loop  1 :  loglikelihood = -2944.9284
VGLM    linear loop  2 :  loglikelihood = -2944.9212
VGLM    linear loop  3 :  loglikelihood = -2944.9212
            loglink(shape1.a) loglink(scale) loglink(shape3.q)
(Intercept)          1.062291       1.953652       -0.05679876
 shape1.a     scale  shape3.q 
2.8929907 7.0544043 0.9447842 

Call:
vglm(formula = y ~ 1, family = sinmad(lss = FALSE, ishape1.a = exp(1)), 
    data = sdata, trace = TRUE)

Pearson residuals:
                     Min      1Q   Median     3Q    Max
loglink(shape1.a) -6.228 -0.3902  0.34219 0.7442 0.8747
loglink(scale)    -2.069 -0.9770 -0.05003 1.0063 1.4397
loglink(shape3.q) -9.256 -0.1508  0.06832 0.3409 4.4608

Coefficients: 
              Estimate Std. Error z value Pr(>|z|)    
(Intercept):1  1.06229    0.04872  21.803   <2e-16 ***
(Intercept):2  1.95365    0.06171  31.659   <2e-16 ***
(Intercept):3 -0.05680    0.11296  -0.503    0.615    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Names of linear predictors: loglink(shape1.a), loglink(scale), 
loglink(shape3.q)

Log-likelihood: -2944.921 on 2997 degrees of freedom

Number of Fisher scoring iterations: 3 

No Hauck-Donner effect found in any of the estimates

Taking a modified step...
VGLM    linear loop  2 :  loglikelihood = -1609.9604
Taking a modified step....
VGLM    linear loop  3 :  loglikelihood = -1585.3258
Taking a modified step.....
VGLM    linear loop  4 :  loglikelihood = -1576.4835
Taking a modified step......
VGLM    linear loop  5 :  loglikelihood = -1575.826
Taking a modified step.........
VGLM    linear loop  6 :  loglikelihood = -1575.2006
Taking a modified step..........
VGLM    linear loop  7 :  loglikelihood = -1575.089
Taking a modified step............
VGLM    linear loop  8 :  loglikelihood = -1575.0724
Taking a modified step...............
VGLM    linear loop  9 :  loglikelihood = -1575.061
Taking a modified step................
VGLM    linear loop  10 :  loglikelihood = -1575.0576
Taking a modified step.................
Warning message:
In vglm.fitter(x = x, y = y, w = w, offset = offset, Xm2 = Xm2,  :
  iterations terminated because half-step sizes are very small
            loglink(shape1.a) loglink(scale) logloglink(shape3.q)
(Intercept)          1.541345      -1.538693            -16.90979
 shape1.a     scale  shape3.q 
4.6708680 0.2146615 1.0000000 

VGAM documentation built on Jan. 16, 2021, 5:21 p.m.