Douglas M. Hultstrand wrote:
Hello,
I am using the lmomco package (lmom.ub and pargev) to compute the GEV
parameters (location, scale, and shape), which are used to estimate
return values. I was wondering how/if I can calculate upper and lower
confidence (CI_u, CI_l) intervals for each return frequency using the
GEV parameters to fill-in the table below?
Xi (location) = 35.396
Alpha (scale) = 1.726
Kappa (shape) = 0.397
Year, Return value, CI_u, CL_l
2, 36.0,
5, 37.3,
10, 38.0,
25, 38.5,
50, 38.8,
100, 39.0,
Thank you for your help/suggestions,
Doug
As far as I am aware, neither lmomco nor any of the other L-moment
packages (lmom, Lmoments, ...) provide confidence interval calculations
directly.
An algebraic expression for the asymptotic covariance matrix
of the parameter estimates is given by Hosking, Wallis and Wood
(Technometrics, vol.27 (1985), pp.251-261, eq. (16)), and with a bit
more algebra this can be made to yield asymptotic standard errors
of estimated quantiles and thence normal-theory confidence intervals
for the quantiles.
Alternatively you could use some kind of bootstrap approach. Perhaps
the simplest method, assuming that your estimates were obtained from a
single sample of data, is to treat the sample as a "region" containing
a single site and to use the methods of regional frequency analysis in
package lmomRFA. This is roughly equivalent to a parametric bootstrap.
An example follows.
library(lmomRFA)
# A data sample
set.seed(1234)
zz <- quagev(runif(30), c(35.396,1.726,0.397))
# Compute L-moments of the sample, considered as a 1-site region
rdat <- regsamlmu(zz)
# Fit GEV distribution to the regional L-moments
rfit <- regfit(rdat,"gev")
# Generate simulations of an artificial 1-site region whose
# frequency distribution is the one fitted to the actual data
sim <- regsimq(rfit$qfunc, nrec = rdat$n,
f = 1 - 1 / c(2,5,10,25,50,100))
# Compute error bounds for quantiles of the site's
# frequency distribution
sitequantbounds(sim, rfit)
J. R. M. Hosking
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
