lmrp: L-moments of a general probability distribution

Description

Computes the $L$-moments of a probability distribution given its cumulative distribution function (for function lmrp) or quantile function (for function lmrq).

Usage

lmrp(pfunc, ..., bounds = c(-Inf, Inf), symm = FALSE, order = 1:4,
     acc = 1e-6, subdiv = 100, verbose = FALSE)

lmrq(qfunc, ..., symm = FALSE, order = 1:4, acc = 1e-6,
     subdiv = 100, verbose = FALSE)

Value

If verbose=FALSE, a numeric vector containing the $L$-moments (of orders 1 and 2) and $L$-moment ratios (of orders 3 and higher). If verbose=TRUE, a data frame with columns as follows:
value$L$-moments (for orders 1 and 2) and $L$-moment ratios (for orders 3 and higher).
abs.errorEstimate of the absolute error in the computed value.
message"OK" or a character string giving the error message resulting from the numerical integration.

Arguments of cumulative distribution functions and quantile functions

pfunc and qfunc can be either the standard Rform of cumulative distribution function or quantile function (i.e. for a distribution with $r$ parameters, the first argument is the variate $x$ or the probability $p$ and the next $r$ arguments are the parameters of the distribution) or the cdf... or qua... forms used throughout the lmom package (i.e. the first argument is the variate $x$ or probability $p$ and the second argument is a vector containing the parameter values). Even for the Rform, however, starting values for the parameters are supplied as a vector start. If bounds is a function, its arguments must match the distribution parameter arguments of pfunc: either a single vector, or a separate argument for each parameter.

Warning

Arguments bounds, symm, order, acc, subdiv, and verbose cannot be abbreviated and must be specified by their full names (if abbreviated, the names would be matched to the arguments of pfunc or qfunc).

Details

$L$-moments are computed by numerical integration, using the Rfunction integrate. lmrq uses the expression $\lambda_r = \int_0^1 Q(u) P^{*}_{r-1}(u) du$ where $Q(u)$ is the quantile function and $P^{*}_{r-1}(u)$ is a shifted Legendre polynomial of order $r$ (Hosking, 1990, eq. (2.2)). lmrp uses the expression $\lambda_r = \int Z_r(F(x)) dx$ where $F(x)$ is the cumulative distribution function and $Z_r$ is a polynomial of degree $r$ defined in Hosking (2007, eq. (3.5)).

References

Hosking, J. R. M. (1990). $L$-moments: analysis and estimation of distributions using linear combinations of order statistics. Journal of the Royal Statistical Society, Series B, 52, 105-124. Hosking, J. R. M. (2007). Distributions with maximum entropy subject to constraints on their $L$-moments or expected order statistics. Journal of Statistical Planning and Inference, 137, 2870-2891.

Examples

Run this code

## Generalized extreme-value (GEV) distribution
## - three ways to get its L-moments
lmrp(cdfgev, c(2,3,-0.2))
lmrq(quagev, c(2,3,-0.2))
lmrgev(c(2,3,-0.2), nmom=4)

## GEV bounds specified as a vector
##
lmrp(cdfgev, c(2,3,-0.2), bounds=c(-13,Inf))
##
## GEV bounds specified as a function -- single vector of parameters
##
gevbounds <- function(para) {
  k <- para[3]
  b <- para[1]+para[2]/k
  c(ifelse(k<0, b,- Inf), ifelse(k>0, b, Inf))
}
lmrp(cdfgev, c(2,3,-0.2), bounds=gevbounds)
##
## GEV bounds specified as a function -- separate parameters
##
pgev <- function(x, xi, alpha ,k)
  pmin(1, pmax(0, exp(-((1-k*(x-xi)/alpha)^(1/k)))))
pgevbounds <- function(xi,alpha,k) {
  b <- xi+alpha/k
  c(ifelse(k<0, b, -Inf), ifelse(k>0, b, Inf))
}
lmrp(pgev, xi=2, alpha=3, k=-0.2, bounds=pgevbounds)


## Normal distribution
lmrp(pnorm)
lmrp(pnorm, symm=0)
lmrp(pnorm, mean=2, sd=3, symm=2)
# For comparison, the exact values
lmrnor(c(2,3), nmom=4)

# Many L-moment ratios of the exponential distribution
# This may warn that "the integral is probably divergent"
lmrq(qexp, order=3:20)

# ... nonetheless the computed values seem accurate:
# compare with the exact values, tau_r = 2/(r*(r-1)):
cbind(exact=2/(3:20)/(2:19), lmrq(qexp, order=3:20, verbose=TRUE))

# Of course, sometimes the integral really is divergent
lmrq(function(p) (1-p)^(-1.5))

# And sometimes the integral is divergent but that's not what
# the warning says (at least on the author's system)
lmrq(qcauchy)

# Works for discrete distributions too, but often requires
# a larger-than-default value of subdiv
lmrp(pbinom, size=20, prob=0.8, bounds=c(0,20), subdiv=1000)

Run the code above in your browser using DataLab