lmoms.bootbarvar: The Exact Bootstrap Mean and Variance of L-moments

Description

This function computes the exact bootstrap mean and variance of L-moments using Using exact analytic expressions for the bootstrap mean and variance of any L-estimator from Hutson and Ernst (2000). The approach by those authors is to use the bootstrap distribution of the single order statistic in conjunction with the joint distribution of two water statistics. The key component is the bootstrap mean vector as well as the variance-covariance matrix of all the order statistics and then performing specific linear combinations of a basic L-estimator combined with the proportion weights used in the computation of L-moments (Lcomoment.Wk, see those examples and division by $n$). Reasonably complex algorithms are used; however, what makes those authors' contribution so interesting is that neither simulation, resampling, or numerical methods are needed.

This function provides a uniquely independent method to compute the L-moments of a sample from the vector of exact bootstrap order statistics. It is anticipated that several of the intermediate computations of this function would be of interest in further computations or graphical visualization. Therefore, this function returns many more numerical values and other L-moment function of the lmomco package. The variance-covariance matrix for large samples requires considerable CPU time; as the matrix is filled, status output is generated.

The example section of this function contains the verification of the implementation as well as provides to additional computations of variance through resampling with replacement and simulation from the parent distribution that generated the sample vector shown in the example.

Usage

lmoms.bootbarvar(x, nmom=6, verbose=TRUE)

Arguments

A vector of data values.

nmom

The number of moments to compute. Default is 6 and can not be less than 3.

verbose

A logical switch on the verbosity of the construction of the variance-covariance matrix of the order statisitics. This operation is the most time consuming of those inside the function and is provided at default of verbose=TRUE to make a gene

Value

An R list is returned.
lambdasVector of the exact bootstrap L-moments. First element is $\hat{\lambda}_1$, second element is $\hat{\lambda}_2$, and so on. This vector is from equation 1.3 and 2.4 of Hutson and Ernst (2000).
ratiosVector of the exact bootstrap L-moment ratios. Second element is $\hat{\tau}$, third element is $\hat{\tau}_3$ and so on.
lamdavarsThe exact bootstrap variances of the L-moments from equation 1.4 of Hutson and Ernst (2000) via crossprod matrix operations.
ratiovarsThe exact bootstrap variances of the L-moment ratios with NA inserted for $r=1,2$ because $r=1$ is the mean and $r=2$ for L-CV is unknown to this author.
varcovar.lambdasThe variance-covariance matrix of the L-moments from which the diagonal are the values lambdavars.
varcovar.lambdas.and.ratiosThe variance-covariance matrix of the first two L-moments and for the L-moment ratios (if nmom$>=3$) from which select diagonal are the values ratiovars.
bootstrap.orderstatisticsThe exact bootstrap estimate of the order statistics from equation 2.2 of Hutson and Ernst (2000).
varcovar.orderstatisticsThe variance-covariance matrix of the order statistics from equations 3.1 and 3.2 of Hutson and Ernst (2000). The diagonal of this matrix represents the variances of each order statistic.
inverse.varcovar.tau23The inversion of the variance-covariance matrix of $\tau_2$ and $\tau_3$ by Cholesky decomposition. This matrix may be used to estimate a joint confidence region of ($\tau_2, \tau_3$) based on asymptotic normality of L-moments.
inverse.varcovar.tau34The inversion of the variance-covariance matrix of $\tau_3$ and $\tau_4$ by Cholesky decomposition. This matrix may be used to estimate a joint confidence region of ($\tau_3, \tau_4$) based on asymptotic normality of L-moments; these two L-moment ratios likely represent the most common ratios used in general L-moment ratio diagrams.
inverse.varcovar.tau46The inversion of the variance-covariance matrix of $\tau_4$ and $\tau_6$ by Cholesky decomposition. This matrix may be used to estimate a joint confidence region of ($\tau_3, \tau_4$) based on asymptotic normality of L-moments; these two L-moment ratios represent those ratios used in L-moment ratio diagrams of symmetrical distributions.
sourceAn attribute identifying the computational source of the results: lmoms.bootbarvar.

References

Hutson, A.D., and Ernst, M.D., 2000, The exact bootstrap mean and variance of an L-estimator: Journal Royal Statistical Society B, v. 62, part 1, pp. 89--94.

