qua2ci.cov: Estimate a Confidence Interval for Quantiles of a Parent Distribution using Sample Variance-Covariances of L-moments

Description

This function estimates the lower and upper limits of a specified confidence interval for aribitrary quantile values for a sample \(x\) and a specified distribution form. The estimation is based on the sample variance-covariance structure of the L-moments (lmoms.cov) through a Monte Carlo approach. The quantile values, actually the nonexceedance probabilities (\(F\) for \(0 \le F \le 1\)), are specified by the user. The user provides type of parent distribution distribution and this form which will be fitted internal to the function.

Usage

qua2ci.cov(x,f, type=NULL, nsim=1000,
                interval=c("confidence", "none"), level=0.90,
                asnorm=FALSE, altlmoms=NULL, flip=NULL, dimless=TRUE,
                usefastlcov=TRUE, nmom=5, getsimlmom=FALSE, verbose=FALSE, ...)

Arguments

A real value vector.

Nonexceedance probabilities (\(0 \le F \le 1\)) of the quantiles for which the confidence interval is needed.

type

Three character distribution type (for example, type='gev').

nsim

The number of simulations to perform. Large numbers produce more refined confidence limit estimates at the cost of CPU time. The default is anticipated to be large enough to semi-quantitatively interpret results without too much computational delay. Larger simulation numbers are recommended.

interval

The type of interval to compute. If "none", then the simulated quantiles are returned at which point only the first value in \(f\) or f[1] will be considered but a warning will be issued to remind the user. This option is nice for making boxplots of the quantile distribution.

level

The confidence interval (\(0 \le \) level \( < 1\)). The interval is specified as the size of the interval for which 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 simulated quantile distribution as specified by the nonexceedance probability argument.

asnorm

Use the mean and standard deviation of the simulated quantiles as parameters of the Normal distribution to estimate the confidence interval. Otherwise, a Bernstein polynomial approximation (dat2bernqua) to the empirical distribution of the simulated quantile distribution is used.

altlmoms

Alternative L-moments to rescale the simulated L-moments from the variance-covariance structure of the sample L-moments in x. These L-moments need to be an lmomco package L-moment object (e.g. lmoms). The presence of alternative L-moments will result in dimless=TRUE.

flip

A flipping or reflection value denoted as \(\eta\). The values in x are flipped by this value (\(y = \eta - x\)) and analysis proceeds with flipped information, and then results are flipped back just prior to returning values with the exception that if getsimlmom=TRUE then the simultated L-moments are in “flipped space.”

dimless

Perform the simulations in dimensionless space meaning that values in x are converted by \(y = (x-\lambda_1)/\lambda_2\) and simulation based on \(y\) and scale is returned on output according to the L-moments of x or the alternative L-moments in altlmoms. Scale is returned to the simulated L-moments, if returned by getsimlmom=TRUE, which is not fully parallel with the returned behavior when flipping is involved.

usefastlcov

A logical to use the function Lmomcov() from the Lmoments package to compute the sample variance-covariance matrices and not the much slower function lmoms.cov in the lmomco package.

nmom

The number of L-moments involved. This argument needs to be high enough to permit parameterization of the distribution in type but computational effort increases as nmom gets large. This option is provided in conjunction with getsimlmom=TRUE to be able to get a “wider set” of simulated L-moments returned than precisely required by the distribution. Also, some distributions might as part of their specific fitting algorithms, require inspection of higher L-moments than seemingly required than their numer of parameters suggests.

getsimlmom

A logical controlling whether the simulated L-moment matrix having nsim rows and nmom columns is returned instead of confidence limits.

verbose

The verbosity of the operation of the function.

...

Additional arguments to pass such as to lmom2par.

Value

An R data.frame is returned.

lwr

The lower value of the confidence interval having nonexceedance probability equal to \((1-\)level\()/2\).

fit

The fit of the quantile based on the L-moments of x and possibly by reflection controlled by flip or based on the alternative L-moments in altlmoms and again by the reflection controlled by flip.

upr

The upper value of the confidence interval having nonexceedance probability equal to \(1-(1-\)level\()/2\).

qua_med

The median of the simulated quantiles.

qua_mean

The mean of the simulated quantiles for which the median and mean should be very close if the simulation size is large enough and the quantile distribution is symmetrical.

qua_var

The variance (\(\sigma^2(F)\)) of the simulated quantiles.

qua_lam2

The L-scale (\(\lambda_2(F)\)) of the simulated quantiles for which \(\sigma^2(F) \approx \pi\times\lambda^2_2(F)\).

Examples

Run this code

# NOT RUN {
samsize <- 128; nsim <- 2000; f <- 0.999
wei <- parwei(vec2lmom(c(100,75,-.3)))
set.seed(1734); X <- rlmomco(samsize, wei); set.seed(1734)
tmp <- qua2ci.cov(X, f, type="wei", nsim=nsim)
print(tmp) # show results of one 2000 replicated Monte Carlo
# nonexceed     lwr    fit    upr  qua_med  qua_mean  qua_var  qua_lam2
#     0.999   310.4  333.2  360.2    333.6     334.3    227.3    8.4988
set.seed(1734)
qf <- qua2ci.cov(X, f, type="wei", nsim=nsim, interval="none") # another
boxplot(qf)
message(" quantile variance: ", round(tmp$qua_var,  digits=2),
        " compared to ", round(var(qf, na.rm=TRUE), digits=2))
set.seed(1734)
genci.simple(wei, n=samsize, f=f)
# nonexceed     lwr    fit    upr  qua_med  qua_mean  qua_var  qua_lam2
#     0.999   289.7  312.0  337.7    313.5     313.6    213.5    8.2330

#----------------------------------------
# Using X from above example, demonstrate that using dimensionless
# simulation that the results are the same.
set.seed(145); qua2ci.cov(X, 0.1, type="wei") # both outputs same
set.seed(145); qua2ci.cov(X, 0.1, type="wei", dimless=TRUE)
# nonexceed     lwr    fit    upr  qua_med  qua_mean  qua_var  qua_lam2
#       0.1  -78.62 -46.01 -11.39   -43.58    -44.38   416.04     11.54

#----------------------------------------
# Using X again, demonstration application of the flip and notice that just
# simple reversal is occurring and that the Weibull is a reversed GEV.
eta <- 0
set.seed(145); qua2ci.cov(X, 0.9, type="wei", nsim=nsim)
# nonexceed     lwr    fit    upr  qua_med  qua_mean  qua_var  qua_lam2
#       0.9   232.2  244.2  255.9    244.3     244.1    51.91    4.0635
set.seed(145); qua2ci.cov(X, 0.9, type="gev", nsim=nsim, flip=eta)
# nonexceed     lwr    fit    upr  qua_med  qua_mean  qua_var  qua_lam2
#       0.9   232.2  244.2  256.2    244.2     244.3    53.02    4.1088
# The values are slightly different, which likely represents a combination
# of numerics of the variance-covariance matrix because the Monte Carlo
# is seeded the same.

#----------------------------------------
# Using X again, removed dimension and have the function add it back.
lmr <- lmoms(X); Y <- (X - lmr$lambdas[1])/lmr$lambdas[2]
set.seed(145); qua2ci.cov(Y, 0.9, type="wei", altlmoms=lmr, nsim=nsim)
# nonexceed     lwr    fit    upr  qua_med  qua_mean  qua_var  qua_lam2
#       0.9   232.2  244.2  255.9    244.3     244.1    51.91   4.0635
# }

Run the code above in your browser using DataLab

Description

Usage

Arguments

Value

See Also

Examples