VGAM (version 1.0-2)

gev: Generalized Extreme Value Distribution Family Function

Description

Maximum likelihood estimation of the 3-parameter generalized extreme value (GEV) distribution.

Usage

gev(llocation = "identitylink", lscale = "loge", lshape = logoff(offset = 0.5), percentiles = c(95, 99), ilocation = NULL, iscale = NULL, ishape = NULL, imethod = 1, gprobs.y = (1:9)/10, gscale.mux = exp((-5:5)/6), gshape = (-5:5) / 11 + 0.01, iprobs.y = NULL, tolshape0 = 0.001, type.fitted = c("percentiles", "mean"), zero = c("scale", "shape"))
gevff(llocation = "identitylink", lscale = "loge", lshape = logoff(offset = 0.5), percentiles = c(95, 99), ilocation = NULL, iscale = NULL, ishape = NULL, imethod = 1, gprobs.y = (1:9)/10, gscale.mux = exp((-5:5)/6), gshape = (-5:5) / 11 + 0.01, iprobs.y = NULL, tolshape0 = 0.001, type.fitted = c("percentiles", "mean"), zero = c("scale", "shape"))

Arguments

llocation, lscale, lshape

Parameter link functions for $mu$, $sigma$ and $xi$ respectively. See Links for more choices.

For the shape parameter, the default logoff link has an offset called $A$ below; and then the linear/additive predictor is $log(xi+A)$ which means that $xi > -A$. For technical reasons (see Details) it is a good idea for $A = 0.5$.

percentiles

Numeric vector of percentiles used for the fitted values. Values should be between 0 and 100. This argument is ignored if type.fitted = "mean".

type.fitted

See CommonVGAMffArguments for information. The default is to use the percentiles argument. If "mean" is chosen, then the mean $mu + sigma * (gamma(1-xi)-1)/xi$ is returned as the fitted values, and these are only defined for $xi<1$.< p="">

ilocation, iscale, ishape

Numeric. Initial value for the location parameter, $sigma$ and $xi$. A NULL means a value is computed internally. The argument ishape is more important than the other two. If a failure to converge occurs, or even to obtain initial values occurs, try assigning ishape some value (positive or negative; the sign can be very important). Also, in general, a larger value of iscale tends to be better than a smaller value.

imethod

Initialization method. Either the value 1 or 2. If both methods fail then try using ishape. See CommonVGAMffArguments for information.

gshape

Numeric vector. The values are used for a grid search for an initial value for $xi$. See CommonVGAMffArguments for information.

gprobs.y, gscale.mux, iprobs.y

Numeric vectors, used for the initial values. See CommonVGAMffArguments for information.

tolshape0

Passed into dgev when computing the log-likelihood.

zero

A specifying which linear/additive predictors are modelled as intercepts only. The values can be from the set {1,2,3} corresponding respectively to $mu$, $sigma$, $xi$. If zero = NULL then all linear/additive predictors are modelled as a linear combination of the explanatory variables. For many data sets having zero = 3 is a good idea. See CommonVGAMffArguments for information.

Value

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

Warning

Currently, if an estimate of $xi$ is too close to 0 then an error may occur for gev() with multivariate responses. In general, gevff() is more reliable than gev(). Fitting the GEV by maximum likelihood estimation can be numerically fraught. If $1 + xi*(y-mu)/sigma <= 0$="" then="" some="" crude="" evasive="" action="" is="" taken="" but="" the="" estimation="" process="" can="" still="" fail.="" this="" particularly="" case="" if="" vgam with s is used; then smoothing is best done with vglm with regression splines (bs or ns) because vglm implements half-stepsizing whereas vgam doesn't (half-stepsizing helps handle the problem of straying outside the parameter space.)

Details

The GEV distribution function can be written $$G(y) = \exp( -[ (y-\mu)/ \sigma ]_{+}^{- 1/ \xi}) $$ where $sigma > 0$, $-Inf < mu < Inf$, and $1 + xi*(y-mu)/sigma > 0$. Here, $x_+ = max(x,0)$. The $mu$, $sigma$, $xi$ are known as the location, scale and shape parameters respectively. The cases $xi>0$, $xi<0$, $xi="0$" correspond="" to="" the="" frechet,="" weibull,="" and="" gumbel="" types="" respectively.="" it="" can="" be="" noted="" that="" (or="" type="" i)="" distribution="" accommodates="" many="" commonly-used="" distributions="" such="" as="" normal,="" lognormal,="" logistic,="" gamma,="" exponential="" weibull.<="" p="">

For the GEV distribution, the $k$th moment about the mean exists if $xi < 1/k$. Provided they exist, the mean and variance are given by $mu + sigma { Gamma(1-xi)-1} / xi$ and $sigma^2 { Gamma(1-2 xi) - Gamma^2 (1- xi) } / xi^2$ respectively, where $Gamma$ is the gamma function.

Smith (1985) established that when $xi > -0.5$, the maximum likelihood estimators are completely regular. To have some control over the estimated $xi$ try using lshape = logoff(offset = 0.5), say, or lshape = extlogit(min = -0.5, max = 0.5), say.

References

Yee, T. W. and Stephenson, A. G. (2007) Vector generalized linear and additive extreme value models. Extremes, 10, 1--19.

Tawn, J. A. (1988) An extreme-value theory model for dependent observations. Journal of Hydrology, 101, 227--250.

Prescott, P. and Walden, A. T. (1980) Maximum likelihood estimation of the parameters of the generalized extreme-value distribution. Biometrika, 67, 723--724.

Smith, R. L. (1985) Maximum likelihood estimation in a class of nonregular cases. Biometrika, 72, 67--90.

See Also

rgev, gumbel, gumbelff, guplot, rlplot.gevff, gpd, weibullR, frechet, extlogit, oxtemp, venice, CommonVGAMffArguments.

Examples

## Not run: 
# # Multivariate example
# fit1 <- vgam(cbind(r1, r2) ~ s(year, df = 3), gev(zero = 2:3),
#              data = venice, trace = TRUE)
# coef(fit1, matrix = TRUE)
# head(fitted(fit1))
# par(mfrow = c(1, 2), las = 1)
# plot(fit1, se = TRUE, lcol = "blue", scol = "forestgreen",
#      main = "Fitted mu(year) function (centered)", cex.main = 0.8)
# with(venice, matplot(year, depvar(fit1)[, 1:2], ylab = "Sea level (cm)",
#      col = 1:2, main = "Highest 2 annual sea levels", cex.main = 0.8))
# with(venice, lines(year, fitted(fit1)[,1], lty = "dashed", col = "blue"))
# legend("topleft", lty = "dashed", col = "blue", "Fitted 95 percentile")
# 
# # Univariate example
# (fit <- vglm(maxtemp ~ 1, gevff, data = oxtemp, trace = TRUE))
# head(fitted(fit))
# coef(fit, matrix = TRUE)
# Coef(fit)
# vcov(fit)
# vcov(fit, untransform = TRUE)
# sqrt(diag(vcov(fit)))  # Approximate standard errors
# rlplot(fit)
# ## End(Not run)