Estimate Parameters of a Generalized Extreme Value Distribution
Estimate the location, scale and shape parameters of a generalized extreme value distribution, and optionally construct a confidence interval for one of the parameters.
egevd(x, method = "mle", pwme.method = "unbiased", tsoe.method = "med", plot.pos.cons = c(a = 0.35, b = 0), ci = FALSE, ci.parameter = "location", ci.type = "two-sided", ci.method = "normal.approx", information = "observed", conf.level = 0.95)
x |
numeric vector of observations. |
method |
character string specifying the method of estimation. Possible values are
|
pwme.method |
character string specifying what method to use to compute the
probability-weighted moments when |
tsoe.method |
character string specifying the robust function to apply in the second stage of
the two-stage order-statistics estimator when |
plot.pos.cons |
numeric vector of length 2 specifying the constants used in the formula for the
plotting positions when |
ci |
logical scalar indicating whether to compute a confidence interval for the
location, scale, or shape parameter. The default value is |
ci.parameter |
character string indicating the parameter for which the confidence interval is
desired. The possible values are |
ci.type |
character string indicating what kind of confidence interval to compute. The
possible values are |
ci.method |
character string indicating what method to use to construct the confidence interval
for the location or scale parameter. Currently, the only possible value is
|
information |
character string indicating which kind of Fisher information to use when
computing the variance-covariance matrix of the maximum likelihood estimators.
The possible values are |
conf.level |
a scalar between 0 and 1 indicating the confidence level of the confidence interval.
The default value is |
If x
contains any missing (NA
), undefined (NaN
) or
infinite (Inf
, -Inf
) values, they will be removed prior to
performing the estimation.
Let \underline{x} = (x_1, x_2, …, x_n) be a vector of
n observations from a generalized extreme value distribution with
parameters location=
η, scale=
θ, and
shape=
κ.
Estimation
Maximum Likelihood Estimation (method="mle"
)
The log likelihood function is given by:
L(η, θ, κ) = -n \, log(θ) - (1 - κ) ∑^n_{i=1} y_i - ∑^n_{i=1} e^{y_i}
where
y_i = -\frac{1}{κ} log[\frac{1 - κ(x_i - η)}{θ}]
(see, for example, Jenkinson, 1969; Prescott and Walden, 1980; Prescott and Walden, 1983; Hosking, 1985; MacLeod, 1989). The maximum likelihood estimators (MLE's) of η, θ, and κ are those values that maximize the likelihood function, subject to the following constraints:
θ > 0
κ ≤ 1
x_i < η + \frac{θ}{κ} \; if κ > 0
x_i > η + \frac{θ}{κ} \; if κ < 0
Although in theory the value of κ may lie anywhere in the interval (-∞, ∞) (see GEVD), the constraint κ ≤ 1 is imposed because when κ > 1 the likelihood can be made infinite and thus the MLE does not exist (Castillo and Hadi, 1994). Hence, this method of estimation is not valid when the true value of κ is larger than 1. Hosking (1985) and Hosking et al. (1985) note that in practice the value of κ tends to lie in the interval -1/2 < κ < 1/2.
The value of -L is minimized using the R function nlminb
.
Prescott and Walden (1983) give formulas for the gradient and Hessian. Only
the gradient is supplied in the call to nlminb
. The values of
the PWME (see below) are used as the starting values. If the starting value of
κ is less than 0.001 in absolute value, it is reset to
sign(k) * 0.001
, as suggested by Hosking (1985).
Probability-Weighted Moments Estimation (method="pwme"
)
The idea of probability-weighted moments was introduced by Greenwood et al. (1979).
Landwehr et al. (1979) derived probability-weighted moment estimators (PWME's) for
the parameters of the Type I (Gumbel) extreme value distribution.
Hosking et al. (1985) extended these results to the generalized extreme value
distribution. See the abstract for Hosking et al. (1985)
for details on how these estimators are computed.
Two-Stage Order Statistics Estimation (method="tsoe"
)
The two-stage order statistics estimator (TSOE) was introduced by
Castillo and Hadi (1994) as an alternative to the MLE and PWME. Unlike the
MLE and PWME, the TSOE of κ exists for all combinations of sample
values and possible values of κ. See the
abstract for Castillo and Hadi (1994) for details
on how these estimators are computed. In the second stage,
Castillo and Hadi (1984) suggest using either the median or the least median of
squares as the robust function. The function egevd
allows three options
for the robust function: median (tsoe.method="med"
; see the R help file for
median
), least median of squares (tsoe.method="lms"
;
see the help file for lmsreg
in the package MASS),
and least trimmed squares (tsoe.method="lts"
; see the help file for
ltsreg
in the package MASS).
Confidence Intervals
When ci=TRUE
, an approximate (1-α)100% confidence intervals
for η can be constructed assuming the distribution of the estimator of
η is approximately normally distributed. A two-sided confidence
interval is constructed as:
[\hat{η} - t(n-1, 1-α/2) \hat{σ}_{\hat{η}}, \, \hat{η} + t(n-1, 1-α/2) \hat{σ}_{\hat{η}}]
where t(ν, p) is the p'th quantile of Student's t-distribution with ν degrees of freedom, and the quantity
\hat{σ}_{\hat{η}}
denotes the estimated asymptotic standard deviation of the estimator of η.
Similarly, a two-sided confidence interval for θ is constructed as:
[\hat{θ} - t(n-1, 1-α/2) \hat{σ}_{\hat{θ}}, \, \hat{θ} + t(n-1, 1-α/2) \hat{σ}_{\hat{θ}}]
and a two-sided confidence interval for κ is constructed as:
[\hat{κ} - t(n-1, 1-α/2) \hat{σ}_{\hat{κ}}, \, \hat{κ} + t(n-1, 1-α/2) \hat{σ}_{\hat{κ}}]
One-sided confidence intervals for η, θ, and κ are computed in a similar fashion.
Maximum Likelihood Estimator (method="mle"
)
Prescott and Walden (1980) derive the elements of the Fisher information matrix
(the expected information). The inverse of this matrix, evaluated at the values
of the MLE, is the estimated asymptotic variance-covariance matrix of the MLE.
This method is used to estimate the standard deviations of the estimated
distribution parameters when information="expected"
. The necessary
regularity conditions hold for κ < 1/2. Thus, this method of
constructing confidence intervals is not valid when the true value of
κ is greater than or equal to 1/2.
Prescott and Walden (1983) derive expressions for the observed information matrix
(i.e., the Hessian). This matrix is used to compute the estimated asymptotic
variance-covariance matrix of the MLE when information="observed"
.
In computer simulations, Prescott and Walden (1983) found that the variance-covariance matrix based on the observed information gave slightly more accurate estimates of the variance of MLE of κ compared to the estimated variance based on the expected information.
Probability-Weighted Moments Estimator (method="pwme"
)
Hosking et al. (1985) show that these estimators are asymptotically multivariate
normal and derive the asymptotic variance-covariance matrix. See the
abstract for Hosking et al. (1985) for details on how
this matrix is computed.
Two-Stage Order Statistics Estimator (method="tsoe"
)
Currently there is no built-in method in EnvStats for computing confidence
intervals when method="tsoe"
. Castillo and Hadi (1994) suggest
using the bootstrap or jackknife method.
a list of class "estimate"
containing the estimated parameters and other information.
See estimate.object
for details.
Two-parameter extreme value distributions (EVD) have been applied extensively since the 1930's to several fields of study, including the distributions of hydrological and meteorological variables, human lifetimes, and strength of materials. The three-parameter generalized extreme value distribution (GEVD) was introduced by Jenkinson (1955) to model annual maximum and minimum values of meteorological events. Since then, it has been used extensively in the hydological and meteorological fields.
The three families of EVDs are all special kinds of GEVDs. When the shape
parameter κ=0, the GEVD reduces to the Type I extreme value (Gumbel)
distribution. (The function zTestGevdShape
allows you to test
the null hypothesis H_0: κ=0.) When κ > 0, the GEVD is
the same as the Type II extreme value distribution, and when κ < 0
it is the same as the Type III extreme value distribution.
Hosking et al. (1985) compare the asymptotic and small-sample statistical properties of the PWME with the MLE and Jenkinson's (1969) method of sextiles. Castillo and Hadi (1994) compare the small-sample statistical properties of the MLE, PWME, and TSOE. Hosking and Wallis (1995) compare the small-sample properties of unbaised L-moment estimators vs. plotting-position L-moment estimators. (PWMEs can be written as linear combinations of L-moments and thus have equivalent statistical properties.) Hosking and Wallis (1995) conclude that unbiased estimators should be used for almost all applications.
Steven P. Millard (EnvStats@ProbStatInfo.com)
Castillo, E., and A. Hadi. (1994). Parameter and Quantile Estimation for the Generalized Extreme-Value Distribution. Environmetrics 5, 417–432.
Forbes, C., M. Evans, N. Hastings, and B. Peacock. (2011). Statistical Distributions. Fourth Edition. John Wiley and Sons, Hoboken, NJ.
Greenwood, J.A., J.M. Landwehr, N.C. Matalas, and J.R. Wallis. (1979). Probability Weighted Moments: Definition and Relation to Parameters of Several Distributions Expressible in Inverse Form. Water Resources Research 15(5), 1049–1054.
Hosking, J.R.M. (1984). Testing Whether the Shape Parameter is Zero in the Generalized Extreme-Value Distribution. Biometrika 71(2), 367–374.
Hosking, J.R.M. (1985). Algorithm AS 215: Maximum-Likelihood Estimation of the Parameters of the Generalized Extreme-Value Distribution. Applied Statistics 34(3), 301–310.
Hosking, J.R.M., J.R. Wallis, and E.F. Wood. (1985). Estimation of the Generalized Extreme-Value Distribution by the Method of Probability-Weighted Moments. Technometrics 27(3), 251–261.
Jenkinson, A.F. (1969). Statistics of Extremes. Technical Note 98, World Meteorological Office, Geneva.
Johnson, N. L., S. Kotz, and N. Balakrishnan. (1995). Continuous Univariate Distributions, Volume 2. Second Edition. John Wiley and Sons, New York.
Landwehr, J.M., N.C. Matalas, and J.R. Wallis. (1979). Probability Weighted Moments Compared With Some Traditional Techniques in Estimating Gumbel Parameters and Quantiles. Water Resources Research 15(5), 1055–1064.
Macleod, A.J. (1989). Remark AS R76: A Remark on Algorithm AS 215: Maximum Likelihood Estimation of the Parameters of the Generalized Extreme-Value Distribution. Applied Statistics 38(1), 198–199.
Prescott, P., and A.T. Walden. (1980). Maximum Likelihood Estimation of the Parameters of the Generalized Extreme-Value Distribution. Biometrika 67(3), 723–724.
Prescott, P., and A.T. Walden. (1983). Maximum Likelihood Estimation of the Three-Parameter Generalized Extreme-Value Distribution from Censored Samples. Journal of Statistical Computing and Simulation 16, 241–250.
# Generate 20 observations from a generalized extreme value distribution # with parameters location=2, scale=1, and shape=0.2, then compute the # MLE and construct a 90% confidence interval for the location parameter. # (Note: the call to set.seed simply allows you to reproduce this example.) set.seed(498) dat <- rgevd(20, location = 2, scale = 1, shape = 0.2) egevd(dat, ci = TRUE, conf.level = 0.9) #Results of Distribution Parameter Estimation #-------------------------------------------- # #Assumed Distribution: Generalized Extreme Value # #Estimated Parameter(s): location = 1.6144631 # scale = 0.9867007 # shape = 0.2632493 # #Estimation Method: mle # #Data: dat # #Sample Size: 20 # #Confidence Interval for: location # #Confidence Interval Method: Normal Approximation # (t Distribution) based on # observed information # #Confidence Interval Type: two-sided # #Confidence Level: 90% # #Confidence Interval: LCL = 1.225249 # UCL = 2.003677 #---------- # Compare the values of the different types of estimators: egevd(dat, method = "mle")$parameters # location scale shape #1.6144631 0.9867007 0.2632493 egevd(dat, method = "pwme")$parameters # location scale shape #1.5785779 1.0187880 0.2257948 egevd(dat, method = "pwme", pwme.method = "plotting.position")$parameters # location scale shape #1.5509183 0.9804992 0.1657040 egevd(dat, method = "tsoe")$parameters # location scale shape #1.5372694 1.0876041 0.2927272 egevd(dat, method = "tsoe", tsoe.method = "lms")$parameters #location scale shape #1.519469 1.081149 0.284863 egevd(dat, method = "tsoe", tsoe.method = "lts")$parameters # location scale shape #1.4840198 1.0679549 0.2691914 #---------- # Clean up #--------- rm(dat)
Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.