lmoms.cov | R Documentation |
Compute the distribution-free, variance-covariance matrix (\widehat{\mathrm{var}}(\lambda)
) of the sample L-moments (\hat\lambda_r
) or alternatively the sample probability-weighted moments (\hat\beta_k
, Elamir and Seheult, 2004, sec. 5). The \widehat{\mathrm{var}}(\lambda)
is defined by the matrix product
\widehat{\mathrm{var}}(\lambda) = \mathbf{C}\,\mathbf{\hat\Theta}\,\mathbf{C}^{\mathrm{T}}\mbox{,}
where the r \times r
matrix \mathbf{C}
for number of moments r
represents the coefficients of the linear combinations converting \beta_k
to \lambda_r
and the r
th row in the matrix is defined as
\mathbf{C}[r,]_{k{=}0:(r-1)} = (-1)^{(r-1-k)} {r-1 \choose k} {r-1+k \choose k}\mbox{,}
where the row is padded from the right with zeros for k < r
to form the required lower triangular structure. Elamir and Seheult (2004) list the \mathbf{C}
matrix for r = 4
.
Letting the falling factorial be defined (matching Elamir and Seheult's nomenclature) as
a^{(b)} = \Gamma(b+1) {a \choose b}\mbox{,}
and letting an entry in the \mathbf{\hat\Theta}
matrix denoted as \hat\theta_{kl}
be defined as
\hat\theta_{kl} = \hat\beta_k\hat\beta_l - \frac{A}{n^{(k+l+2)}}\mbox{,}
where \hat\beta_k
are again the sample probability-weighted moments and are computed by pwm
, and finally A
is defined as
A = \sum_{i=1}^{n-1}\sum_{j=i+1}^{n}
\bigl[ (i-1)^{(k)} (j-k-2)^{(l)} +
(i-1)^{(l)} (i-l-2)^{(k)}
\bigr] X_{i:n}X_{j:n}\mbox{,}
where X_{i:n}
are the sample order statistics for a sample of size n
.
Incidentally, the matrix \mathbf{\hat\Theta}
is the variance-covariance structure (\widehat{\mathrm{var}}
) of the \hat\beta
, thus \widehat{\mathrm{var}}(\beta) = \mathbf{\hat\Theta}
, which can be returned by a logical function argument (as.pwm=TRUE
) instead of \widehat{\mathrm{var}}(\lambda)
. The last example in Examples provides a demonstration.
lmoms.cov(x, nmom=5, as.pwm=FALSE, showC=FALSE,
se=c("NA", "lamse", "lmrse", "pwmse"), ...)
x |
A vector of data values. |
nmom |
The number of moments to compute. Default is 5. |
as.pwm |
A logical controlling whether the distribution-free, variance-covariance of sample probability-weighted moments ( |
showC |
A logical controlling whether the matrix |
se |
Compute standard errors ( |
... |
Other arguments to pass should they be needed (none were at first implementation). |
An R matrix
is returned. In small samples and substantially sized r
, one or more \hat\theta_{kl}
will be NaN
starting from the lower right corner of the matrix. The function does not test for this nor reduce the number of moments declared in nmom
itself. To reiterate, the square roots along the \widehat{\mathrm{var}}(\lambda)
diagonal are SE
for the respective L-moments.
Function lmoms.cov
was developed as a double check on the evidently separately developed r \le 4
(nmom
) implementations of \widehat{\mathrm{var}}(\lambda)
in packages Lmoments and nsRFA. Also the internal structure closely matches the symbolic mathematics by Elamir and Seheult (2004), but this practice comes at the expense of more than an order of magnitude slower execution times than say either of the functions Lmomcov()
(package Lmoments) or varLmoments()
(package nsRFA). For a high speed and recommended implementation, please use the Lmoments package by Karvanen (2016)—Karvanen extended this implementation to larger r
for the lmomco package.
For se="lmrse"
, the Taylor-series-based approximation is suggested by Elamir and Seheult (2004, p. 348) to estimate the variance of an L-moment ratio (\tau_r
for r \ge 3
) is based on structure of the variance of the ratio of two uniform variables in which the numerator is the r
th L-moment and the denominator is \lambda_2
:
\mathrm{var}(\tau_r) \cong
\biggl[ \frac{\mathrm{var}(\lambda_r)}{\mathrm{E}(\lambda_r)^2} + \frac{\mathrm{var}(\lambda_2)}{\mathrm{E}(\lambda_2)^2} - \frac{2\mathrm{cov}(\lambda_r,\lambda_2)}{\mathrm{E}(\lambda_r)\mathrm{E}(\lambda_2)} \biggr]
\biggl[\frac{\mathrm{E}(\lambda_r)}{\mathrm{E}(\lambda_2)} \biggr]^2\mbox{,}
where \mathrm{var}(\cdots)
are the along the diagonal of \widehat{\mathrm{var}}(\lambda)
and \mathrm{cov}(\cdots)
are the off-diagonal covariances. The expectations \mathrm{E}(\cdots)
are replaced with the sample estimates. Only for se="lmrse"
the SE
of the coefficient of L-variation (\tau_2^{SE}
) is computed but retained as an attribute (attr()
function) of the returned vector and not housed within the vector—the \lambda_2^{SE}
continues to be held in the 2nd position of the returned vector.
W.H. Asquith
Elamir, E.A.H., and Seheult, A.H., 2004, Exact variance structure of sample L-moments: Journal of Statistical Planning and Inference, v. 124, pp. 337–359.
Karvanen, Juha, 2016, Lmoments—L-moments and quantile mixtures: R package version 1.2-3, accessed February 22, 2016 at https://cran.r-project.org/web/packages/Lmoments/index.html
lmoms
, pwm
## Not run:
nsim <- 1000; n <- 10 # Let us compute variance of lambda_3
VL3sample <- mean(replicate(nsim, { zz <- lmoms.cov(rexp(n),nmom=3); zz[3,3] }))
falling.factorial <- function(a, b) gamma(b+1)*choose(a,b)
VL3exact <- ((4*n^2 - 3*n - 2)/30)/falling.factorial (10, 3) # Exact variance is from
print(c(VL3sample, VL3exact)) # Elamir and Seheult (2004, table 1, line 8)
#[1] 0.01755058 0.01703704 # the values obviously are consistent
## End(Not run)
## Not run:
# Data considered by Elamir and Seheult (2004, p. 348)
library(MASS); data(michelson); Light <- michelson$Speed
lmoms(Light, nmom=4)$lambdas # 852.4, 44.3, 0.83, 6.5 # matches those authors
lmoms.cov(Light) # [1, ] ==> 62.4267, 0.7116, 2.5912, -3.9847 # again matches
# The authors report standard error of L-kurtosis as 0.03695, which matches
lmoms.cov(Light, se="lmrse")[4] # 0.03695004
## End(Not run)
## Not run:
D <- rnorm(100) # Check results of Lmoments package.
lmoms.cov(D, rmax=5)[,5]
# lam1 lam2 lam3 lam4 lam5
#3.662721e-04 3.118812e-05 5.769509e-05 6.574662e-05 1.603578e-04
Lmoments::Lmomcov(D, rmax=5)[,5]
# L1 L2 L3 L4 L5
#3.662721e-04 3.118812e-05 5.769509e-05 6.574662e-05 1.603578e-04
## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.