disfitqua: Fit a Distribution using Minimization of Available Quantiles

Description

This function fits a distribution to available quantiles (or irregular quantiles) through $n$-dimensional minimization using the optim function. The objective function forms are either root mean-square error (RMSE) or mean absolute deviation (MAD), and the objective functions are expected to result in slightly different estimates of distribution parameters. The RMSE form ($\sigma_{\mathrm{RMSE}}$) is defined as $$\sigma_{\mathrm{RMSE}} = \biggl[ \frac{1}{m}\,\sum_{i=1}^m [x_o(f_i) - \hat{x}(f_i)]^2\biggr]^{1/2}\mbox{,}$$ where $m$ is the length of the vector of $o$bserved quantiles $x_o(f_i)$ for nonexceedance probability $f_i$ for $i \in 1, 2, \cdots, m$, and $\hat{x}(f_i)$ for $i \in 1, 2, \cdots, m$ are quantile estimates based on the “current” iteration of the parameters for the selected distribution having $n$ parameters for $n \le m$. Similarly, the MAD form ($\sigma_{\mathrm{MAD}}$) is defined as $$\sigma_{\mathrm{MAD}} = \frac{1}{m}\,\sum_{i=1}^m | x_o(f_i) - \hat{x}(f_i) | \mbox{.}$$

The disfitqua function is not intended to be an implementation of the method of percentiles but rather is intended for circumstances in which the available quantiles are restricted to either the left or right tails of the distribution. It is evident that a form of the method of percentiles however could be pursued by disfitqua when the length of $x(f)$ is equal to the number of distribution parameters ($n = m$). The situation of $n < m$ however is thought to be the most common application.

The right-tail restriction is the general case in flood-peak hydrology in which the median and select quantiles greater than the median can be available from empirical studies (e.g. Asquith and Roussel, 2009) or rainfall-runoff models. The available quantiles suit engineering needs and thus left-tail quantiles simply are not available. This circumstance might appear quite unusual to users from most statistical disciplines but quantile estimates can exist from regional study of observed data. The Examples section provides further motivation and discussion.

Usage

disfitqua(x, f, objfun=c("rmse", "mad"),
                init.lmr=NA, init.para=NA, type=NA, verbose=FALSE, ... )

Arguments

The quantiles $x_o(f)$ for the nonexceedance probabilities in f.

The nonexceedance probabilities $f$ of the quantiles $x_o(f)$ in x.

objfun

The form of the objective function as previously described.

init.lmr

Optional initial values for the L-moments from which the initial starting parameters for the optimization will be determined. The optimizations by this function are not performed on the L-moments during the optimization. The form of init.lmr is that of an L-moment object from the lmomco package (e.g. lmoms).

init.para

Optional initial values for the parameters used for starting values for the optim function. If this argument is not set nor is init.lmr, then unrigorous estimates of the mean $\lambda_1$ and L-scale $\lambda_2$ are made from the available quantiles, higher L-moment ratios $\tau_r$ for $r \ge 3$ are set to zero, and the L-moments converted to the initial parameters.

type

The distribution type specified by the abbreviations listed under dist.list.

verbose

A logical switch on the verbosity of output.

...

Additional arguments to pass to the optim function.

Value

An R list is returned, and this list contains at least the following items:

type

The type of distribution in character format (see dist.list).

para

The parameters of the distribution.

source

Attribute specifying source of the parameters---“disfitqua”.

init.para

A vector of the initial parameters actually passed to the optim function to serve only as a reminder.

disfitqua

The returned list from the optim function. This list contains a repeat of the parameters, the value of the objective function ($\sigma_{\mathrm{RMSE}}$ or $\sigma_{\mathrm{MAD}}$), the interation count, and convergence status.

References

Asquith, W.H., and Roussel, M.C., 2009, Regression equations for estimation of annual peak-streamflow frequency for undeveloped watersheds in Texas using an L-moment-based, PRESS-minimized, residual-adjusted approach: U.S. Geological Survey Scientific Investigations Report 2009--5087, 48 p., http://pubs.usgs.gov/sir/2009/5087

Examples

Run this code

# NOT RUN {
# Suppose the following quantiles are estimated using eight equations provided by
# Asquith and Roussel (2009) for some watershed in Texas:
Q <- c(1480, 3230, 4670, 6750, 8700, 11000, 13600, 17500)
# These are real estimates from a suite of watershed properties but the watershed
# itself and location are not germane to demonstrate this function.
LQ <- log10(Q) # transform to logarithms of cubic feet per second
# Convert the averge annual return periods for the quantiles into probability
P <- T2prob(c(2, 5, 10, 25, 50, 100, 200, 500)); qP <- qnorm(P) # std norm variates
# The log-Pearson Type III (LPIII) is immensely popular for flood-risk computations.
# Let us compute LPIII parameters to the available quantiles and probabilities for
# the watershed. The log-Pearson Type III is "pe3" in the lmomco with logarithms.
par1 <- disfitqua(LQ, P, type="pe3", objfun="rmse") # root mean square error
par2 <- disfitqua(LQ, P, type="pe3", objfun="mad" ) # mean absolute deviation
# Now express the fitted distributions in forms of an LPIII.
LQfit1 <- qlmomco(P, par1); LQfit2 <- qlmomco(P, par2)

plot( qP, LQ, pch=5, xlab="STANDARD NORMAL VARIATES",
                     ylab="FLOOD QUANTILES, CUBIC FEET PER SECOND")
lines(qP, LQfit1, col=2); lines(qP, LQfit2, col=4) # red and blue lines
# }
# NOT RUN {
# Now demonstrate how a Wakeby distribution can be fit. This is an example of how a
# three parameter distribution might be fit and then the general L-moments secured for
# an alternative fit using a far more complicated distribution. The Wakeby for the
# above situation does not fit "out of the box." The types "gld", "aep4", and "kap"
# all with four parameters work with some serious CPUs burned for gld.
lmr1 <- theoLmoms(par1) # need five L-moments but lmompe3() only gives four,
# therefore must compute the L-moment by numerical integration provided by theoLmoms().
par3 <- disfitqua(LQ, P, type="wak", objfun="rmse", init.lmr=lmr1)
lines(qP, par2qua(P, par3), col=6, lty=2) # dashed line, par2qua alternative to qlmomco

# Finally, the initial L-moment equivalents and then the L-moments of the fitted
# distribution can be computed and compared.
par2lmom(vec2par(par3$init.para, type="wak"))$ratios # initial L-moments
par2lmom(vec2par(par3$para,      type="wak"))$ratios # final L-moments
# }

Run the code above in your browser using DataLab