[R] lmomRFA: update "regfit" object

2017-05-04 Thread Douglas Hultstrand
Hello,

I am creating error bounds based on simulated data across multiple data 
durations.  I was wondering if there is a way to update an object class 
from "regfit" from lmomRFA package?  The reason is for consistency 
across durations.

Example below:

library(lmom); library(lmomRFA)

data <- c(0.42, 0.13, 0.59, 0.12, 0.78, 0.17, 0.3, 0.41, 0.28, 0.79)
# random data

reg <- regsamlmu(data)# calc l-moments

*org_gev <- regfit(reg,"gev")*# original gev fit

# UPDATE *"org_gev*" values (below) and save as "update_gev"

# xi = 0.65

# alpha = 0. 51

# k = -0.023


Thank you for the help,
Doug




-- 
-
Douglas M. Hultstrand, MS
Senior Hydrometeorologist
Applied Weather Associates
Monument, Colorado
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.


[R] lmom and lmomRFA Upper and Lower Bounds Simulation Questions

2015-05-10 Thread Douglas Hultstrand
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.