Description Usage Arguments Value Note Author(s) See Also Examples

This function estimates the lower and upper limits of a specified confidence interval for an aribitrary quantile value of a specified parent distribution [quantile function *Q(F,θ)* with parameters *θ*] using Monte Carlo simulation. The quantile value, actually the nonexceedance probability (*F* for *0 ≤ F ≤ 1*) of the value, is specified by the user. The user also provides the parameters of the parent distribution (see `lmom2par`

). This function does consider an estimate of the variance-covariance structure of the sample data (for that see `qua2ci.cov`

). The `qua2ci.simple`

is the original implementation and dates close to the initial releases of lmomco and was originally named `qua2ci`

. That name is now deprecated but retained as an alias, which will be removed at some later release.

For `nsim`

simulation runs (ideally a large number), samples of size *n* are drawn from *Q(F,θ)*. The L-moments of each simulated sample are computed using `lmoms`

and a distribution of the same type is fit. The *F*-quantile of the just-fitted distribution is computed and placed into a vector. The process of simulating the sample, computing the L-moments, computing the parameters, and solving for the *F*-quantile is repeated for the specified number of simulation runs.

To estimate the confidence interval, the L-moments of the vector simulated quantiles are computed. Subsequently, the parameters of a user-specified distribution “error” distribution (`edist`

) are computed. The two quantiles of this error distribution for the specified confidence interval are computed. These two quantiles represent the estimated lower and upper limits for the confidence interval of the parent distribution for samples of size *n*. The error distribution defaults to the Generalized Normal (see `pargno`

) because this distribution has the Normal as a special case but extends the fit to the 3rd L-moment (*τ_3*) for exotic situations in which some asymmetry in the quantile distribution might exist.

Finally, it is often useful to have vectors of lower and upper limits for confidence intervals for a vector of *F* values. The function `genci.simple`

does just that and uses `qua2ci.simple`

as the computational underpinning.

1 2 | ```
qua2ci.simple(f,para,n, level=0.90, edist="gno", nsim=1000, showpar=FALSE,
empdist=TRUE, verbose=FALSE, maxlogdiff=6, ...)
``` |