Wang, D., and Hutson, A.D., 2013, Joint confidence region estimation of L-moments with an extension to survival data: Journal of Applied Statistics, in press

Examples

Run this code

para <- vec2par(c(0,1), type="gum") # Parameters of Gumbel
   n <- 10; nmom <- 6; nsim <- 2000
   # X <- rlmomco(n, para) # This is commented out because
   # the sample below is from the Gumbel distribution as in para.
   # However, the seed for the random number generator was not recorded.
   X <- c( -1.4572506, -0.7864515, -0.5226538,  0.1756959,  0.2424514,
            0.5302202,  0.5741403,  0.7708819,  1.9804254,  2.1535666)
   EXACT.BOOTLMR <- lmoms.bootbarvar(X, nmom=nmom)
   LA <- EXACT.BOOTLMR$lambdavars
   LB <- LC <- rep(NA, length(LA))
   set.seed(n)
   for(i in 1:length(LB)) {
     LB[i] <- var(replicate(nsim,
                  lmoms(sample(X, n, replace=TRUE), nmom=nmom)$lambdas[i]))
   }
   set.seed(n)
   for(i in 1:length(LC)) {
     LC[i] <- var(replicate(nsim,
                  lmoms(rlmomco(n, para), nmom=nmom)$lambdas[i]))
   }
   print(LA) # The exact bootstrap variances of the L-moments.
   print(LB) # Bootstrap variances of the L-moments by actual resampling.
   print(LC) # Simulation of the variances from the parent distribution.

   # The variances for this example are as follows:
   #> print(LA)
   #[1] 0.115295563 0.018541395 0.007922893 0.010726508 0.016459913 0.029079202
   #> print(LB)
   #[1] 0.117719198 0.018945827 0.007414461 0.010218291 0.016290100 0.028338396
   #> print(LC)
   #[1] 0.17348653 0.04113861 0.02156847 0.01443939 0.01723750 0.02512031
   # The variances, when using simulation of parent distribution,
   # appear to be generally larger than those based only on resampling
   # of the available sample of only 10 values.

   # Interested users may inspect the exact bootstrap estimates of the
   # order statistics and the variance-covariance matrix.
   # print(EXACT.BOOTLMR$bootstrap.orderstatistics)
   # print(EXACT.BOOTLMR$varcovar.orderstatistics)

   # The output for these two print functions is not shown, but what follows
   # are the numerical confirmations from A.D. Hutson (personnal commun., 2012)
   # using his personnal algorithms (outside of R).
   # Date: Jul 2012, From: ahutson, To: Asquith
   # expected values the same
   # -1.174615143125091, -0.7537760316881618, -0.3595651823632459,
   # -0.028951905838698,  0.2360931764028858,  0.4614289985084462,
   #  0.713957210869635,  1.0724040932920058,  1.5368435379648948,
   #  1.957207045977329
   # and the first two values on the first row of the matrix are
   # 0.1755400544274771,  0.1306634198810892
# Wang and Hutson (2013): Attempt to reproduce first entry of
# row 9 (n=35) in Table 1 of the reference, which is 0.878.
Xsq  <- qchisq(1-0.05, 2); n <- 35; nmom <- 4; nsim <- 1000
para <- vec2par(c(0,1), type="gum") # Parameters of Gumbel
eta  <- as.vector(lmorph(par2lmom(para))$ratios[3:4])
h <- 0
for(i in 1:nsim) {
   X <- rlmomco(n,para); message(i)
   EB <- lmoms.bootbarvar(X, nmom=nmom, verbose=FALSE)
   lmr    <- lmoms(X); etahat <- as.vector(lmr$ratios[c(3,4)])
   Pinv   <- EB$inverse.varcovar.tau34
   deta   <- (eta - etahat)
   LHS <- t(deta)   if(LHS > Xsq) { # Comparison to Chi-squared distribution
     h <- h + 1 # increment because outside ellipse
     message("Outside: ",i, " ", h, " ", round(h/i, digits=3))
   }
}
message("Empirical Coverage Probability with Alpha=0.05 is ",
        round(1 - h/nsim, digits=3), " and count is", h)
# I have run this loop and recorded an h=123 for the above
# settings. I compute a coveraged probability of 0.877,
# which agrees with Wang and Hutson (2013) within 0.001.
# Hence "very down the line" computations of
# lmoms.bootbarvar() appear to be verified.