This function generates random samples of specified size from a specified parent distribution. Subsequently, the type of parent distribution is fit to the L-moments of the generated sample. The fitted distribution is then plotted. It is the user's responsibility to have an active plot already drawn; unless the callplot option is TRUE. This function is useful to demonstration of sample size on the uncertainty of a fitted distribution---a motivation for this function is as a classroom exercise.
gen.freq.curves(n, para, F=NULL, nsim=10, callplot=TRUE, aslog=FALSE,
asprob=FALSE, showsample=FALSE, showparent=FALSE,
lowerCI=NA, upperCI=NA, FCI=NA, ...)This function is largely used for its graphical side effects, but if estimates of the lower and upper confidence limits are known (say from genci.simple) then this function can be used to evaluate the counts of simulations at nonexceedance probability FCI outside the limits provided in lowerCI and upperCI.
Sample size to draw from parent as specified by para.
The parameters from lmom2par or vec2par.
The nonexceedance probabilities for horizontal axis---defaults to nonexceeds when the argument is NULL.
The number of simulations to perform (frequency curves to draw)---the default is 10.
Calls plot to acquire a graphics device---default is TRUE, but the called plot is left empty.
Compute log10 of quantiles---note that
NaNs produced in: log(x, base)will be produced for less than zero values. Otherwise this is a harmless message.
The qnorm function is used to convert nonexceedance probabilities, which are produced by nonexceeds, to Standard Normal variates. The Normal distribution will be a straight line when this argument is TRUE and aslog=FALSE.
Each simulated sample is drawn through plotting positions (pp).
The curve for the parent distribution is plotted on exit from the function if TRUE. Further plotting options can not be controlled---unlike the situation with the drawing of the simulated frequency curves.
An optional estimate of the lower confidence limit for the FCI nonexceedance probability.
An optional estimate of the upper confidence limit for the FCI nonexceedance probability.
The nonexeedance probability of interest for the confidence limits provided in lowerCI and upperCI.
Additional parameters are passed to the lines call within the function---except for the drawing of the parent distribution (see argument showparent).
W.H. Asquith
lmom2par, nonexceeds, rlmomco, lmoms
if (FALSE) {
# 1-day rainfall Travis county, Texas
para <- vec2par(c(3.00, 1.20, -.0954), type="gev")
F <- .99 # the 100-year event
n <- 46 # sample size derived from 75th percentile of record length distribution
# for Edwards Plateau from Figure 3 of USGS WRIR98-4044 (Asquith, 1998)
# Argument for 75th percentile is that the contours of distribution parameters
# in that report represent a regionalization of the parameters and hence
# record lengths such as the median or smaller for the region seem too small
# for reasonable exploration of confidence limits of precipitation.
nsim <- 5000 # simulation size
seed <- runif(1, min=1, max=10000)
set.seed(seed)
CI <- genci.simple(para, n, F=F, nsim=nsim, edist="nor")
lo.nor <- CI$lower; hi.nor <- CI$upper
set.seed(seed)
CI <- genci.simple(para, n, F=F, nsim=nsim, edist="aep4")
lo.aep4 <- CI$lower; hi.aep4 <- CI$upper
message("NORMAL ERROR DIST: lowerCI = ",lo.nor, " and upperCI = ",hi.nor)
message(" AEP4 ERROR DIST: lowerCI = ",lo.aep4," and upperCI = ",hi.aep4)
qF <- qnorm(F)
# simulated are grey, parent is black
set.seed(seed)
counts.nor <- gen.freq.curves(n, para, nsim=nsim,
asprob=TRUE, showparent=TRUE, col=rgb(0,0,1,0.025),
lowerCI=lo.nor, upperCI=hi.nor, FCI=F)
set.seed(seed)
counts.aep4 <- gen.freq.curves(n, para, nsim=nsim,
asprob=TRUE, showparent=TRUE, col=rgb(0,0,1,0.025),
lowerCI=lo.aep4, upperCI=hi.aep4, FCI=F)
lines( c(qF,qF), c(lo.nor, hi.nor), lwd=2, col=2)
points(c(qF,qF), c(lo.nor, hi.nor), pch=1, lwd=2, col=2)
lines( c(qF,qF), c(lo.aep4,hi.aep4), lwd=2, col=2)
points(c(qF,qF), c(lo.aep4,hi.aep4), pch=2, lwd=2, col=2)
percent.nor <- (counts.nor$count.above.upperCI +
counts.nor$count.below.lowerCI) /
counts.nor$count.valid.simulations
percent.aep4 <- (counts.aep4$count.above.upperCI +
counts.aep4$count.below.lowerCI) /
counts.aep4$count.valid.simulations
percent.nor <- 100 * percent.nor
percent.aep4 <- 100 * percent.aep4
message("NORMAL ERROR DIST: ",percent.nor)
message(" AEP4 ERROR DIST: ",percent.aep4)
# Continuing on, we are strictly focused on F being equal to 0.99
# Also we are no restricted to the example using the GEV distribution
# The vargev() function is from Handbook of Hydrology
"vargev" <-
function(para, n, F=c("F080", "F090", "F095", "F099", "F998", "F999")) {
F <- as.character(F)
if(! are.pargev.valid(para)) return()
F <- match.arg(F)
A <- para$para[2]
K <- para$para[3]
AS <- list(F080=c(-1.813, 3.017, -1.4010, 0.854),
F090=c(-2.667, 4.491, -2.2070, 1.802),
F095=c(-3.222, 5.732, -2.3670, 2.512),
F098=c(-3.756, 7.185, -2.3140, 4.075),
F099=c(-4.147, 8.216, -0.2033, 4.780),
F998=c(-5.336, 10.711, -1.1930, 5.300),
F999=c(-5.943, 11.815, -0.6300, 6.262))
AS <- as.environment(AS); CO <- get(F, AS)
varx <- A^2 * exp( CO[1] + CO[2]*exp(-K) + CO[3]*K^2 + CO[4]*K^3 ) / n
names(varx) <- NULL
return(varx)
}
sdx <- sqrt(vargev(para, n, F="F099"))
VAL <- qlmomco(F, para)
lo.vargev <- VAL + qt(0.05, df=n) * sdx # minus covered by return of qt()
hi.vargev <- VAL + qt(0.95, df=n) * sdx
set.seed(seed)
counts.vargev <- gen.freq.curves(n, para, nsim=nsim,
xlim=c(0,3), ylim=c(3,15),
asprob=TRUE, showparent=TRUE, col=rgb(0,0,1,0.01),
lowerCI=lo.vargev, upperCI=hi.vargev, FCI=F)
percent.vargev <- (counts.vargev$count.above.upperCI +
counts.vargev$count.below.lowerCI) /
counts.vargev$count.valid.simulations
percent.vargev <- 100 * percent.vargev
lines(c(qF,qF), range(c(lo.nor, hi.nor,
lo.aep4, hi.aep4,
lo.vargev,hi.vargev)), col=2)
points(c(qF,qF), c(lo.nor, hi.nor), pch=1, lwd=2, col=2)
points(c(qF,qF), c(lo.aep4, hi.aep4), pch=3, lwd=2, col=2)
points(c(qF,qF), c(lo.vargev,hi.vargev), pch=2, lwd=2, col=2)
message("NORMAL ERROR DIST: ",percent.nor)
message(" AEP4 ERROR DIST: ",percent.aep4)
message("VARGEV ERROR DIST: ",percent.vargev)
}
Run the code above in your browser using DataLab