lmom.ub, lmom2pwm),
Probability-Weighted Moment estimation (pwm.ub, pwm.pp, pwm.gev), parameter estimation for numerous familiar and not-so-familiar distributions (see section Distributions), and L-moment estimation for the same distributions from the parameters (lmom2par). In words, L-moments are derived from the expectations of order statistics and are linear with respect to the probability-weighted moments. L-moments are directly analogous to the well-known product moments; however, L-moments have many advantages including unbiasedness, robustness, and consistency with respect to the conventional product (central) moments (mean, standard deviation, skew, kurtosis, ...). L-moments have particular internal relations to themselves and boundness (see
are.lmom.valid). In general, this package is oriented around the FORTRAN algorithms of J.R.M. Hosking, and the nomenclature for many of the functions parallels that of the Hosking library. Extensions are made. Additionally, recent developments by Robert Serfling and Peng Xiao have extended L-moments into multivariate space---the L-comoments. The sample L-comoments (Lcomoment.Lk12) are implemented here for an arbitrary number of random variables and moment order value. The L-comoments are considered experimental in this package; however, for algorithm evaluation the diagonal of the L-comoment matrix (Lcomoment.matrix) produces conventional L-moments (lmom.ub) of the corresponding order. A complimentary addition to L-moment theory are the trimmed L-moments by Elsayed Elamir and Allan Seheult. The sample trimmed L-moments are supported by TLmoms and are useful tools for parameter estimation for some distributions including the Cauchy; the trimmed L-moments have further robust properties compared to the usual or ordinary L-moments.$$\lambda_r \equiv \frac{1}{r}\sum^{r-1}_{k=0} (-1)^k { r-1 \choose k } E[X_{r-k:r}]$$
where $r \ge 1$ is the order of the L-moment, and $E[X_{r-k:r}]$ is the expectation of the $r-k$ order statisic of a sample of size $r$. The expectation of an order statistic is
$$E[X_{j:r}] = \frac{r!}{(j-1)!(r-j)!}\int^1_0 X(F)\times F^{j-1}(1-F)^{r-j} \,\mathrm{d}F \mbox{.}$$
The theoretical L-moments for the supported distributions are computing from the above equations using numerical integration and the theoLmom function. However, each distribution has is own system of L-moment equations; these are used for parameter estimation by this package. The first four L-moments in terms of integrals are
$$\lambda_1 = \int^1_0 X(F) \,\mathrm{d}F \mbox{,}$$ $$\lambda_2 = \int^1_0 X(F) \times (2F - 1) \,\mathrm{d}F \mbox{,}$$ $$\lambda_3 = \int^1_0 X(F) \times (6F^2 - 6F + 1) \,\mathrm{d}F \mbox{, and}$$ $$\lambda_4 = \int^1_0 X(F) \times (20F^3 - 30F^2 + 12F - 1) \,\mathrm{d}F \mbox{.}$$
The L-moment ratios are the dimensionless quantities
$$\tau = \mbox{LCV} = \lambda_2/\lambda_1 = \mbox{coefficient of L-variation, }$$ $$\tau_3 = \mbox{TAU3} = \lambda_3/\lambda_2 = \mbox{L-skew, }$$ $$\tau_4 = \mbox{TAU4} = \lambda_4/\lambda_2 = \mbox{L-kurtosis, and}$$ $$\tau_5 = \mbox{TAU5} = \lambda_5/\lambda_2 = \mbox{not named.}$$
pwm2lmom or directly. In either case, the L-moments are computed from the sample order statistics $x_{1:n} \le x_{2:n} \le \dots x_{n:n}$. In the most compact notation, the sample L-moments are$$\hat{\lambda}_r = \frac{1}{n}\sum^n_{i=1}w^{(r)}_{i:n} x_{i:n} \mbox{,}$$
where $w^{(r)}_{i:n}$ are weight factors. The weight factors are well illustrated in figure 2.6 of Hosking and Wallis (1997).
$$w^{(r)}_{i:n} = \sum_{j=0}^{min{i-1,r-1}} (-1)^{r-1-j} \frac{{r-1 \choose j}{r-1+j \choose j}{i-1 \choose j}} {{n-1 \choose j}} \mbox{.}$$
theoTLmoms). The $r \ge 1$ TL-moment is defined for trim level $t \ge 0$ as follows$$\lambda^{(t)}_r = \frac{1}{r} \sum^{r-1}_{k=0} \frac{(-1)^k {r-1 \choose k}(r+2t)!}{(r+t-k-1)!(t+k)!} \int^1_0 X(F) \times F^{r+t-k-1}(1-F)^{t+k} \,\mathrm{d}F \mbox{.}$$
It is important to note that when $t=0$ the ordinary $r$th L-moment is formed. The first four theoretical TL-moments for $t=1$ are
$$\lambda^{(1)}_1 = 6 \int^1_0 X(F) \times F \times (1-F) \,\mathrm{d}F \mbox{,}$$ $$\lambda^{(1)}_2 = 6 \int^1_0 X(F) \times F \times (1-F)(2F - 1) \,\mathrm{d}F \mbox{,}$$ $$\lambda^{(1)}_3 = \frac{20}{3} \int^1_0 X(F) \times F \times (1-F)(5F^2 - 5F + 1) \,\mathrm{d}F \mbox{,}$$ $$\lambda^{(1)}_4 = \frac{15}{2} \int^1_0 X(F) \times F \times (1-F)(14F^3 - 21F^2 + 9F - 1) \,\mathrm{d}F \mbox{,}$$
$$\hat{\lambda}^{(t)}_r = \frac{1}{r}\sum^{n-t}_{i=t+1} \left[ \frac{\sum\limits^{r-1}_{k=0}{ (-1)^k {r-1 \choose k} {i-1 \choose r+t-1-k} {n-i \choose t+k} }}{{n \choose r+2t}} \right] x_{i:n}$$
where $t$ represents the trimming level of the $t$-largest of $t$-smallest values, $r$ represents the order of the L-moments, $n$ represents the sample size, and $x_{i:n}$ is the $i$th sample order statistics ($x_{1:n} \le x_{2:n} \le \dots x_{n:n}$). Sample asymmetrical TL-moments are supported by the TLmoms function.
$$\lambda_{[12]} = \frac{1}{k} \sum^{k-1}_{j=0} (-1)^k {k-1 \choose j} E[X^{(12)}_{[k-j:k]}] \mbox{,}$$ where $E[X^{(12)}_{[k-j:k]}]$ is the expected value of the $X_{[k-j:n]}$ concomitant of $X2$ with respect to $X1$.
$$E[X^{(12)}_{[k-j:k]}] = r {n \choose r} \int^1_0 \int^1_0 X(F_1) \times (F_2)^{r-1}(1-F_2)^{n-r} \,\mathrm{d}F_1 \,\mathrm{d}F_2 \mbox{,}$$
where $X(F_1)$ is quantile function of random variable $X1$ and $F_2$ is nonexceedance probability of random variable $X2$.
The sample L-comoment (Lcomoment.Lk12) of random variable $X1$ are computed from the concomitants of $X2$. The concomitants inturn are used in a weighted summation and expectation calculation to compute the L-comoment of $X1$ with respect to $X2$.
$$\hat{\lambda}_{k[12]} = n^{-1}\sum_{r=1}^{n} w^{(k)}_{r:n} x^{(12)}_{[r:n]} \mbox{,}$$
where $x^{(12)}_{[r:n]}$ are the sample concomitant statistics and $w^{(k)}_{r:n}$ is
$$w^{(k)}_{r:n} = \sum_{j=0}^{min{r-1,k-1}} (-1)^{k-1-j} \frac{{k-1 \choose j}{r-1+j \choose j}{r-1 \choose j}} {{n-1 \choose j}} \mbox{,}$$
The sample L-comoment matrix (Lcomoment.matrix) for $d$-random variables is
$$\mbox{\boldmath $\Lambda_k$} = (\hat{\lambda}_{k[ij]})$$
computed over the pairs ($X^{(i)},X^{(j)}$) where $1 \le i \le j \le d$.
CCC in the documentation). At present (2006), 13 distributions are supported by stable algorithms for parameter estimation using L-moments (parCCC, such as parexp), L-moment estimation using parameters (lmomCCC, such as lmomexp), cumulative distribution function (nonexceedance probability as a function of the variable), and quantile distribution function (variable as a function of nonexceedance probability). A dispatcher for parameter estimation from the L-moments is lmom2par. A dispatcher for L-moment estimation from the parameters is par2lmom. The cumulative distribution functions are cdfCCC, such as cdfexp; a dispatcher to the cumulative distribution functions is par2cdf. The quantile functions are quaCCC, such as quaexp; a dispatcher to the quantile functions is par2qua. For TL-moments, a TL is prepended to CCC.The distributions supported are the Cauchy, Exponential, Gamma, Generalized Extreme Value, Generalized Lambda, Generalized Logistic, Generalized Normal, Generalized Pareto, Gumbel, Kappa, Normal, Pearson Type III, and Wakeby. Some of these distributions (Exponential, Gamma, and Normal) have functional implementation within the standard R-package distribution. As this package in part mirrors the Hosking FORTRAN libraries in wide spread use by the L-moment user-community (in particular environmental and hydrologic sciences---at least those familiar to the author), these three distributions are implemented here in a parallel function context. However, it is important to note that R functions are used for these three
cau = Cauchy distribution (two parameters)---TL-moment implementation
exp = Exponential distribution (two parameters)---ordinary L-moment implementation
gam = Gamma distribution (two parameters)---ordinary L-moment implementation
gev = Generalized Extreme Value distribution (three parameters)---ordinary L-moment implementation
gld = Generalized Lambda distribution (four parameters)---ordinary L-moment implementation.
TLgld = Generalized Lambda distribution (four parameters)---TL-moment implementation.
glo = Generalized Logistic distribution (three parameters)---ordinary L-moment implementation
gno = Generalized Normal (log-Normal) distribution (three parameters)---ordinary L-moment implementation
gpa = Generalized Pareto distribution (three parameters)---ordinary L-moment implementation
TLgpa = Generalized Pareto distribution (three parameters)---TL-moment implementation
gum = Gumbel distribution (two parameters)---ordinary L-moment implementation
kap = Kappa distribution (four parameters)---ordinary L-moment implementation
nor = Normal distribution (two parameters)---ordinary L-moment implementation
pe3 = Pearson Type III distribution (three parameters)---ordinary L-moment implementation
wak = Wakeby distribution (five parameters)---ordinary L-moment implementation
wei = Weibull distribution (three parameters)---ordinary L-moment implementation
Parameters for the distribution are placed into a particular object format (see vec2par, lmom2par, or parexp). The parameter object (simply an R list) in turn can be passed as an argument to the distribution functions. A broader intent of this package is to support modular code design when users are heavily involved in distributional analysis. Therefore, this package contains a number of ancillary functions such as are.par.valid or is.CCC, where CCC is the three character syntax as previously shown to assist users in building sophisticated tools in R.
lmrdia and
plotlmrdia)---namely, the construction of the L-skew and L-kurtosis diagram. On the diagram the theoretical trajectories of most of the aforementioned distributions are shown. These diagrams are difficult to explain here, but are well documented in the literature (lmom.references).lmom.diff function computes the difference between the L-moments derived from a parameterized distribution and the L-moments as computed from the data. From the author's experience, construction of freq.curve.CCC are available for ease of use. There also is a freq.curve.all function that computes the frequency curve for all the distributions given an L-moment object. The curves require vectors of nonexceedance probabilities (see nonexceeds). Related to nonexceedance probabilities, note that for this documentation nonexceedance probability is shown as $0 \le F \le 1$; however, some distributions might not be valid at $F = 0$ or $F = 1$. Finally, the functions lmom.test.all, which dispatches to distribution-specific functions following the lmom.test.CCC naming convention and provides user-level output to help evaluate the algorithms of this package.# One has the following peak streamflow values in cubic meters per second
data <- c(123,2250,543,178,67,5100,248,1500,342,329,
543,980,1020,4502,3406,856,297)
# Compute the unbiased L-moments of the data--high L-skew.
# This data clearly is not Normally distributed.
lmr <- lmom.ub(data)
# One style of parameter estimation for a Kappa distribution
Kappapar <- lmom2par(lmr,type='kap')
# Another style of parameter estimation for Normal distribution
Normalpar <- parnor(lmr)
# Vector of useful nonexceedance probabilities
F <- nonexceeds()
# Generate Kappa frequency curve
Qk <- freq.curve.kap(nonexceeds(),Kappapar)
# Generate Normal frequency curve
Qn <- freq.curve.nor(nonexceeds(),Normalpar)
# Plot them up
plot(F,Qk,type="n",ylim=c(-6000,24000))
lines(F,Qk)
lines(F,Qn,col=2)
X1 <- data
# Generate some related data
X2 <- abs(rnorm(length(data)))*data
L2 <- Lcomoment.matrix(data.frame(RanVar=X1,AnotherRanVar=X2),k=2)
L2
# Compute the convential L-moments of variable 1 and 2
X1lmr <- lmom.ub(X1)
X2lmr <- lmom.ub(X2)
# Show that the diagonal of the Lcomoment matrix equals the
# conventional moments of same order (2nd order in this case).
print(c(X1lmr$L2,L2$matrix[1,1]))
print(c(X2lmr$L2,L2$matrix[2,2]))
# Compute the L-correlation values
Lrho <- Lcomoment.correlation(L2)
Lrho
# Compare the off-diagonal terms to the conventional
# correlation coefficient. The off-diagonal terms will
# not be equal or equal in value to the conventional
# correlation coefficient.
cor(X1,X2)Run the code above in your browser using DataLab