if (FALSE) {
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
}
if (FALSE) {
# 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
# coverage 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.
}
Run the code above in your browser using DataLab