levy.RdEstimates the scale parameter of the Levy distribution by maximum likelihood estimation.
levy(location = 0, lscale = "loglink", iscale = NULL)Location parameter. Must have a known value. Called \(a\) below.
Parameter link function for the (positive) scale parameter
\(b\). See Links for more choices.
Initial value for the \(b\) parameter. By default, an initial value is chosen internally.
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.
An object of class "vglmff"
(see vglmff-class).
The object is used by modelling functions
such as vglm,
and vgam.
Nolan, J. P. (2005). Stable Distributions: Models for Heavy Tailed Data.
The Nolan article was at
http://academic2.american.edu/~jpnolan/stable/chap1.pdf.
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
#>
#> No Hauck-Donner effect found in any of the estimates
#>
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