lmrp.RdComputes the \(L\)-moments or trimmed \(L\)-moments
of a probability distribution
given its cumulative distribution function (for function lmrp)
or quantile function (for function lmrq).
lmrp(pfunc, ..., bounds=c(-Inf,Inf), symm=FALSE, order=1:4,
ratios=TRUE, trim=0, acc=1e-6, subdiv=100, verbose=FALSE)
lmrq(qfunc, ..., symm=FALSE, order=1:4, ratios=TRUE, trim=0,
acc=1e-6, subdiv=100, verbose=FALSE)Cumulative distribution function.
Quantile function.
Arguments to pfunc or qfunc.
Either a vector of length 2, containing the lower and upper bounds of the distribution, or a function that calculates these bounds given the distribution parameters as inputs.
For lmrq, a logical value indicating whether
the distribution is symmetric about its median.
For lmrp, either the logical value FALSE or NA
to indicate that the distribution is not symmetric,
or a numeric value to indicate that the distribution is symmetric
and that the specified value is the center of symmetry.
If the distribution is symmetric, odd-order \(L\)-moments are exactly zero and the symmetry is used to slightly speed up the computation of even-order \(L\)-moments.
Orders of the \(L\)-moments and \(L\)-moment ratios to be computed.
Logical. If FALSE, \(L\)-moments are computed;
if TRUE (the default), \(L\)-moment ratios are computed.
Degree of trimming. If a single value, symmetric trimming of the specified degree will be used. If a vector of length 2, the two values indicate the degrees of trimming at the lower and upper ends of the “conceptual sample” (Elamir and Seheult, 2003) of order statistics that is used to define the trimmed \(L\)-moments.
Requested accuracy. The function will try to achieve this level of accuracy, as relative error for \(L\)-moments and absolute error for \(L\)-moment ratios.
Maximum number of subintervals used in numerical integration.
Logical. If FALSE, only the values of the \(L\)-moments
and \(L\)-moment ratios are returned. If TRUE, more details of the
numerical integration are returned: see “Value” section below.
Computations use expressions in Hosking (2007):
eq. (7) for lmrp, eq. (5) for lmrq.
Integrals in those expressions are computed by numerical integration.
If verbose is FALSE and ratios is FALSE,
a numeric vector containing the \(L\)-moments.
If verbose is FALSE and ratios is TRUE,
a numeric vector containing the \(L\)-moments (of orders 1 and 2)
and \(L\)-moment ratios (of orders 3 and higher).
If verbose is TRUE, a data frame with columns as follows:
\(L\)-moments (if ratios is FALSE),
or \(L\)-moments and \(L\)-moment ratios (if ratios is TRUE).
Estimate of the absolute error in the computed value.
"OK" or a character string giving the error message
resulting from the numerical integration.
pfunc and qfunc can be either the standard R form 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 R form, 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.
In package lmom versions 1.6 and earlier, the “Details” section stated that
“Integrals in those expressions are computed by numerical integration,
using the R function integrate”.
As of version 2.0, numerical integration uses an internal function that directly calls
(slightly modified versions of) Fortran routines in QUADPACK (Piessens et al. 1983).
R's own integrate function uses C code “based on” the QUADPACK routines,
but in R versions 2.12.0 through 3.0.1 did not in every case reproduce the results
that would have been obtained with the Fortran code (this is R bug PR#15219).
Elamir, E. A. H., and Seheult, A. H. (2003). Trimmed L-moments. Computational Statistics and Data Analysis, 43, 299-314.
Hosking, J. R. M. (2007). Some theory and practical uses of trimmed L-moments. Journal of Statistical Planning and Inference, 137, 3024-3039.
Piessens, R., deDoncker-Kapenga, E., Uberhuber, C., and Kahaner, D. (1983). Quadpack: a Subroutine Package for Automatic Integration. Springer Verlag.
Arguments bounds, symm, order,
ratios, trim, 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).
## Generalized extreme-value (GEV) distribution
## - three ways to get its L-moments
lmrp(cdfgev, c(2,3,-0.2))
#> lambda_1 lambda_2 tau_3 tau_4
#> 4.4634457 2.5967856 0.3050929 0.2180272
lmrq(quagev, c(2,3,-0.2))
#> lambda_1 lambda_2 tau_3 tau_4
#> 4.4634457 2.5967856 0.3050929 0.2180272
lmrgev(c(2,3,-0.2), nmom=4)
#> lambda_1 lambda_2 tau_3 tau_4
#> 4.4634457 2.5967856 0.3050929 0.2180272
## GEV bounds specified as a vector
lmrp(cdfgev, c(2,3,-0.2), bounds=c(-13,Inf))
#> lambda_1 lambda_2 tau_3 tau_4
#> 4.4634457 2.5967856 0.3050929 0.2180272
## 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)
#> lambda_1 lambda_2 tau_3 tau_4
#> 4.4634457 2.5967856 0.3050929 0.2180272
## 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)
#> lambda_1 lambda_2 tau_3 tau_4
#> 4.4634457 2.5967856 0.3050929 0.2180272
## Normal distribution
lmrp(pnorm)
#> lambda_1 lambda_2 tau_3 tau_4
#> 3.905026e-17 5.641896e-01 5.869438e-17 1.226017e-01
lmrp(pnorm, symm=0)
#> lambda_1 lambda_2 tau_3 tau_4
#> 0.0000000 0.5641896 0.0000000 0.1226017
lmrp(pnorm, mean=2, sd=3, symm=2)
#> lambda_1 lambda_2 tau_3 tau_4
#> 2.0000000 1.6925688 0.0000000 0.1226017
# For comparison, the exact values
lmrnor(c(2,3), nmom=4)
#> lambda_1 lambda_2 tau_3 tau_4
#> 2.0000000 1.6925688 0.0000000 0.1226017
# Many L-moment ratios of the exponential distribution
# This may warn that "the integral is probably divergent"
lmrq(qexp, order=3:20)
#> Warning: in computation of L-moment of order 14: the integral is probably divergent
#> Warning: in computation of L-moment of order 15: the integral is probably divergent
#> Warning: in computation of L-moment of order 16: the integral is probably divergent
#> Warning: in computation of L-moment of order 20: the integral is probably divergent
#> tau_3 tau_4 tau_5 tau_6 tau_7 tau_8
#> 0.333333333 0.166666667 0.100000000 0.066666667 0.047619048 0.035714286
#> tau_9 tau_10 tau_11 tau_12 tau_13 tau_14
#> 0.027777778 0.022222222 0.018181818 0.015151515 0.012820513 0.010989011
#> tau_15 tau_16 tau_17 tau_18 tau_19 tau_20
#> 0.009523810 0.008333333 0.007352941 0.006535948 0.005847953 0.005263158
# ... 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))
#> exact value abs.error message
#> tau_3 0.333333333 0.333333333 4.220333e-07 OK
#> tau_4 0.166666667 0.166666667 7.542545e-08 OK
#> tau_5 0.100000000 0.100000000 4.527430e-08 OK
#> tau_6 0.066666667 0.066666667 3.037147e-08 OK
#> tau_7 0.047619048 0.047619048 2.235705e-08 OK
#> tau_8 0.035714286 0.035714286 1.850671e-08 OK
#> tau_9 0.027777778 0.027777778 1.829525e-08 OK
#> tau_10 0.022222222 0.022222222 2.254327e-08 OK
#> tau_11 0.018181818 0.018181818 3.338318e-08 OK
#> tau_12 0.015151515 0.015151515 5.467703e-08 OK
#> tau_13 0.012820513 0.012820513 9.293015e-08 OK
#> tau_14 0.010989011 0.010989011 1.590551e-07 the integral is probably divergent
#> tau_15 0.009523810 0.009523810 2.717775e-07 the integral is probably divergent
#> tau_16 0.008333333 0.008333333 4.649678e-07 the integral is probably divergent
#> tau_17 0.007352941 0.007352941 3.013771e-08 OK
#> tau_18 0.006535948 0.006535948 4.385433e-08 OK
#> tau_19 0.005847953 0.005847953 6.413495e-08 OK
#> tau_20 0.005263158 0.005263158 9.704056e-08 the integral is probably divergent
# Of course, sometimes the integral really is divergent
if (FALSE) { # \dontrun{
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)
lmrp(pcauchy)
#> Warning: in computation of L-moment of order 2: maximum number of subdivisions reached
#> Warning: in computation of L-moment of order 4: maximum number of subdivisions reached
#> lambda_1 lambda_2 tau_3 tau_4
#> 4.812893e-14 3.573871e+01 1.334866e-15 9.674570e-01
# Trimmed L-moments for Cauchy distribution are finite
lmrp(pcauchy, symm=0, trim=1)
#> lambda(1,1)_1 lambda(1,1)_2 tau(1,1)_3 tau(1,1)_4
#> 0.0000000 0.6978272 0.0000000 0.3428084
# Works for discrete distributions too, but often requires
# a larger-than-default value of 'subdiv'
lmrp(ppois, lambda=5, subdiv=1000)
#> lambda_1 lambda_2 tau_3 tau_4
#> 5.00000007 1.24548015 0.07612672 0.11879552