Levy Distribution Family Function

Estimates the scale parameter of the Levy distribution by maximum likelihood estimation.

Usage

levy(location = 0, lscale = "loglink", iscale = NULL)

Arguments

location

Location parameter. Must have a known value. Called $a$ below.

lscale

Parameter link function for the (positive) scale parameter $b$. See Links for more choices.

iscale

Initial value for the $b$ parameter. By default, an initial value is chosen internally.

Details

The Levy distribution is one of three stable distributions whose density function has a tractable form. The formula for the density is $$f(y;b) = \sqrt{\frac{b}{2\pi}} \exp \left( \frac{-b}{2(y - a)} \right) / (y - a)^{3/2} $$ where $a<y<\infty$ and $b>0$. Note that if $a$ is very close to min(y) (where y is the response), then numerical problem will occur. The mean does not exist. The median is returned as the fitted values.

Value

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

References

Nolan, J. P. (2005). Stable Distributions: Models for Heavy Tailed Data.

Author

T. W. Yee

Examples

nn <- 1000; loc1 <- 0; loc2 <- 10
myscale <- 1  # log link ==> 0 is the answer
ldata <-
  data.frame(y1 = loc1 + myscale/rnorm(nn)^2,  # Levy(myscale, a)
             y2 = rlevy(nn, loc = loc2, scale = exp(+2)))
# Cf. Table 1.1 of Nolan for Levy(1,0)
with(ldata, sum(y1 > 1) / length(y1))  # Should be 0.6827
#> [1] 0.678
with(ldata, sum(y1 > 2) / length(y1))  # Should be 0.5205
#> [1] 0.521

fit1 <- vglm(y1 ~ 1, levy(location = loc1), ldata, trace = TRUE)
#> Iteration 1: loglikelihood = -3287.6921
#> Iteration 2: loglikelihood = -3251.2061
#> Iteration 3: loglikelihood = -3249.8083
#> Iteration 4: loglikelihood = -3249.8064
#> Iteration 5: loglikelihood = -3249.8064
coef(fit1, matrix = TRUE)
#>             loglink(scale)
#> (Intercept)  -0.0006090238
Coef(fit1)
#>     scale 
#> 0.9993912 
summary(fit1)
#> 
#> Call:
#> vglm(formula = y1 ~ 1, family = levy(location = loc1), data = ldata, 
#>     trace = TRUE)
#> 
#> Coefficients: 
#>              Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -0.000609   0.044721  -0.014    0.989
#> 
#> Name of linear predictor: loglink(scale) 
#> 
#> Log-likelihood: -3249.806 on 999 degrees of freedom
#> 
#> Number of Fisher scoring iterations: 5 
#> 
head(weights(fit1, type = "work"))
#>      [,1]
#> [1,]  0.5
#> [2,]  0.5
#> [3,]  0.5
#> [4,]  0.5
#> [5,]  0.5
#> [6,]  0.5

fit2 <- vglm(y2 ~ 1, levy(location = loc2), ldata, trace = TRUE)
#> Iteration 1: loglikelihood = -5427.7446
#> Iteration 2: loglikelihood = -5387.0866
#> Iteration 3: loglikelihood = -5385.3411
#> Iteration 4: loglikelihood = -5385.338
#> Iteration 5: loglikelihood = -5385.338
coef(fit2, matrix = TRUE)
#>             loglink(scale)
#> (Intercept)       2.017616
Coef(fit2)
#>    scale 
#> 7.520376 
c(median = with(ldata, median(y2)),
  fitted.median = head(fitted(fit2), 1))
#>        median fitted.median 
#>      26.96287      26.53061