lmomco (version 2.2.5)

genci.simple: Generate (Estimate) Confidence Intervals for Quantiles of a Parent Distribution

Description

This function estimates the lower and upper limits of a specified confidence interval for a vector of nonexceedance probabilities $F$ of a specified parent distribution [quantile function $Q(F,\theta)$ with parameters $\theta$] using Monte Carlo simulation. The $F$ are specified by the user. The user also provides $\Theta$ of the parent distribution (see lmom2par). This function is a wrapper on qua2ci.simple; please consult the documentation for that function for further details of the simulations.

Usage

genci.simple(para, n, f=NULL, level=0.90, edist="gno", nsim=1000, expand=FALSE, verbose=FALSE, showpar=FALSE, quiet=FALSE)

Arguments

para

The parameters from lmom2par or similar.

n

The sample size for each Monte Carlo simulation will use.

f

Vector of nonexceedance probabilities ($0 \le f \le 1$) of the quantiles for which the confidence interval are needed. If NULL, then the vector as returned by nonexceeds is used.

level

The confidence interval ($0 \le $ level $ < 1$). The interval is specified as the size of the interval. The default is 0.90 or the 90th percentile. The function will return the 5th ($(1-0.90)/2$) and 95th ($1-(1-0.90)/2$) percentile cumulative probability of the error distribution for the parent quantile as specified by the nonexceedance probability argument (f). This argument is passed unused to qua2ci.simple.

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 (lmomgno). For example, if the L-skew ($\tau_4$) or L-kurtosis ($\tau_4$) values depart considerably from those of the Normal ($\tau_3 = 0$ and $\tau_4 = 0.122602$), then the Generalized Normal or some alternative distribution would likely provide more reliable confidence interval estimation. This argument is passed unused to qua2ci.simple.

nsim

The number of simulations (replications) for the sample size n to perform. Much larger simulation numbers are recommended---see discussion about qua2ci.simple. This argument is passed unused to qua2ci.simple. Users are encouraged to experiment with qua2ci.simple to get a feel for the value of edist and nsim.

expand

Should the returned values be expanded to include information relating to the distribution type and L-moments of the distribution at the corresponding nonexceedance probabilities---in other words the information necessary to reconstruct the reported confidence interval. The default is FALSE. If expand=FALSE then a single data.frame of the lower and upper limits along with the true quantile value of the parent is returned. If expand=TRUE, then a more complicated list containing multiple data.frames is returned.

verbose

The verbosity of the operation of the function. This argument is passed unused to qua2ci.simple.

showpar

The parameters of the edist for each simulation for each $F$ value passed to qua2ci.simple are printed. This argument is passed unused to qua2ci.simple.

quiet

Suppress incremental counter for a count down of the $F$ values.

Value

An R data.frame or list is returned (see discussion of argument expand). The following elements could be available.

See Also

genci, gen.freq.curves

Examples

## Not run: 
# # For all these examples, nsim is way too small.
# mean   <- 0; sigma <- 100
# parent <- vec2par(c(mean,sigma), type='nor') # make parameter object
# f      <- c(0.5, 0.8, 0.9, 0.96, 0.98, 0.99) # nonexceed probabilities
# # nsim is small for speed of example not accuracy.
# CI     <- genci.simple(parent, n=10, f=f, nsim=20); FF <- CI$nonexceed
# plot( FF, CI$true, type='l', lwd=2)
# lines(FF, CI$lwr, col=2); lines(FF, CI$upr, col=3)
# 
# pdf("twoCIplots.pdf")
# # The qnorm() call has been added to produce "normal probability"
# # paper on the horizonal axis. The parent is heavy-tailed.
# GEV  <- vec2par(c(10000,1500,-0.3), type='gev') # a GEV distribution
# CI   <- genci.simple(GEV, n=20, nsim=200, edist='gno')
# ymin <- log10(min(CI$lwr[! is.na(CI$lwr)]))
# ymax <- log10(max(CI$upr[! is.na(CI$upr)]))
# qFF  <- qnorm(CI$nonexceed) 
# plot( qFF, log10(CI$true), type='l', ylim=c(ymin,ymax),lwd=2)
# lines(qFF, log10(CI$lwr), col=2); lines(qFF, log10(CI$upr), col=3)
# # another error distribution model
# CI   <- genci.simple(GEV, n=20, nsim=200, edist='aep4')
# lines(qFF,log10(CI$lwr),col=2,lty=2); lines(qFF,log10(CI$upr),col=3,lty=2)
# dev.off() # ## End(Not run)