`f` |
Nonexceedance probability ( |

`para` |
The parameters from |

`n` |
The sample size for each Monte Carlo simulation will use. |

`level` |
The confidence interval ( |

`edist` |
The model for the error distribution. Although the Normal (the default) commonly is assumed in error analyses, it need not be, as support for other distributions supported by lmomco is available. The default is the Generalized Normal so the not only is the Normal possible but asymmetry is also accomodated ( |

`nsim` |
The number of simulations (replications) for the sample size |

`showpar` |
The parameters of the |

`empdist` |
If |

`verbose` |
The verbosity of the operation of the function. |

`maxlogdiff` |
The maximum permitted difference in log10 space between a simulated quantile and the true value. It is possible that a well fit simulated sample to the parent distribution type provides crazy quantile estimates in the far reaches of either tail. The default value of 6 was chosen based on experience with the Kappa distribution fit to a typical heavy-right tail flood magnitude data set. The concern motivating this feature is that as the number of parameters increases, it seems progressively there is more chance for a distribution tail to swing wildy into regions for which an analyst would not be comfortable with given discipline-specific knowledge. The choice of 6-log cycles is |

`...` |
Additional arguments to pass such as to |

An **R** `list`

is returned. The `lwr`

and `upr`

match the nomenclature of `qua2ci.cov`

but because `qua2ci.simple`

is provided the parent, the `true`

value is returned, whereas `qua2ci.cov`

returns the `fit`

.

`lwr` |
The lower value of the confidence interval having nonexceedance probability equal to |

`true` |
The value returned by |

`upr` |
The upper value of the confidence interval having nonexceedance probability equal to |

`elmoms` |
The L-moments from |

`epara` |
The parameters of the error distribution fit using the |

`empdist` |
An |

`ifail` |
A diagnostic value. A value of zero means that successful exit was made. |

`ifailtext` |
A descriptive message related to the |

`nsim` |
An echoing of the |

`sim.attempts` |
The number of executions of the |

The `empdist`

element in the returned `list`

is an **R** `environment`

that contains:

`simquas` |
A |

`empir.dist.lwr` |
The |

`empir.dist.upr` |
The |

`bern.smooth.lwr` |
The |

`bern.smooth.upr` |
The |

`epmoms` |
The product moments of the simulated quantiles from |

This function relies on a `while`

loop that runs until `nsim`

have successfully completed. Some reasons for an early `next`

in the loop include invalid L-moments by `are.lmom.valid`

of the simluated data or invalid fitted parameters by `are.par.valid`

to simulated L-moments. See the source code for more details.

W.H. Asquith

`lmoms`

, `pmoms`

, `par2qua`

, `genci.simple`

, `qua2ci.cov`

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 | ```
## Not run:
# It is well known that standard deviation (sigma) of the
# sample mean is equal to sigma/sample_size. Let is look at the
# quantile distribution of the median (f=0.5)
mean <- 0; sigma <- 100
parent <- vec2par(c(mean,sigma), type='nor')
CI <- qua2ci.simple(0.5, parent, n=10, nsim=20)
# Theoretrical sample mean sigma = 100/10 = 10
# L-moment theory: L-scale * sqrt(pi) = sigma
# Thus, it follows that the quantity
CI$elmoms$lambdas[2]/sqrt(pi)
# approaches 10 as nsim --> Inf.
## End(Not run)
# Another example.
D <- c(123, 34, 4, 654, 37, 78, 93, 95, 120) # fake sample
lmr <- lmoms(D) # compute the L-moments of the data
WEI <- parwei(lmr) # estimate Weibull distribution parameters
CI <- qua2ci.simple(0.75,WEI,20, nsim=20, level=0.95)
# CI contains the estimate 95-percent confidence interval for the
# 75th-percentile of the parent Weibull distribution for 20 sample size 20.
## Not run:
pdf("Substantial_qua2ci_example.pdf")
level <- 0.90; cilo <- (1-level)/2; cihi <- 1 - cilo
para <- lmom2par(vec2lmom(c(180,50,0.75)), type="gev")
A <- qua2ci.simple(0.98, para, 30, edist="gno", level=level, nsim=3000)
Apara <- A$epara; Aenv <- A$empdist
Bpara <- lmom2par(A$elmoms, type="aep4")
lo <- log10(A$lwr); hi <- log10(A$upr)
xs <- 10^(seq(lo-0.2, hi+0.2, by=0.005))
lo <- A$lwr; hi <- A$upr; xm <- A$true; sbar <- mean(Aenv$simquas)
dd <- density(Aenv$simquas, adjust=0.5)
pk <- max(dd$y, dlmomco(xs, Apara), dlmomco(xs, Bpara))
dx <- dd$x[dd$x >= Aenv$empir.dist.lower & dd$x <= Aenv$empir.dist.upper]
dy <- dd$y[dd$x >= Aenv$empir.dist.lower & dd$x <= Aenv$empir.dist.upper]
dx <- c(dx[1], dx, dx[length(dx)]); dy <- c(0, dy, 0)
plot(c(0), c(0), type="n", xlim=range(xs), ylim=c(0,pk),
xlab="X VALUE", ylab="PROBABILITY DENSITY")
polygon(dx, dy, col=8)
lines(xs, dlmomco(xs, Apara)); lines(xs, dlmomco(xs, Bpara), col=2, lwd=2)
lines(dd, lty=2, lwd=2, col=8)
lines(xs, dlmomco(xs, para), col=6); lines(c(xm,xm), c(0,pk), lty=4, lwd=3)
lines(c(lo,lo,NA,hi,hi), c(0,pk,NA,0,pk), lty=2)
xlo <- qlmomco(cilo, Apara); xhi <- qlmomco(cihi, Apara)
points(c(xlo, xhi), c(dlmomco(xlo, Apara), dlmomco(xhi, Apara)), pch=16)
xlo <- qlmomco(cilo, Bpara); xhi <- qlmomco(cihi, Bpara)
points(c(xlo, xhi), c(dlmomco(xlo, Bpara), dlmomco(xhi, Bpara)), pch=16, col=2)
lines(rep(Aenv$empir.dist.lwr, 2), c(0,pk), lty=3, lwd=2, col=3)
lines(rep(Aenv$empir.dist.upr, 2), c(0,pk), lty=3, lwd=2, col=3)
lines(rep(Aenv$bern.smooth.lwr,2), c(0,pk), lty=3, lwd=2, col=4)
lines(rep(Aenv$bern.smooth.upr,2), c(0,pk), lty=3, lwd=2, col=4)
cat(c( "F(true) = ", round(plmomco(xm, Apara), digits=2),
"; F(mean(sim), edist) = ", round(plmomco(sbar, Apara), digits=2), "\n"), sep="")
dev.off()
## End(Not run)
## Not run:
ty <- "nor" # try running with "glo" (to get the L-skew "fit", see below)
para <- lmom2par(vec2lmom(c(-180,70,-.5)), type=ty)
f <- 0.99; n <- 41; ns <- 1000; Qtrue <- qlmomco(f, para)
Qsim1 <- replicate(ns, qlmomco(f, lmom2par(lmoms(rlmomco(n, para)), type=ty)))
Qsim2 <- qua2ci.simple(f, para, n, nsim=ns, edist="gno")
Qbar1 <- mean(Qsim1); Qbar2 <- mean(Qsim2$empdist$simquas)
epara <- Qsim2$epara; FT <- plmomco(Qtrue, epara)
F1 <- plmomco(Qbar1, epara); F2 <- plmomco(Qbar2, epara)
cat(c( "F(true) = ", round(FT, digits=2),
"; F(via sim.) = ", round(F1, digits=2),
"; F(via edist) = ", round(F2, digits=2), "\n"), sep="")
# The given L-moments are highly skewed, but a Normal distribution is fit so
# L-skew is ignored. The game is deep tail (f=0.99) estimation. The true value of the
# quantile has a percentile on the error distribution 0.48 that is almost exactly 0.5
# (median = mean = symmetrical error distribution). A test run shows nice behavior:
# F(true) = 0.48; F(via sim.) = 0.49; F(via edist) = 0.5
# But another run with ty <- "glo" (see how 0.36 << [0.52, 0.54]) has
# F(true) = 0.36; F(via sim.) = 0.54; F(via edist) = 0.52
# So as the asymmetry becomes extreme, the error distribution becomes asymmetrical too.
## End(Not run)
``` |

lmomco documentation built on March 18, 2018, 1:45 p.m.

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.