[R] lmomco package and confidence limits?

J. R. M. Hosking JRMH001 at gmail.com
Tue Nov 17 18:36:25 CET 2009


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




More information about the R-help mailing list