Hello,
I am using the lmom and lmomRFA to compute the return frequencies using
the GEV distribution.Iam trying to generate upper and lower bound
frequency estimates.
I provided a working example of the code that I am using to estimate the
upper and lower bounds. Specific questions I have are:
1) Is the methodology appropriate and applied correctly?
2) In this example, are the simulated qhat values better than the actual
fitted data estimates (GEV)?
_*Code Example Below:*_
library(lmom); library(lmomRFA)
set.seed(1234)
# Example data
max_data -
c(11.14,10.95,10.21,9.88,9.85,9.74,9.74,9.73,9.62,8.95,8.38,8.20)
moments = samlmu(max_data, sort.data = TRUE)
parGEV - pelgev(moments) # GEV
x - c(10,25,50,100,500,1000)
Fx - (1-(1/x))
GEV = Fx
PDS - 1/(log(x)-log(x-1))
for (i in seq_along(Fx)) {
GEV[i] = round(quagev(Fx[i], parGEV), 3)
}
BOUNDS ###
# A data sample
zz_gev - quagev(runif(500), parGEV)
# Compute L-moments of the sample, considered as a 1-site region
rdat_gev - regsamlmu(zz_gev)
# Fit GEV distribution to the regional L-moments
rfit_gev - regfit(rdat_gev,gev)
# Generate simulations of an artificial 1-site region with frequency
distribution fitted to the actual data
sim_gev - regsimq(rfit_gev$qfunc, nrec = rdat_gev$n, f = 1 - 1 / x )
# Compute error bounds for quantiles of the site's frequency distribution
CI_gev - sitequantbounds(sim_gev, rfit_gev)
# 90% Error Bounds/CI
Clwr - round(CI_gev$bound.0.05, 3)
Cupr - round(CI_gev$bound.0.95, 3)
qhat - round(CI_gev$Qhat, 3)
# print data
data.frame(Fx, GEV, Clwr, Cupr, qhat)
END BOUNDS ##
#Plot and visualize the data
plot(x, GEV, log=x, ylim=c(10.5,12.5), col=blue)
lines(PDS,GEV, col=blue)
lines(PDS,Clwr, lty=3)
lines(PDS,Cupr, lty=3)
lines(PDS,qhat, col=red, lty=5)
Thanks for your help,
Doug
--
-
Douglas M. Hultstrand, MS
Senior Hydrometeorologist
Applied Weather Associates
Monument, Colorado
voice: 970-682-2462
mobile: 720-771-5840
www.appliedweatherassociates.com
dhultstr...@appliedweatherassociates.com
-
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.