Automated Selection of r for the r Largest Order Statistics
Approach with Adjustment for Sequential Testing

arXiv:1604.01984v1 [stat.ME] 7 Apr 2016

Brian Bader · Jun Yan · Xuebin Zhang

Abstract The r largest order statistics approach is
widely used in extreme value analysis because it may
use more information from the data than just the block
maxima. In practice, the choice of r is critical. If r is
too large, bias can occur; if too small, the variance of
the estimator can be high. The limiting distribution of
the r largest order statistics, denoted by GEVr , extends
that of the block maxima. Two specification tests are
proposed to select r sequentially. The first is a score test
for the GEVr distribution. Due to the special characteristics of the GEVr distribution, the classical chi-square
asymptotics cannot be used. The simplest approach is
to use the parametric bootstrap, which is straightforward to implement but computationally expensive. An
alternative fast weighted bootstrap or multiplier procedure is developed for computational efficiency. The
second test uses the difference in estimated entropy between the GEVr and GEVr−1 models, applied to the
r largest order statistics and the r − 1 largest order
statistics, respectively. The asymptotic distribution of
the difference statistic is derived. In a large scale simulation study, both tests held their size and had substantial
power to detect various misspecification schemes. A new
This research was partially supported by an NSF grant (DMS
1521730), a University of Connecticut Research Excellence
Program grant, and a contract with Environment and Climate Change Canada.
Brian Bader
University of Connecticut
E-mail: brian.bader@uconn.edu
Jun Yan
University of Connecticut
E-mail: jun.yan@uconn.edu
Xuebin Zhang
Environment and Climate Change Canada
E-mail: xuebin.zhang@canada.ca

approach to address the issue of multiple, sequential hypotheses testing is adapted to this setting to control the
false discovery rate or familywise error rate. The utility of the procedures is demonstrated with extreme sea
level and precipitation data.
Keywords entropy · generalized extreme value ·
goodness-of-fit · multiplier bootstrap · score test ·
sequential testing

1 Introduction
The largest order statistics approach is an extension
of the block maxima approach that is often used in extreme value modeling. The focus of this paper is (Smith,
1986, p.28–29): “Suppose we are given, not just the
maximum value for each year, but the largest ten (say)
values. How might we use this data to obtain better estimates than could be made just with annual maxima?”
The r largest order statistics approach may use more
information than just the block maxima in extreme
value analysis, and is widely used in practice when
such data are available for each block. The approach
is based on the limiting distribution of the r largest
order statistics which extends the generalized extreme
value (GEV) distribution (e.g., Weissman, 1978). This
distribution, denoted by GEVr , has the same parameters as the GEV distribution, which makes it useful to
estimate the GEV parameters when the r largest values
are available for each block. The approach was investigated by Smith (1986) for the limiting joint Gumbel
distribution and extended to the more general limiting joint GEVr distribution by Tawn (1988). Because
of the potential gain in efficiency relative to the block
maxima only, the method has found many applications
such as corrosion engineering (e.g., Scarf and Laycock,

2

1996), hydrology (e.g., Dupuis, 1997), coastal engineering (e.g., Guedes Soares and Scotto, 2004), and wind
engineering (e.g., An and Pandey, 2007).
In practice, the choice of r is a critical issue in extreme value analysis with the r largest order statistics
approach. In general r needs to be small relative to the
block size B (not the number of blocks n) because as r
increases, the rate of convergence to the limiting joint
distribution decreases sharply (Smith, 1986). There is
a trade-off between the validity of the limiting result
and the amount of information required for good estimation. If r is too large, bias can occur; if too small,
the variance of the estimator can be high. Finding the
optimal r should lead to more efficient estimates of the
GEV parameters without introducing bias. A much related but different problem is the selection of threshold or fraction of a sample extreme value analysis (see
Scarrott and MacDonald, 2012, for a review). Our focus
here is the selection of r for situations where a number
of largest values are available each of n blocks. In contrast, the methods for threshold or fraction selection
reviewed in Scarrott and MacDonald (2012) deal with
a single block (n = 1) of a large size B.
The selection of r has not been as actively researched
as the threshold selection problem in the one sample
case. Smith (1986) and Tawn (1988) used probability
(also known as PP) plots for the marginal distribution
of the rth order statistic to assess its goodness of fit.
The probability plot provides a visual diagnosis, but
different viewers may reach different conclusions in the
absence of a p-value. Further, the probability plot is
only checking the marginal distribution for a specific r
as opposed to the joint distribution. Tawn (1988) suggested an alternative test of fit using a spacings results
in Weissman (1978). Let Dn:i be the spacing between
the ith and (i + 1)th largest value in a sample of size
B from a distribution in the domain of attraction of
the Gumbel distribution. Then {iDi : i = 1, . . . , r − 1}
is approximately a set of independent and identically
distributed exponential random variables as B → ∞.
The connections among the three limiting forms of the
GEV distribution (e.g., Embrechts, Klüppelberg, and
Mikosch, 1997, p.123) can be used to transform from
the Fréchet and the Weibull distribution to the Gumbel
distribution. Testing the exponentiality of the spacings
on the Gumbel scale provides an approximate diagnosis
of the joint distribution of the r largest order statistics
when B is large. A limitation of this method, however,
is that prior knowledge of the domain of attraction of
the distribution is needed. Lastly, Dupuis (1997) proposed a robust estimation method, where the weights
can be used to detect inconsistencies with the GEVr
distribution and assess the fit of the data to the joint

Brian Bader, Jun Yan, Xuebin Zhang

Gumbel model. The method can be extended to general
GEVr distributions but the construction of the estimating equations is computing intensive with Monte Carlo
integrations.
In this paper, two specification tests are proposed to
select r through a sequence of hypothesis testing. The
first is the score test (e.g., Rao, 2005), but because of
the nonstandard setting of the GEVr distribution, the
usual χ2 asymptotic distribution is invalid. A parametric bootstrap can be used to assess the significance of
the observed statistic, but is computationally demanding. A fast, large sample alternative to parametric bootstrap based on the multiplier approach (Kojadinovic
and Yan, 2012) is developed. The second test uses the
difference in estimated entropy between the GEVr and
GEVr−1 models, applied to the r largest order statistics
and the r − 1 largest order statistics, respectively. The
asymptotic distribution is derived with the central limit
theorem. Both tests are intuitive to understand, easy to
implement, and have substantial power as shown in the
simulation studies. Each of the two tests is carried out
to test the adequacy of the GEVr model for a sequence
of r values. The very recently developed stopping rules
for ordered hypotheses in G’Sell, Wager, Chouldechova,
and Tibshirani (2015) are adapted to control the false
discovery rate (FDR), the expected proportion of incorrectly rejected null hypotheses among all rejections,
or familywise error rate (FWER), the probability of at
least one type I error in the whole family of tests. All
the methods are available in an R package eva (Bader
and Yan, 2015).
The rest of the article is organized as follows. The
problem is set up in Section 2 with the GEVr distribution, observed data, and the hypothesis to be tested.
The score test is proposed in Section 3 with two implementations: parametric bootstrap and multiplier bootstrap. The entropy difference (ED) test is proposed and
the asymptotic distribution of the testing statistic is
derived in Section 4. A large scale simulation study on
the empirical size and power of the tests are reported in
Section 5. In Section 6, the multiple, sequential testing
problem is addressed by adapting recent developments
on this application. The tests are applied to sea level
and precipitation datasets in Section 7. A discussion
concludes in Section 8. The Appendix contains the details of random number generation from the GEVr distribution and a sketch of the proof of the asymptotic
distribution of the ED test statistic.
2 Model and Data Setup
The limit joint distribution of the r largest order statistics of a random sample of size B as B → ∞ is the

Automated Selection of r for the r Largest Order Statistics Approach with Adjustment for Sequential Testing

GEVr distribution with density function (Weissman,
1978)
fr (x1 , x2 , ..., xr |µ, σ, ξ)
X

r
n
o
1
1
+1
=σ −r exp − (1 + ξzr )− ξ −
log(1 + ξzj )
ξ
j=1
(1)
for some location parameter µ, scale parameter σ > 0
and shape parameter ξ, where x1 > · · · > xr , zj =
(xj − µ)/σ, and 1 + ξzj > 0 for j = 1, . . . , r. When
r = 1, this distribution is exactly the GEV distribution. The parameters θ = (µ, σ, ξ)> remain the same for
j = 1, . . . , r, r  B, but the convergence rate to the
limit distribution reduces sharply as r increases. The
conditional distribution of the rth component given the
top r − 1 variables in (1) is the GEV distribution right
truncated by xr−1 , which facilitates simulation from the
GEVr distribution; see Appendix A.
The r largest order statistics approach is an extension of the block maxima approach in extreme value
analysis when a number of largest order statistics are
available for each one of a collection of independent
blocks (Smith, 1986; Tawn, 1988). Specifically, let (xi1 ,
. . . , xir ) be the observed r largest order statistics from
block i for i = 1, . . . , n. Assuming independence across
blocks, the GEVr distribution is used in place of the
GEV distribution in the block maxima approach to
(r)
make likelihood-based inference about θ. Let li (θ) =
l(r) (xi1 , . . . , xir |θ), where
l(r) (x1 , . . . , xr |θ)
= − r log σ − (1 + ξzr )

− ξ1


−

X
r
1
+1
log(1 + ξzj )
ξ
j=1
(2)

is the contribution to the log-likelihood from a single
block (x1 , . . . , xr ). The maximum likelihood estimator
(r)
(MLE) of θ using the r largest order statistics is θ̂n =
Pn (r)
arg max i=1 li (θ).
Model checking is a necessary part of statistical analysis. The rationale of choosing a larger value of r is to
use as much information as possible, but not set r too
high so that the GEVr approximation becomes poor
due to the decrease in convergence rate. Therefore, it is
critical to test the goodness-of-fit of the GEVr distribution with a sequence of null hypotheses
(r)

H0 : the GEVr distribution fits the sample of the r
largest order statistics well
for r = 1, . . . , R, where R is the maximum, predetermined number of top order statistics to test. Two test

3

(r)

procedures for H0 are developed for a fixed r first
to help choose r ≥ 1 such that the GEVr model still
adequately describes the data. The sequential testing
process and the multiple testing issue are investigated
in Section 6.

3 Score Test
A score statistic for testing goodness-of-fit hypothesis
(r)
H0 is constructed in the usual way with the score
function and the Fisher information matrix (e.g., Rao,
2005). For ease of notation, the superscript (r) is dropped.
Define the score function
S(θ) =

n
X
i=1

Si (θ) =

n
X

∂li (θ)/∂θ

i=1

and Fisher information matrix I(θ), which have been
derived in Tawn (1988). The behaviour of the maximum likelihood estimator is the same as that derived
for the block maxima approach (Smith, 1985; Tawn,
1988), which requires ξ > −0.5. The score statistic is
Vn =

1 >
S (θ̂n )I −1 (θ̂n )S(θ̂n ).
n

Under standard regularity conditions, Vn would asymptotically follow a χ2 distribution with 3 degrees of freedom. The GEVr distribution, however, violates the regularity conditions for the score test (e.g., Casella and
Berger, 2002, pp. 516-517), as its support depends on
the parameter values unless ξ = 0. For illustration, Figure 1 presents a visual comparison of the empirical distribution of Vn with n = 5000 from 5000 replicates,
overlaid with the χ2 (3) distribution, for ξ ∈ {−0.25,
0.25} and r ∈ {1, 2, 5}. The sampling distribution of
Vn appears to be much heavier tailed than χ2 (3), and
the mismatch increases as r increases as a result of the
reduced convergence rate.
[Fig. 1 about here.]
Although the regularity conditions do not hold, the
score statistic still provides a measure of goodness-offit since it is a quadratic form of the score, which has
expectation zero under the null hypothesis. Extremely
large values of Vn relative to its sampling distribution
would suggest lack of fit, and, hence, possible misspec(r)
ification of H0 . So the key to applying the score test
is to get an approximation of the sampling distribution of Vn . Two approaches for the approximation are
proposed.

4

Brian Bader, Jun Yan, Xuebin Zhang
(k)

3.1 Parametric Bootstrap
The first solution is parametric bootstrap. For hypoth(r)
esis H0 , the test procedure goes as follows:
1. Compute θ̂n under H0 with the observed data.
2. Compute the testing statistic Vn .
3. For every k ∈ {1, ..., L} with a large number L, repeat:
(a) Generate a bootstrap sample of size n for the
r largest statistics from GEVr with parameter
vector θ̂n .
(k)
(b) Compute the θ̂n under H0 with the bootstrap
sample.
(k)
(c) Compute the score test statistic Vn .
4. Return
p-value of Vn as
PLan approximate
(k)
L−1 k=1 1(Vn > Vn ).
Straightforward as it is, the parametric bootstrap
approach involves sampling from the null distribution
and computing the MLE for each bootstrap sample,
which can be very computationally expensive. This is
especially true as the sample size n and/or the number
of order statistics r included in the model increases.

3.2 Multiplier Bootstrap
Multiplier bootstrap is a fast, large sample alternative
to parametric bootstrap in goodness-of-fit testing (e.g.,
Kojadinovic and Yan, 2012). The idea is to approximate
the asymptotic distribution of n−1/2 I −1/2 (θ)S(θ) using
its asymptotic representation
n

1 X
n−1/2 I −1/2 (θ)S(θ) = √
φi (θ),
n i=1
where φi (θ) = I −1/2 (θ)Si (θ). Its asymptotic distribution is the same as the asymptotic distribution of
n

1 X
Wn (Z, θ) = √
(Zi − Z̄)φi (θ),
n i=1
conditioning on the observed data, where Z = (Z1 , ...,
Zn ) is a set of independent and identically distributed
multipliers (independent of thePdata), with expectation
n
0 and variance 1, and Z̄ = n1 i=1 Zi . The multipliers
R∞
1
must satisfy 0 {Pr(|Z1 | > x)} 2 dx < ∞. An example
of a possible multiplier distribution is N (0, 1).
The multiplier bootstrap test procedure is summarized as follows:
1. Compute θ̂n under H0 with the observed data.
2. Compute the testing statistic Vn .
3. For every k ∈ {1, ..., L} with a large number L, repeat:

(k)

(a) Generate Z (k) = (Z1 , . . . , Zn ) from N (0, 1).
(b) Compute a realization from the approximate distribution of Wn (Z, θ) with Wn (Z (k) , θ̂n ).
(k)
(c) Compute Vn (θ̂n ) = Wn> (Z (k) , θ̂n )Wn (Z (k) , θ̂n ).
4. Return an approximate p-value of Vn as
PL
(k)
L−1 k=1 1(Vn > Vn ).
This multiplier bootstrap procedure is much faster
than parametric bootstrap procedure because, for each
sample, it only needs to generate Z and compute Wn (Z,
θ̂n ). The MLE only needs to be obtained once from the
observed data.

4 Entropy Difference Test
Another specification test for the GEVr model is derived based on the difference in entropy for the GEVr
and GEVr−1 models. The entropy for a continuous random variable with density f is (e.g., Singh, 2013)
Z

∞

E[− ln f (y)] = −

f (y) log f (y)dy.
−∞

It is essentially the expectation of negative log-likelihood.
The expectation can be approximated with the sample
average of the contribution to the log-likelihood from
the observed data, or simply the log-likelihood scaled
by the sample size n. Assuming that the r − 1 top order statistics fit the GEVr−1 distribution, the difference
in the log-likelihood between GEVr−1 and GEVr pro(r)
vides a measure of deviation from H0 . Its asymptotic
distribution can be derived. Large deviation from the
(r)
expected difference under H0 suggests a possible mis(r)
specification of H0 .
From the log-likelihood contribution in (2), the difference in log-likelihood for the ith block, Yir (θ) =
(r)
(r−1)
li − li
, is
1

1

Yir (θ) = − log σ − (1 + ξzir )− ξ + (1 + ξzir−1 )− ξ

1
−
+ 1 log(1 + ξzir ).
ξ

(3)

Pn
Pn
Let Y¯r = n1 i=1 Yir and SY2r = i=1 (Yir − Y¯r )2 /(n−1)
be the sample mean and sample variance, respectively.
Consider a standardized version of Ȳr as
Tn(r) (θ) =

√

n(Ȳr − ηr )/SYr ,

(4)

where ηr = − log σ − 1 + (1 + ξ)ψ(r), and ψ(x) =
d log Γ (x)/dx is the digamma function. The asymptotic
(r)
distribution of Tn is summarized by Theorem 1 whose
proof is relegated to Appendix B.

Automated Selection of r for the r Largest Order Statistics Approach with Adjustment for Sequential Testing
(r)

Theorem 1 Let Tn (θ) be the quantity computed based
on a random sample of size n from the GEVr distri(r−1)
bution with parameters θ and assume that H0
is
(r)
true. Then Tn converges in distribution to N (0, 1) as
n → ∞.
(r)

Note that in Theorem 1, Tn is computed from a
random sample of size n from a GEVr distribution.
If the random sample were from a distribution in the
domain of attraction of a GEV distribution, the quality
of the approximation of the GEVr distribution to the
r largest order statistics depends on the size of each
block B → ∞ with r  B. The block size B is not to
be confused with the sample size n. Assuming ξ > −0.5,
(r)
(r)
the proposed ED statistic for H0 is Tn (θ̂n ), where θ̂n
is the MLE of θ with the r largest order statistics for
the GEVr distribution. Since θ̂n is consistent for θ with
(r)
ξ > −0.5, Tn (θ̂n ) has the same limiting distribution
(r)
(r)
as Tn (θ) under H0 .
(r)
To assess the convergence of Tn (θ̂n ) to N (0, 1),
1000 GEVr replicates were simulated under configurations of r ∈ {2, 5, 10}, ξ ∈ {−0.25, 0, 0.25}, and n ∈ {50,
100}. Their quantiles are compared with those of N (0,
1) via quantile-quantile plots (not presented). It appears that a larger sample size is needed for the normal
approximation to be good for larger r and negative ξ.
This is expected because larger r means higher dimension of the data, and because the MLE only exists for
ξ > −0.5 (Smith, 1985). For r less than 5 and ξ ≥ 0,
the normal approximation is quite good; it appears satisfactory for sample size as small as 50. For r up to 10,
sample size 100 seems to be sufficient.

5 Simulation Results
5.1 Size
The empirical sizes of the tests are investigated first. For
the score test, the parametric bootstrap version and the
multiplier bootstrap version are equivalent asymptotically, but may behave differently for finite samples. It
is of interest to know how large a sample size is needed
for the two versions of the score test to hold their levels. Random samples of size n were generated from the
GEVr distribution with r ∈ {1, 2, 3, 4, 5, 10}, µ = 0,
σ = 1, and ξ ∈ {−0.25, 0, 0.25}. All three parameters
(µ, σ, ξ) were estimated.
When the sample size is small, there can be numerical difficulty in obtaining the MLE. For the multiplier bootstrap score and ED test, the MLE only needs
to obtained once, for the dataset being tested. However, in addition, the parametric bootstrap score test
must obtain a new sample and obtain the MLE for each

5

bootstrap replicate. To assess the severity of this issue,
10,000 datasets were simulated for ξ ∈ {−0.25, 0, 0.25},
r ∈ {1, 2, 3, 4, 5, 10}, n ∈ {25, 50}, and the MLE was
attempted for each dataset. Failure never occurred for
ξ ≥ 0. With ξ = −0.25 and sample size 25, the highest
failure rate of 0.69% occurred for r = 10. When the
sample size is 50, failures only occurred when r = 10,
at a rate of 0.04%.
For the parametric bootstrap score test with sample
size n ∈ {25, 50, 100}, Table 1 summarizes the empirical size of the test at nominal levels 1%, 5%, and 10%
obtained from 1000 replicates, each carried out with
bootstrap sample size L = 1000. Included only are the
cases that converged successfully. Otherwise, the results
show that the agreement between the empirical levels
and the nominal level is quite good for samples as small
as 25, which may appear in practice when long record
data is not available.
[Table 1 about here.]
For the multiplier bootstrap score test, the results
for sample sizes n ∈ {25, 50, 100, 200, 300, 400} are summarized in Table 2. When the sample size is less than
100, it appears that there is a large discrepancy between
the empirical and nominal level. For ξ ∈ {0, 0.25}, there
is reasonable agreement between the empirical level and
the nominal levels for sample size at least 100. For
ξ = −0.25 and sample size at least 100, the agreement
is good except for r = 1, in which case, the empirical
level is noticeably larger than the nominal level. This
may be due to different rates of convergence for various
ξ values as ξ moves away from −0.5. It is also interesting
to note that, everything else being held, the agreement
becomes better as r increases. This may be explained
by the more information provided by larger r for the
same sample size n, as can be seen directly in the fisher
information matrix (Tawn, 1988, pp. 247–249). For the
most difficult case with ξ = −0.25 and r = 1, the agreement gets better as sample size increases and becomes
acceptable when sample size was 1000 (not reported).
[Table 2 about here.]
(r)

To assess the convergence of Tn (θ̂n ) to N (0, 1),
10,000 replicates of the GEVr distribution were simulated with µ = 0 and σ = 1 for each configuration
of r ∈ {2, 5, 10}, ξ ∈ {−0.25, 0, 0.25}, and n ∈ {50,
100}. A rejection for nominal level α, is denoted if
(r)
|Tn (θ̂n )| > |Z α2 |, where Z α2 is the α/2 percentile of
the N(0,1) distribution. Using this result, the empirical
size of the ED test can be summarized, and the results
are presented in Table 3.
[Table 3 about here.]

6

Brian Bader, Jun Yan, Xuebin Zhang

For sample size 50, the empirical size is above the
nominal level for all configurations of r and ξ. As the
sample size increases from 50 to 100, the empirical size
stays the same or decreases in every setting. For sample
size 100, the agreement between nominal and observed
size appears to be satisfactory for all configurations of
r and ξ. For sample size 50, the empirical size is slightly
higher than the nominal size, but may be acceptable to
some practitioners. For example, the empirical size for
nominal size 10% is never above 12%, and for nominal
size 5%, empirical size is never above 7%.
In summary, the multiplier bootstrap procedure of
the score test can be used as a fast, reliable alternative to the parametric bootstrap procedure for sample
size 100 or more when ξ ≥ 0. When only small samples
are available (less than 50 observations), the parametric bootstrap procedure is most appropriate since the
multiplier version does not hold its size and the ED test
relies upon samples of size 50 or more for the central
limit theorem to take effect.

5.2 Power
The powers of the score tests and the ED test are studied with two data generating schemes under the alternative hypothesis. In the first scheme, 4 largest order
statistics were generated from the GEV4 distribution
with µ = 0, σ = 1, and ξ ∈ {−0.25, 0, 0.25}, and
the 5th one was generated from a KumGEV distribution right truncated by the 4th largest order statistic.
The KumGEV distribution is a generalization of the
GEV distribution (Eljabri, 2013) with two additional
parameters a and b which alter skewness and kurtosis. Defining Gr (x) to be the distribution function of
the GEVr (µ, σ, ξ) distribution, the distribution function of the KumGEVr (µ, σ, ξ, a, b) is given by Fr (x) =
1 − {1 − [Gr (x)]a }b for a > 0, b > 0. The score test and
the ED test were applied to the top 5 order statistics
with sample size n ∈ {100, 200}. When a = b = 1, the
null hypothesis of GEV5 is true. Larger difference from
1 of parameters a and b means larger deviation from
the null hypothesis of GEV5 .
[Table 4 about here.]
Table 4 summarizes the empirical rejection percentages obtained with nominal size 5%, for a sequence
value of a = b from 0.4 to 2.0, with increment 0.2. Both
tests hold their sizes when a = b = 1 and have substantial power in rejecting the null hypothesis for other values of a = b. Between the two tests, the ED test demonstrated much higher power than the score test in the
more difficult cases where the deviation from the null

hypothesis is small; for example, the ED test’s power
almost doubled the score test’s power for a = b ∈ {0.8,
1.2}. As expected, the powers of both tests increase as
a = b moves away from 1 or as the sample sizes increases.
In the second scheme, top 6 order statistics were
generated from the GEV6 distribution with µ = 0,
σ = 1, and ξ ∈ {−0.25, 0, 0.25}, and then the 5th order
statistic was replaced from a mixture of the 5th and 6th
order statistics. The tests were applied to the sample of
first 5 order statistics with sample sizes n ∈ {100, 200}.
The mixing rate p of the 5th order statistic took values
in {0.00, 0.10, 0.25, 0.50, 0.75, 0.90, 1.00}. When p = 1
the null hypothesis of GEV5 is true. Smaller values of
p indicate larger deviations from the null. Again, both
tests hold their sizes when p = 1 and have substantial
power for other values of p, which increases as p decreases or as the sample sizes increases. The ED test
again outperforms the score test with almost doubled
power in the most difficult cases with p ∈ {0.75, 0.90}.
For sample size 100 with p = 0.50, for instance, the ED
test has power above 93% while the score test only has
power above 69%.
[Table 5 about here.]

6 Automated Sequential Testing Procedure
(r)

As there are R hypotheses H0 , r = 1, . . . , R, to be
tested in a sequence in the methods proposed, the sequential, multiple testing issue needs to be addressed.
Most methods for error control assume that all the tests
can be run first and then a subset of tests are chosen
to be rejected (e.g., Benjamini, 2010a,b). The errors to
be controlled are either the FWER (Shaffer, 1995), or
the FDR (Benjamini and Hochberg, 1995; Benjamini
and Yekutieli, 2001). In contrast to the usual multiple
testing procedures, however, a unique feature in this
setting is that the hypotheses must be rejected in an
(r)
(k)
ordered fashion: if H0 is rejected, r < R, then H0
will be rejected for all r < k ≤ R. Despite the extensive
literature on multiple testing and the more recent developments on FDR control and its variants, no definitive
procedure has been available for error control in ordered
tests until the recent work of G’Sell et al (2015).
Consider a sequence of null hypotheses H1 , . . . , Hm .
An ordered test procedure must reject H1 , . . . , Hk for
some k ∈ {0, 1, . . . , m}, which rules out the classical
methods for FDR control (Benjamini and Hochberg,
1995). Let p1 , . . . , pm ∈ [0, 1] be the corresponding pvalues of the m hypotheses such that pj is uniformly
distributed over [0, 1] when Hj is true. The methods of
G’Sell et al (2015) transform the sequence of p-values

Automated Selection of r for the r Largest Order Statistics Approach with Adjustment for Sequential Testing

to a monotone sequence and then apply the original
Benjamini–Hochberg procedure on the monotone sequence. They proposed two rejections rules, each returning a cutoff k̂ such that H1 , . . . , Hk̂ are rejected.
The first is called ForwardStop,
(
)
k
1X
k̂F = max k ∈ {1, . . . , m} : −
log(1 − pi ) ≤ α ,
k i=1
and the second is called StrongStop,


m

X
log pj  αk 
≤
,
k̂S = max k ∈ {1, . . . , m} : exp

j
m
j=k

where α is a pre-specified level. Both rules were shown
to control the FDR at level α under the assumption of
independent p-values. ForwardStop sets the rejection
threshold at the largest k at which the average of first
k transformed p-values is small enough. As it does not
depend on those p-values with later indices, this rule
is robust to potential misspecification at later indices.
StrongStop offers a stronger guarantee than ForwardStop. If the non-null p-values indeed precede the null
p-values, it controls the FWER at level α in addition
to the FDR. Thus, for ForwardStop, this α refers to the
FDR and for StrongStop, α refers to the FWER. As the
decision to stop at k depends on all the p-values after
k, its power may be harmed if, for example, the very
last p-values are slightly higher than expected under
the null hypotheses.
To apply the two rules to our setting, note that our
objective is to give a threshold r̂ such that the first r̂
of m = R hypotheses are accepted instead of rejected.
Therefore, we put the p-values in reverse order: let the
ordered set of p-values {p1 , . . . , pR } correspond to hy(R)
(1)
potheses {H0 , . . . , H0 }. The two rules give a cut(R)
off k̂ ∈ {1, . . . , R} such that the hypotheses H0 , . . . ,
(R−k̂+1)

H0
are rejected. If no k̂ ∈ {1, . . . , R} exists, then
no rejection is made.
A caveat is that, unlike the setting of G’Sell et al
(2015), the p-values of the sequential tests are dependent. Nonetheless, the ForwardStop and StrongStop procedures may still provide some error control. For example, in the non-sequential multiple testing scenario Benjamini and Yekutieli (2001) show that their procedure
controls the FDR under certain positive dependency
conditions, while Blanchard and Roquain (2009) implement adaptive versions of step-up procedures that provably control the FDR under unspecified dependence
among p-values.
The empirical properties of the two rules for the
tests in this paper are investigated in simulation studies. To check the empirical FWER of the StrongStop

7

rule, only data under the null hypotheses are needed.
With R = 10, ξ ∈ {−0.25, 0.25}, n ∈ {30, 50, 100, 200},
µ = 0, and σ = 1, 1000 GEV10 samples were generated.
For the ED, multiplier bootstrap score, and parametric
bootstrap score test, the observed FWER is compared
to the expected rates at various nominal α control levels. The StrongStop procedure is used, as well as no
error control (i.e. a rejection occurs any time the raw
p-value is below the nominal level). The results of this
simulation are presented in Figure 2.
It is clear that the StrongStop reasonably controls
the FWER for the ED test and the agreement between
the observed and expected rate increases as the sample
size increases. For both the parametric and multiplier
bootstrap versions of the score test however, the observed FWER is above the expected rate, at times 10%
higher. Regardless, it is apparent that using no error
control results in an inflated FWER, and this inflation
can only increase as the number of tests increase.
[Fig. 2 about here.]
To check the empirical FDR of the ForwardStop
rule, data need to be generated from a non-null model.
To achieve this, consider the sequence of specification
tests of GEVr distribution with r ∈ {1, . . . , 6}, where
the 5th and 6th order statistics are misspecified. Specifically, data from the GEV7 distribution with µ = 0 and
σ = 1 were generated for n blocks; then the 5th order
statistic is replaced with a 50/50 mixture of the 5th
and 6th order statistics, and the 6th order statistic is
replaced with a 50/50 mixture of the 6th and 7th order
statistics. This is replicated 1000 times for each value of
ξ ∈ {−0.25, 0.25} and n ∈ {30, 50, 100, 200}. For nominal level α, the observed FDR is defined as the number
of false rejections (i.e. any rejection of r ≤ 4) divided
by the number of total rejections.
The results are presented in Figure 3. The plots
show that the ForwardStop procedure controls the FDR
for the ED test, while for both versions of the score test,
the observed FDR is slightly higher than the expected
at most nominal rates. Here, sample size does not appear to effect the observed rates.
[Fig. 3 about here.]
Similarly, the observed FWER rate in this particular simulation setting can be found by taking the number of simulations with at least one false rejection (here,
any rejection of r ≤ 4) and dividing that number by
the total number of simulations. This calculation is
performed for a variety of nominal levels α, using the
StrongStop procedure. The results are presented in Figure 3. In this particular simulation setting, the StrongStop

8

Brian Bader, Jun Yan, Xuebin Zhang

procedure controls the FWER for the ED test and both
versions of the score test at all sample sizes investigated.
It is of interest to investigate the performance of
the ForwardStop and StrongStop in selecting r for the
r largest order statistics method. In the simulation setting described in the last paragraph, the correct choice
of r should be 4, and a good testing procedure should
provide a selection r̂ close to 4. The choice r̂ ∈ {0,
. . . , 6} is recorded using the ED test and bootstrap
score tests with both ForwardStop and StrongStop. Due
to space constraints, we choose to present one setting,
where ξ = 0.25 and n = 100. The non-adjusted sequential procedure is also included, testing in an ascending
manner from r = 1 and r̂ is chosen by the first rejection
found (if any). The results are summarized in Table 6.
[Table 6 about here.]
In general, larger choices of α lead to a higher percentage of r̂ = 4 being correctly chosen with ForwardStop or StrongStop. Intuitively, this is not surprising
since a smaller α makes it more difficult to reject the
‘bad’ hypotheses of r ∈ {5, 6}. A larger choice of α
also leads to a higher probability of rejecting too many
tests; i.e. choosing r too small. From the perspective of
model specification, this is more desirable than accepting true negatives. A choice of 6, 5, or 0 is problematic,
but choosing 1, 2, or 3 is acceptable, although some information is lost. When no adjustment is used and an
ascending sequential procedure is used, both tests have
reasonable classification rates. When α = 0.05, the ED
test achieves the correct choice of r 79.9% of the time,
with the parametric bootstrap and multiplier bootstrap
score tests achieving 68.1% and 59.6% respectively. Of
course, as the number of tests (i.e., R) increase, with
no adjustment the correct classification rates will go
down and the ForwardStop/StrongStop procedures will
achieve better rates. This may not be too big an issue here as R is typically small. In the case where rich
data are available and R is big, the ForwardStop and
StrongStop becomes more useful as they are designed
to handle a large number of ordered hypothesis.

7 Illustrations
7.1 Lowestoft Sea Levels
Sea level readings in 60 and 15 minute intervals from a
gauge at Lowestoft off the east coast of Britain during
the years 1964–2014 are available from the UK Tide
Gauge Network website. The readings are hourly from
1964–1992 and in fifteen minute intervals from 1993 to
present. Accurate estimates of extreme sea levels are of

great interest. The current data are of better quality
and with longer record than those used in Tawn (1988)
— annual maxima during 1953–1983 and hourly data
during 1970–78 and 1980–82.
Justification of the statistical model was considered
in detail by Tawn (1988). The three main assumptions
needed to justify use of the GEVr model are: (1) The
block size B is large compared to the choice of r; (2)
Observations within each block and across blocks are
approximately independent; and (3) The distribution
of the block maxima follows GEV1 . The first assumption is satisfied, by letting R = 125, and noting that
the block size for each year is B = 365 × 24 = 8760
from 1964–1992 and B = 365 × 96 = 35040 from 1993–
2014. This ensures that r  B. The third assumption
is implicitly addressed in the testing procedure; if the
goodness-of-fit test for the block maxima rejects, all
subsequent tests for r > 1 are rejected as well.
The second assumption can be addressed in this
setting by the concept of independent storms (Tawn,
1988). The idea is to consider each storm as a separate
event, with each storm having some storm length, say
τ . Thus, when selecting the r largest values from each
block, only a single contribution can be obtained from
each storm, which can be considered the r largest independent annual events. By choosing τ large enough, this
ensures both approximate independence of observations
within each block and across blocks. The procedure to
extract the independent r largest annual events is as
follows:
1. Pick out the largest remaining value from the year
(block) of interest.
2. Remove observations within a lag of τ /2 from both
sides of the value chosen in step 1.
3. Repeat (within each year) until the r largest are
extracted.
A full analysis is performed on the Lowestoft sea
level data using τ = 60 as the estimated storm length
(Tawn, 1988). Using R = 125, both the parametric
bootstrap score (with bootstrap sample size L = 10,
000) and ED test are applied sequentially on the data.
The p-values of the sequential tests (adjusted and unadjusted) can be seen in Figure 4. Due to the large number
of tests, the adjustment for multiplicity is desired and
thus, ForwardStop is used to choose r. For this dataset,
the score test is more powerful than the ED test. With
ForwardStop and the score test, Figure 4 suggests that
r = 33. The remainder of this analysis proceeds with
the choice of r = 33. The estimated parameters and
corresponding 95% profile confidence intervals for r = 1
through r = 40 are shown in Figure 5.
[Fig. 4 about here.]

Automated Selection of r for the r Largest Order Statistics Approach with Adjustment for Sequential Testing

[Fig. 5 about here.]
When r = 33, the parameters are estimated as µ̂ =
3.462 (0.023), σ̂ = 0.210 (0.013), and ξˆ = −0.017 (0.023),
with standard errors in parenthesis. An important risk
measure is the t-year return level zt (e.g., Hosking, 1990;
Ribereau, Guillou, and Naveau, 2008; Singo, Kundu,
Odiyo, Mathivha, and Nkuna, 2012). It can be thought
of here as the sea level that is exceeded once every t
years on average. Specifically, the t-year return level is
the 1 − 1/t quantile of the GEV distribution
(

µ − σξ 1 − [− log(1 − 1t )]−ξ , ξ 6= 0,
zt =
µ − σ log[− log(1 − 1t )],
ξ = 0.

9

are applied sequentially on the data. The p-values of the
sequential tests (ForwardStop, StrongStop, and unadjusted) are shown in Figure 6. The score test does not
pick up anything. The ED test obtains p-values 0.002
and 0.016, respectively, for r = 9 and r = 10, which
translates into a rejection using ForwardStop. Thus,
Figure 6 suggests that r = 8 be used for the analysis.
[Fig. 7 about here.]

With r = 8, the estimated parameters are given
as µ̂ = 6.118 (0.139), σ̂ = 2.031 (0.118), and ξˆ =
0.219 (0.032). This suggests a heavy upper tail for the
estimated distribution (i.e. ξˆ > 0). The progression of
parameters and certain return level estimates can be
The return levels can be estimated with parameter valseen in Figure 7. The 50, 100, and 200 year return level
ues replaced with their estimates, and confidence inter95% confidence intervals for r = 8 are calculated using
vals can be constructed using profile likelihood (e.g.,
the profile likelihood method and are given by (16.019,
Coles, 2001, p.57).
22.411), (18.606, 26.979), and (21.489, 31.136), respecThe 95% profile likelihood confidence intervals for
tively. The advantages of using r = 8 versus the block
the 50, 100, and 200 year return levels (i.e. z50 , z100 ,
maxima for analysis are quite clear from Figure 7. The
z200 ) are given by (4.102, 4.461), (4.210, 4.641) and (4.312,
standard error of the shape parameter decreases from
4.824), respectively. The benefit of using r = 1 versus
0.071 to 0.032, a decrease of over 50%. Similarly, the 50
r = 33 can be seen in the return level confidence interyear return level 95% confidence intervals decreases in
vals in Figure 5. For example, the point estimate of the
width by over 25%.
100 year return level decreases slightly as r increases
and the width of the 95% confidence interval decreases
drastically from 2.061 (r = 1) to 0.432 (r = 33), as
8 Discussion
more information is used. The lower bound of the interval however remains quite stable, shifting from 4.330
We proposed two model specification tests for a fixed
to 4.210 — less than a 3% change. Similarly, the stannumber of largest order statistics as the basis for sedard error of the shape parameter estimate decreases
lecting r for the r largest order statistics approach in
by over two-thirds when using r = 33 versus r = 1.
extreme value analysis. The score test has two ver7.2 Annual Maximum Precipitation: Atlantic City, NJ
The top 10 annual precipitation events (in centimeters)
were taken from the daily records of a rain gauge station
in Atlantic City, NJ from 1874–2015. The year 1989 is
missing, while the remaining records are greater than
98% complete. This provides a total record length of 141
years. The raw data is a part of the Global Historical
Climatology Network (GHCN-Daily), with an overview
given by Menne, Durre, Vose, Gleason, and Houston
(2012). The specific station identification in the dataset
is USW00013724.
[Fig. 6 about here.]
Unlike for the Lowestoft sea level data, a rather
small value is set for R at R = 10 because of the much
lower frequency of the daily data. Borrowing ideas from
Section 7.1, a storm length of τ = 2 is used to ensure approximate independence of observations. Both the parametric bootstrap score (with L = 10, 000) and ED test

sions of bootstrap procedure: the multiplier bootstrap
method providing a fast, large sample alternative to the
parametric bootstrap method, with a speedup of over
100 times. The ED test depends on asymptotic normal approximation of the testing statistic, which becomes acceptable for sample size over 50. It assumes
that the r − 1 top order statistics included already fits
the GEVr−1 distribution. Therefore, the initial hypothesis at r = 1 needs to be tested with the score tests.
Both tests hold their size better when the shape parameter is further away from the lower limit of −0.5
or sample size is larger. When only small samples are
available (50 observations or less), the parametric bootstrap score test is recommended.
Alternative versions of the ED test have been explored. One may define the testing statistics as the difference in entropy between GEV1 and GEVr , instead of
between GEVr −1 and GEVr . Nonetheless, it appeared
to require a larger sample to hold its size from our simu(r)
lation studies (not reported). In the calculation of Tn ,
(1)
the block maxima MLE θ̂n can be used as an estimate

10

Brian Bader, Jun Yan, Xuebin Zhang
(r)

for θ in place of θ̂n . Again, in our simulation studies,
this version of the ED test was too conservative, thus
reducing the power, when the sample size was not large
enough. This may be explained in that the resulting
ŜYr underestimates SYr .
Naively, the tests may be performed sequentially for
each r ∈ {1, . . . , R}, for a prefixed, usually small R, at
(r)
a certain significance level until H0 is rejected. The issue of multiple, sequential testing is addressed in detail
by adapting two very recent stopping rules to control
the FDR and the FWER that are developed specifically
for situations when hypotheses must be rejected in an
ordered fashion (G’Sell et al, 2015). It is shown that
these automated procedures reasonably control the error rate for the tests discussed in this paper. The naive
unadjusted sequential testing procedure also appears to
have decent performance at choosing the correct r and
for this scenario is more conservative (selects a smaller
r) than the stopping rules.
The tests can be extended to allow covariates in
the parameters. For example, extremal precipitation in
a year may be affected by large scale climate indexes
such as the Southern Oscillation Index (SOI), which
may be incorporated as a covariate in the location parameter (e.g., Shang, Yan, and Zhang, 2011). Both tests
can be carried out with additional model parameters.
When the underlying data falls into a rich class of dependence structures (such as time series), this dependence may be incorporated directly instead of using a
procedure to achieve approximate independence (e.g.
the storm length τ in Section 7). For example, take
the GEV-GARCH model (Zhao, Scarrott, Oxley, and
Reale, 2011) when r = 1. It may be extended to the
case where r > 1 and the tests presented here may be
applied to select r under this model assumption.

The resulting vector (x1 , . . . , xr ) is a single observation from
the GEVr distribution.
For ξ → 0, caveat is needed in numerical evaluation. Using function expm1 for exp(1 + x) for x → 0 provides much
improved accuracy in comparison to a few implementations
in existing R packages. For readability, here is a simplified
version of our implemetation in R package eva (Bader and
Yan, 2015).
## Quantile function of a GEVr(loc, scale, shape)
qgev <- function(p, loc = 0, scale = 1, shape = 0,
lower.tail = TRUE, log.p = FALSE) {
if (log.p) p <- exp(p)
if(shape == 0) {
loc - scale * log(-log(p))
} else
loc + scale * expm1(log(-log(p)) * -shape) / shape)
}

## Random number generator of GEVr;
## Returns a matrix of n rows and r columns,
## each row a draw from GEVr
rgevr <- function(n, r, loc = 0, scale = 1, shape = 0) {
umat <- matrix(runif(n * r), n, r)
if (r > 1) {
matrix(qgev(t(apply(umat, 1, cumprod)),
loc, scale, shape),
ncol = r)
} else {
qgev(umat, loc, scale, shape)
}
}

B Asymptotic Distribution of Tn(r) (θ)
Proof (Theorem 1) Consider a random vector (X1 , X2 , ..., Xr )
which follows a GEVr (θ) distribution. The following result
given by Tawn (1988, pg. 248) will be used:

A Generating from GEVr Distribution
The GEVr distribution is closely connected to the GEV distribution. Let X1 > · · · > Xr follow a GEVr distribution (1).
It is obvious that the GEV1 distribution is the GEV distribution with the same parameters, which is the marginal
distribution of X1 . More interestingly, note that, the conditional distribution of X2 given X1 = x1 is simply the GEV
distribution righted truncated by x1 . In general, given (X1 ,
. . . , Xk ) = (x1 , . . . , xk ) for 1 ≤ k < r, the conditional distribution of Xk+1 is the GEV distribution righted truncated at xk .
This property can be exploited to generate the r components
in a realized GEVr observation.
The pseudo algorithm to generate a single observation is
the following:
– Generate the first value x1 from the (unconditional) GEV
distribution.
– For i = 2, . . . , r:
– Generate xi from the GEV distribution right truncated by xi−1 .

=

(−ξ )c−a
Γ (j )

1
−( ξ
+ b)

logc (1 + ξZj )]
a
a
X
(−1)α
Γ (c) (j + bξ − αξ + 1)

h(j|θ, a, b, c) ≡ E [Zja (1 + ξZj )

α=0

α

(5)
where Zj = (Xj − µ)/σ and Γ (c) is the cth derivative of the
gamma function, for a ∈ Z, b ∈ R, and c ∈ Z, such that
(j + bξ − αξ + 1) 6∈ {0, −1, −2, . . .}, α = 0, 1, . . . , a.
Assume that ξ 6= 0 and 1 + ξZj > 0 for j = 1, . . . , r.
The difference in log-likelihoods for a single observation from
the GEVr (θ) and GEVr−1 (θ) distribution is given by (3) in
Section 4. Thus, the first moment of Yir is
E [Y1r ] = − log σ − h(r|θ, 0, 0, 0) + h(r − 1|θ, 0, 0, 0)


−

1
ξ


+ 1 h(r|θ, 0, −ξ −1 , 1)

= − log σ − 1 + (1 + ξ )ψ (r)
where ψ (x) =

Γ (1) (x)
.
Γ (x)

Automated Selection of r for the r Largest Order Statistics Approach with Adjustment for Sequential Testing
To prove that the second moment of Yir is finite, note
that
(
1
−ξ

|Y1r | ≤ 4 max

log σ , (1 + ξZ1r )
1
−ξ

(1 + ξZ1r−1 )

,

1
ξ

,

)



+ 1 log(1 + ξZ1r−1 )

,

which implies
(
Y12r

≤ 16

max

1
−ξ

log σ , (1 + ξZ1r )
1
−ξ

(1 + ξZ1r−1 )

,

1
ξ

,



+ 1 log(1 + ξZ1r−1 )

)!2
.

The bound of E (Y12r ) can be established by applying (5) to
the last three terms in the max operator,
2
−ξ

E [(1 + ξZ1r )

2

−ξ

E [(1 + ξZ1r−1 )

] = h(r|θ, 0, ξ −1 , 0) < ∞,
] = h(r − 1|θ, 0, ξ −1 , 0) < ∞,

E [log2 (1 + ξZ1r−1 )] = h(r − 1|θ, 0, −ξ −1 , 2) < ∞.

The desired result then follows from the central limit theorem
and Slutsky’s theorem.
The case where ξ = 0 in Theorem 1 can easily be derived by taking the limit as ξ → 0 in (3) and in (5) by the
Dominated Convergence Theorem.
Acknowledgements We would like to thank Dr. Zhiyi Chi
for his comments and discussion on the multiple testing issues.

References
An Y, Pandey MD (2007) The r largest order statistics model
for extreme wind speed estimation. Journal of Wind Engineering and Industrial Aerodynamics 95(3):165–182
Bader B, Yan J (2015) eva: Extreme Value Analysis with
Goodness-of-Fit Testing
Benjamini Y (2010a) Discovering the false discovery rate.
Journal of the Royal Statistical Society: Series B (Statistical Methodology) 72(4):405–416
Benjamini Y (2010b) Simultaneous and selective inference:
Current successes and future challenges. Biometrical Journal 52(6):708–721
Benjamini Y, Hochberg Y (1995) Controlling the false discovery rate: A practical and powerful approach to multiple
testing. Journal of the Royal Statistical Society Series B
57(1):289–300
Benjamini Y, Yekutieli D (2001) The control of the false discovery rate in multiple testing under dependency. The Annals of Statistics 29(4):1165–1188
Blanchard G, Roquain É (2009) Adaptive false discovery rate
control under independence and dependence. The Journal
of Machine Learning Research 10:2837–2871
Casella G, Berger RL (2002) Statistical Inference, 2nd edn.
Duxbury Pacific Grove, CA
Coles S (2001) An Introduction to Statistical Modeling of
Extreme Values, 1st edn. Springer
Dupuis DJ (1997) Extreme value theory based on the r largest
annual events: A robust approach. Journal of Hydrology
200(1):295–306

11

Eljabri SSM (2013) New statistical models for extreme values. PhD thesis, The University of Manchester, Manchester, UK
Embrechts P, Klüppelberg C, Mikosch T (1997) Modelling
Extremal Events, vol 33. Springer Science & Business Media
G’Sell MG, Wager S, Chouldechova A, Tibshirani R (2015)
Sequential selection procedures and false discovery rate
control. Journal of the Royal Statistical Society: Series B
(Statistical Methodology) DOI 10.1111/rssb.12122, URL
http://dx.doi.org/10.1111/rssb.12122, forthcoming
Guedes Soares C, Scotto MG (2004) Application of the r
largest-order statistics for long-term predictions of significant wave height. Coastal Engineering 51(5):387–394
Hosking JRM (1990) L-moments: Analysis and estimation of
distributions using linear combinations of order statistics.
Journal of the Royal Statistical Society Series B (Methodological) 52(1):105–124
Kojadinovic I, Yan J (2012) Goodness-of-fit testing based on
a weighted bootstrap: A fast large-sample alternative to
the parametric bootstrap. Canadian Journal of Statistics
40(3):480–500
Menne MJ, Durre I, Vose RS, Gleason BE, Houston TG
(2012) An overview of the global historical climatology network-daily database. Journal of Atmospheric and
Oceanic Technology 29(7):897–910
Rao CR (2005) Score test: Historical review and recent developments. In: Balakrishnan N, Kannan N, Nagaraja HN
(eds) Advances in Ranking and Selection, Multiple Comparisons, and Reliability, Springer, pp 3–20
Ribereau P, Guillou A, Naveau P (2008) Estimating return
levels from maxima of non-stationary random sequences
using the generalized PWM method. Nonlinear Processes
in Geophysics 15(6):1033–1039
Scarf PA, Laycock PJ (1996) Estimation of extremes in corrosion engineering. Journal of applied statistics 23(6):621–
644
Scarrott C, MacDonald A (2012) A review of extreme
value threshold estimation and uncertainty quantification.
REVSTAT–Statistical Journal 10(1):33–60
Shaffer JP (1995) Multiple hypothesis testing. Annual Review
of Psychology 46(1):561–584
Shang H, Yan J, Zhang X (2011) El Nin̈o–Southern Oscillation influence on winter maximum daily precipitation in
California in a spatial model. Water Resources Research
47:W11,507–W11,515, DOI 10.1029/2011WR010415
Singh VP (2013) Entropy Theory and Its Application in Environmental and Water Engineering. John Wiley & Sons
Singo LR, Kundu PM, Odiyo JO, Mathivha FI, Nkuna
TR (2012) Flood frequency analysis of annual maximum stream flows for Luvuvhu river catchment, Limpopo
Province, South Africa. Tech. rep., University of Venda,
Department of Hydrology and Water Resources
Smith RL (1985) Maximum likelihood estimation in a class
of nonregular cases. Biometrika 72(1):67–90
Smith RL (1986) Extreme value theory based on the r largest
annual events. Journal of Hydrology 86(1):27–43
Tawn JA (1988) An extreme-value theory model for dependent observations. Journal of Hydrology 101(1):227–250
Weissman I (1978) Estimation of parameters and large quantiles based on the k largest observations. Journal of the
American Statistical Association 73(364):812–815
Zhao X, Scarrott CJ, Oxley L, Reale M (2011) GARCH dependence in extreme value models with Bayesian inference.
Mathematics and Computers in Simulation 81(7):1430–
1440

12

TABLES

Table 1: Empirical size (in %) for the parametric bootstrap score test under the null distribution GEVr , with µ = 0
and σ = 1 based on 1000 samples, each with bootstrap sample size L = 1000.
Sample Size

r

Nominal Size

25

50

100

1.0

5.0

10.0

1.0

5.0

10.0

1.0

5.0

10.0

ξ = −0.25

1
2
3
4
5
10

0.4
0.1
0.3
0.3
0.4
2.7

2.8
2.6
2.5
1.8
2.4
5.3

6.0
6.0
5.0
5.4
6.7
8.7

1.1
0.8
0.8
0.6
0.4
0.5

4.8
3.4
4.3
3.1
3.3
3.9

9.3
6.5
7.7
6.9
8.3
8.4

0.6
0.6
1.1
1.1
0.6
0.7

4.1
3.6
4.8
5.1
3.1
4.2

8.0
8.1
8.1
8.8
6.5
7.6

ξ=0

1
2
3
4
5
10

1.3
1.4
1.7
1.5
1.6
1.5

5.2
5.1
6.2
4.5
5.8
4.0

8.9
9.4
10.9
8.5
10.4
7.3

1.6
2.0
2.1
1.3
2.4
1.5

5.3
4.9
6.0
6.0
6.2
4.3

9.0
10.0
10.2
10.2
9.9
8.9

0.8
1.0
0.8
1.0
1.2
0.7

4.7
4.3
4.9
4.4
5.0
4.6

9.3
9.9
9.8
9.8
9.7
8.2

ξ = 0.25

1
2
3
4
5
10

1.7
1.8
1.5
1.2
1.7
1.1

4.5
5.1
4.4
3.3
4.4
4.6

9.7
8.7
9.4
8.1
9.4
8.3

2.6
1.8
1.5
1.1
1.1
1.5

7.1
4.4
3.7
4.6
4.2
6.1

11.5
8.5
8.1
9.7
8.6
10.7

1.1
0.5
1.0
1.1
0.6
1.0

4.6
2.9
4.2
4.3
4.8
3.9

9.1
7.5
9.4
9.6
9.6
8.5

TABLES

13

Table 2: Empirical size (in %) for multiplier bootstrap score test under the null distribution GEVr , with µ = 0 and
σ = 1. 1000 samples, each with bootstrap sample size L = 1000 were used. Although not shown, the empirical size
for r = 1 and ξ = −0.25 becomes acceptable when sample size is 1000.
Sample Size

r

Nominal Size
ξ = −0.25

1
2
3
4
5
10
1
2
3
4
5
10
1
2
3
4
5
10

ξ=0

ξ = 0.25

Sample Size

ξ=0

ξ = 0.25

1
2
3
4
5
10
1
2
3
4
5
10
1
2
3
4
5
10

50

100

1.0

5.0

10.0

1.0

5.0

10.0

1.0

5.0

10.0

7.0
2.0
2.1
3.3
3.6
2.0
3.3
2.8
6.1
5.1
4.2
3.1
1.8
5.7
7.1
5.4
4.4
3.3

13.4
6.9
5.8
7.2
9.0
7.0
8.4
8.7
12.1
10.4
9.0
9.2
6.7
12.7
12.2
9.8
10.1
8.9

18.9
13.4
11.7
12.3
14.0
10.3
15.3
14.4
16.5
14.5
14.5
14.4
13.7
17.1
16.5
16.8
15.8
15.3

6.3
1.3
1.1
1.1
2.3
2.6
2.2
1.8
3.0
3.6
2.2
2.4
1.3
4.7
5.3
3.7
3.5
2.4

13.8
6.4
5.9
4.9
6.8
7.4
7.0
7.5
7.2
10.1
8.2
6.4
4.7
9.9
9.4
9.0
8.2
6.6

19.6
12.4
11.1
10.8
11.2
12.8
12.5
13.0
12.2
14.9
12.5
9.8
10.4
14.9
14.8
13.4
13.6
12.3

5.4
1.6
1.1
1.0
1.1
2.1
1.1
0.9
1.5
1.0
1.7
0.6
0.8
3.5
4.2
2.6
2.4
1.6

11.4
6.9
5.0
5.2
6.2
6.4
4.6
5.7
6.0
5.6
6.5
4.6
4.4
7.4
8.4
6.0
7.4
5.8

16.3
13.6
10.8
11.9
10.6
10.1
9.2
10.3
10.4
10.3
12.0
9.0
11.5
11.6
12.5
11.4
11.4
10.9

r

Nominal Size
ξ = −0.25

25

200

300

400

1.0

5.0

10.0

1.0

5.0

10.0

1.0

5.0

10.0

5.4
1.4
1.5
1.1
1.1
1.4
1.3
0.5
1.4
1.1
1.8
1.1
0.9
3.2
1.8
1.2
1.6
1.7

10.5
6.7
5.9
5.6
4.5
6.4
6.1
5.0
4.5
5.4
6.2
3.8
4.9
7.9
6.1
4.9
5.9
6.6

15.2
12.8
11.8
10.6
9.3
11.6
11.2
10.6
9.8
10.6
12.5
9.3
11.4
11.7
10.7
11.2
10.0
12.4

3.6
1.4
1.1
1.0
1.2
1.7
0.8
1.2
1.3
1.3
0.9
0.9
0.9
2.3
2.6
1.2
1.3
0.9

8.2
6.4
5.4
5.6
5.7
6.2
5.2
5.7
6.0
5.2
4.6
5.2
5.0
7.1
7.0
6.0
7.3
4.4

12.6
11.4
10.8
11.5
11.7
11.3
10.0
11.8
9.6
9.9
9.8
12.6
10.8
11.2
11.2
9.9
11.8
9.8

2.8
1.4
1.2
1.0
1.2
0.8
1.0
1.0
0.9
0.9
1.2
1.2
0.7
2.5
1.0
1.2
1.2
1.6

7.1
5.1
6.6
4.7
4.7
5.0
5.1
5.9
4.4
5.0
4.6
4.9
5.2
6.6
4.8
5.8
3.9
5.7

12.5
10.9
11.9
9.0
10.2
10.1
11.4
11.0
8.3
9.1
9.0
9.7
9.2
12.1
10.6
11.8
8.4
10.4

Table 3: Empirical size (in %) for the entropy difference (ED) test under the null distribution GEVr , with µ = 0 and
σ = 1 based on 10,000 samples.
Sample Size

r

Nominal Size

50

100

1.0

5.0

10.0

1.0

5.0

10.0

ξ = −0.25

2
5
10

1.5
2.4
2.3

5.7
6.8
6.8

10.8
11.9
11.7

1.3
1.6
1.9

5.5
5.9
6.0

10.1
10.6
11.1

ξ=0

2
5
10

1.3
1.6
2.3

5.6
5.9
6.5

11.0
11.2
11.8

1.2
1.5
1.6

5.3
5.7
5.9

10.4
10.6
10.7

ξ = 0.25

2
5
10

1.3
1.6
2.0

5.7
5.8
6.6

10.7
11.5
11.9

1.3
1.3
1.4

5.4
5.6
5.5

10.5
10.2
10.4

14

TABLES

Table 4: Empirical rejection rate (in %) of the multiplier score test and the ED test in the first data generating
scheme from 1000 replicates.
Sample Size

ξ

−0.25

100

0
0.25
−0.25

200

0
0.25

Test

Value of a=b
0.4

0.6

0.8

1.0

1.2

1.4

1.6

1.8

2.0

Score
ED
Score
ED
Score
ED

99.9
100.0
100.0
100.0
100.0
100.0

84.8
99.0
87.0
98.8
87.7
97.5

20.4
46.5
21.6
40.0
20.3
37.7

5.4
4.6
7.4
5.2
6.2
4.8

21.0
48.7
24.2
40.6
25.8
34.8

41.2
89.5
48.9
87.2
54.2
78.1

62.3
99.2
67.8
98.5
74.2
96.1

79.0
100.0
79.6
100.0
82.9
99.5

83.0
99.8
89.4
99.7
89.5
99.7

Score
ED
Score
ED
Score
ED

100.0
100.0
100.0
100.0
100.0
100.0

98.6
100.0
99.4
99.9
99.3
100.0

40.7
78.4
44.6
75.0
44.5
71.0

5.2
6.2
6.1
5.5
6.3
5.2

29.8
70.0
34.9
64.6
37.0
57.2

64.7
99.2
75.0
98.1
73.4
95.9

86.4
100.0
92.4
99.8
91.8
100.0

95.9
100.0
97.3
100.0
97.0
100.0

97.5
100.0
98.6
100.0
98.9
100.0

Table 5: Empirical rejection rate (in %) of the multiplier score test and the ED tests in the second data generating
scheme from 1000 replicates.
Sample Size

100

ξ

−0.25
0
0.25

200

−0.25
0
0.25

Test

Mixing Rate p
0.00

0.10

0.25

0.50

0.75

0.90

1.00

Score
ED
Score
ED
Score
ED

99.7
100.0
100.0
100.0
99.9
100.0

99.5
100.0
99.7
100.0
99.7
100.0

95.5
100.0
97.8
100.0
96.6
99.9

69.4
97.7
72.4
96.0
70.8
93.6

24.1
51.8
22.7
47.6
24.7
43.4

7.8
10.9
6.2
10.3
5.8
9.8

5.8
6.2
6.8
5.6
5.3
5.2

Score
ED
Score
ED
Score
ED

99.9
100.0
100.0
100.0
100.0
100.0

100.0
100.0
100.0
100.0
100.0
100.0

99.7
100.0
100.0
100.0
100.0
100.0

95.6
100.0
96.5
100.0
97.2
99.7

43.4
83.6
44.4
79.5
46.9
72.5

11.4
20.0
11.2
20.0
9.2
17.9

5.1
5.8
5.4
5.5
5.5
4.2

6
5
4
3
2
1
0

6
5
4
3
2
1
0

6
5
4
3
2
1
0

PB Score

MB Score

Significance:

r

ED

Test

49.9
2.5
38.3
1.6
2.7
4.2
0.8

35.8
2.5
58.4
0.5
0.6
0.8
1.4

19.0
1.9
76.2
0.8
1.1
1.0
-

0.01

16.9
2.3
59.6
2.8
4.4
8.3
5.7

16.1
1.8
68.1
2.3
3.0
3.7
5.0

3.4
2.1
79.9
4.1
5.3
5.2
-

0.05

6.9
0.7
59.1
4.0
7.1
10.6
11.6

8.5
0.8
63.9
4.5
4.9
6.9
10.5

1.5
1.1
70.2
7.6
9.5
10.1
-

0.1

1.3
0.3
44.0
6.6
11.0
14.0
22.8

1.8
0.5
46.1
6.8
10.3
13.9
20.6

0.5
0.7
50.5
11.8
16.2
20.3
-

0.2

Unadjusted
0.3

0.2
0.1
31.2
7.5
10.0
19.5
31.5

0.6
0.2
31.3
7.1
11.7
18.3
30.8

0.0
0.2
35.2
15.3
19.9
29.4
-

0.0
0.0
18.4
6.0
10.6
20.0
45.0

0.2
0.1
19.8
8.0
12.7
18.5
40.7

0.0
0.1
22.3
15.1
22.8
39.7
-

0.4

71.7
1.3
26.6
0.3
0.1
0.0
0.0

53.5
1.9
42.6
1.4
0.5
0.1
0.0

86.6
1.3
12.0
0.1
0.0
0.0
-

0.01

0.05

40.3
2.0
53.3
2.8
0.6
0.9
0.1

33.0
1.3
57.8
4.2
1.7
1.3
0.7

69.8
1.5
25.1
3.4
0.1
0.1
-

24.7
0.8
59.3
7.4
3.4
2.7
1.7

23.4
0.7
58.9
9.0
3.0
2.3
2.7

58.5
1.1
31.7
6.2
1.9
0.6
-

0.1

12.6
0.5
49.9
15.7
8.6
4.7
8.0

13.9
0.8
51.1
15.5
7.3
4.6
6.8

43.0
0.9
33.4
12.4
5.4
4.9
-

0.2

ForwardStop
0.3

7.8
0.3
40.1
16.0
9.7
7.9
18.2

8.6
0.4
42.0
16.0
10.2
8.3
14.5

30.0
0.2
31.6
17.0
9.4
11.8
-

0.4

5.5
0.4
28.0
17.9
9.5
7.8
30.9

5.9
0.4
31.1
17.1
12.3
9.6
23.6

22.4
0.0
26.3
18.5
11.5
21.3
-

0.01

51.6
39.4
6.2
0.6
0.7
1.5
0.0

40.6
29.8
29.3
0.0
0.1
0.2
0.0

52.4
25.0
22.6
0.0
0.0
0.0
-

27.3
50.7
18.5
1.2
0.8
1.5
0.0

25.1
37.7
36.5
0.4
0.1
0.2
0.0

22.2
18.9
58.9
0.0
0.0
0.0
-

0.05

16.9
42.9
35.0
2.5
1.2
1.5
0.0

18.8
29.2
50.5
1.1
0.2
0.2
0.0

13.1
13.6
72.9
0.3
0.1
0.0
-

0.1

10.4
25.5
55.4
3.6
2.8
2.3
0.0

12.1
17.5
65.7
3.3
1.0
0.4
0.0

5.4
6.7
84.9
2.2
0.7
0.1
-

0.2

StrongStop
0.3

6.2
15.8
62.3
7.0
5.3
3.4
0.0

7.6
10.5
73.5
4.8
2.4
1.2
0.0

1.7
2.9
89.0
4.5
1.8
0.1
-

0.4

4.3
9.4
64.8
7.8
7.0
6.4
0.3

5.6
6.6
74.6
5.6
4.5
2.7
0.4

1.0
1.4
85.7
6.9
4.4
0.6
-

TABLES
15

Table 6: Percentage of choice of r using the ForwardStop and StrongStop rules at various significance levels or FDRs,
under ED, parametric bootstrap (PB) score, and multiplier bootstrap (MB) score tests, with n = 100 and ξ = 0.25
for the simulation setting described in Section 6. Correct choice is r = 4.

16

FIGURES

Shape: −0.25

Shape: 0.25

0.2
r: 1

0.1
0.0

r: 2

Density

0.2
0.1
0.0

0.2
r: 5

0.1
0.0
0

5

10

15

0

5

10

15

Score Statistic
Fig. 1: Comparisons of the empirical distribution based on 5000 replicates of the score test statistic and the χ2 (3)
distribution (red solid curve). The number of blocks used is n = 5000 with parameters µ = 0, σ = 1, and ξ ∈ (−0.25,
0.25).

FIGURES

1.00

0.25
0.00
1.00
0.75
0.50
0.25
0.00

n: 50

n: 100

n: 200

●
●●●●
●●
●●●●●
●●●●●●
●●●●●
●●●●●
●●●●
●●●●
●●●
●●●
●●●
●●●
●●●
●●●
●
●●
●●
●
●●
●●
●●
●●
●●
●
●●●
●●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●●
●
● ●●●
●●
●●●●
●●
●●●
●●
●●
●
●●
●●●
●●●●●●●
●●
●●
●●●●●
●●●
●●
●
●●
●●
●●
●
●●
●● ●●
●●
●●
●
●●
●●●●●●
●
●●
●● ●●●●
●
●●●●●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●●
●
●●●●●
●
●
●●
●●
●
●●
●●
●●
●
●●●●
●●
●
●
●●
●
●● ●●●
●●
●●●
●
●●●●
●
●
●●
●●
●
●●
●● ●●
●●
●●●●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
● ●
● ●● ●●
●
●●
● ●●
● ●●●●
●●
●●
●
● ●●
●●● ●●
●
●●
●●●●●
● ●
●●
●
●●●●
●● ●●
●
●
● ●●
●
● ●
● ●
●
●
●●
●●
●
●
●
●
●
●
●●●●●
●●●●●
●●●●●●
●●●●●●
●●●●
●●●●● ●●●●
●●●
●●●●
●
●●●
●●
●
●
●●●
●●●
●●●
●●●
●●● ●●●●●●●●●
●●●
●●●
●●
●●
●●
● ●
●●
●●
●●●
●●
●●
●●
●
●●●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●●●
●
●●●●
●
●●●
●●
●●
●●
●●
●●
●●
●●
●●
●●
●●
●●
●●
●●
●●●●●●
●
●● ●●● ●●●
●
●
●
●
●
●
●
●●
●●
●
● ●●
●
●●
●
●
●
●● ●●
●
●
●●
●
●●
●●
●●
●●
●
● ●●
●
● ●●●●●●
●●
●
●●
●●
●●
●●
● ●
●
●
●
●
● ●●●
●●●
●
●●
●●
●●●●●●●
●
●
●
●
●
●
●●
●
●
●
●●
●
●
●
●
●●●●●
●
● ●●
●
●
●● ●●
●
●
●
●
●●●●
●●●
●
●
●
●
●
●●●
●
●
●
●
●
●
●
●
●
●
●●
●
●●●
●
● ●●
● ●
●●
●●
●
●
●●
●●
●
●
●
●

Shape: 0.25

Observed FWER

0.50

n: 30

Shape: −0.25

0.75

17

0.05

0.15

0.25

0.05

0.15

0.25

0.05

Target FWER

0.15

0.25

0.05

0.15

0.25

Adjustment ● None StrongStop
Test ● ED ● MB Score ● PB Score
Fig. 2: Observed FWER for the ED, parametric bootstrap (PB) score, and multiplier bootstrap (MB) score tests
(using No Adjustment and StrongStop) versus expected FWER at various nominal levels. The 45 degree line indicates
agreement between the observed and expected rates under H0 .

n: 30

0.5

n: 50

n: 100

n: 200
Shape: −0.25

0.4
0.2
0.1
0.0
0.5

●
●●
●
●
●●●●
●●●
●●
●●●●●
●●●
●●●
●
●●● ●
●
●●●●
●●●●●
●●●●●●●●●●●
●
●●
●●●●●● ●●●
●●●
●●
●●●
●●●●
●●
●●
●●●
●●●
●●●● ●●●●●●
●●
●●●●●●●●●●●
●
●●●
●●●
●●●
●●●●
●●
●
●
●●
●
●
●
●
●
●
●
●
●●
●
●●
●
●
●
●●
●●
●●
●●
●
●●
●●●
●
●
●●
●
●
●
●●●
●●●
●●
●●
●●
●●
●●●●●●●●●●
●●●●●
●●●●●
●●●
●
●●●●●●
●●●●●
●●●●●●
●●
●●
●●
●●●●
●●
●●
●●
●●●
●●●
●●
●●●
●●●●
●●
●●●●●●●●
●●
●●●
●●●
●●
●●
●●●●●●●●●●
●●●
●●●●●●●●●●●
●●●●●●●●●●●
●●●
●●
●
●●
●●
●●●●
●●●●●●●
●●●●●●●
●●●●●●●●●●
●●●
●●●●

Shape: 0.25

Observed Level

0.3

0.4
0.3
0.2
0.1
0.0

●●
●●
●●●
●●●
●●●
●●●●
●●●●●
●●●●●
●●●●●
●
●●
●●●●●
●
●●●●●
●●●●●
●●
●●●●
●●●
●●●●
●●●●
●
●●●●●●●●●●
●
●
●●
●
●
●●●
●
●
●
●
●●
●●
●●
●
●
●
●
●
●
●
●●●●●●●●●●●
●●●
●●● ●●●
●●
●
●●●
●●●
●●●●
●●●●
●●
●●●
●●●
●●●
●●●●●
●●
●●●●●●●●●●●●●●
●●●●●●●●●●●●●●
●●●●● ●●●●●
●●●●
●●
●●●●●●
●●●●●●
●●●●●
●●●
●●●●
●●●●
●●●
●●●●
●●●
●●●
●●●●●
●●●●●●●●
●●●●●●●●●
●●●●
●●●●●●●
●●●●●●●●●
●●
●●●●●●●●●●●●
●●●
●
●●●●●
●●
●●●●●
●●●●●●●
●●●●●●●
●●●●
●●●●●●●●●
●●●●
●●●●●●●●●
●●●

0.05

0.15

0.25

0.05

0.15

0.25

0.05

Target Level

0.15

0.25

0.05

0.15

0.25

Control ForwardStop ● StrongStop
Test ● ED ● MB Score ● PB Score
Fig. 3: Observed FDR (from ForwardStop) and observed FWER (from StrongStop) versus expected FDR and FWER,
respectively, at various nominal levels. This is for the simulation setting described in Section 6, using the ED,
parametric bootstrap (PB) score, and multiplier bootstrap (MB) score tests. The 45 degree line indicates agreement
between the observed and expected rates.

18

FIGURES
●

0.8

●●●●●●●●●●●●●

●●●●

●●●●●●●●●●●●●●●●●●
●

●●

●●●●●

●●●

ForwardStop

●

0.6

●
●●●●

●●●●

0.4

●●

●●●●●●●●●●●●

●
●●●

●●

●●

●●●●●●●●

0.2

●●●

●●●●●●●●

●●●●●

●●

0.0
●●●●●●●●●●●●●●●●●●●●●●●●

●●●●●●●●

●●●●●●

●●●●

●●●●●

0.9

●●● ●●●●●●
●
●●●●●●●●

●●

●

0.6

●●●
●

●●●●
●●●●

●●●

0.3

●

●
●

0.0
1.00

●●

●

●

●

●

0.75

●

●●

●

●

●

●

●

●

●

●

●

●
●

0.00
0

10

●
●

20

●

●
●

●

●

●●

●

●

●

●●
●

●
●

●

●

●

●
●

40

●
●

●

50

60

●

●

●

70

●

80

●

●

●

●

●
●

●

●

● ● ●

●●●

●

●

30

●●●●●●●●

●

●
●

●

●
●

●●

●

●

●

0.25

●
●

●

●

●●

●

●●

●

●●

●

●

●

●

●

●

0.50

●

●

●

●

●

●
●

●

●

●

UnAdjusted

●

●●●●●●

●●●●●

StrongStop

P−Value

●●●●●●●●●●●●●

●
●●●●

●●●●●●●●●●●●●

● ●
● ●●

90

●
●

●

100

●

●
●●●● ●

110

●
● ● ●● ●
● ●●●
●

120

r
Test ED PB Score
●

Fig. 4: P-Values using ForwardStop, StrongStop, and no adjustment for the ED and PB Score tests applied to the
Lowestoft sea level data. The horizontal dashed line represents the 0.05 possible cutoff value.

FIGURES

19

Location

3.52

50 Year Return Level
5.6

3.48
●

3.44

●

● ● ● ● ● ●
● ● ● ● ● ● ● ● ● ● ● ●
● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●
● ●

5.2
4.8

● ●

3.40

● ● ● ●

4.4

● ● ● ● ●

● ● ● ● ●
● ● ● ● ● ● ●

● ● ● ● ● ● ● ● ● ● ● ●
● ● ● ● ● ● ●

3.36
Scale

0.300

6.0

0.275

Value

100 Year Return Level

6.5

0.250

●
●

0.225

● ●
● ● ● ● ●
● ●
● ● ●
● ● ● ● ● ● ●

0.200

5.5
● ● ● ● ● ● ● ● ● ●
● ● ● ● ● ● ●
● ●

5.0

● ● ● ●

4.5

● ● ● ● ●
● ● ● ● ●
● ● ● ● ● ● ●
● ● ● ● ● ● ● ● ● ● ● ● ●
● ● ● ● ● ●

0.175
Shape

200 Year Return Level

7

0.2
0.1

●

6

● ●
●

● ●
● ● ●

0.0

● ● ● ● ●
● ● ● ● ● ● ●

●

● ● ●
● ● ● ● ● ● ●
● ● ● ● ● ●
● ●

● ● ● ●

5

● ● ● ● ●

−0.1
0

5

10

15

20

25

30

35

40

0

5

● ● ● ● ●
● ● ● ● ● ● ●

10

15

20

● ● ● ● ● ● ● ● ● ● ● ● ●
● ● ● ● ● ●

25

30

35

40

r
Fig. 5: Location, scale, and shape parameter estimates, with 95% profile confidence intervals for r = 1, . . . , 40 for the
Lowestoft sea level data. Also included are the estimates and 95% profile likelihood confidence intervals for the 50
and 100 year return levels.

20

FIGURES

ForwardStop

2
●

●
●

1

●

●

●

●

●

●

●

StrongStop

P−Value

0
1.5
●

1.0

●
●

●

●

●

0.5
●

0.0
1.00

●
●
●

●

UnAdjusted

0.75
0.50

●
●

●

0.25

●

0.00
1

2

3

4

5

6

7

8

●

●

9

10

r
Test ED PB Score
●

Fig. 6: P-Values using ForwardStop, StrongStop, and no adjustment for the ED and PB Score tests applied to the
Atlantic City precipitation data. The horizontal dashed line represents the 0.05 possible cutoff value.

FIGURES

21

Location

50 Year Return Level
25.0

6.25
●

6.00

●

●

●

●

●

●

●

22.5

●

20.0

●

●

●

●

●
●

5.75

●

●

●

●

●

●

●

●

●

9

10

17.5

5.50
Scale

Value

100 Year Return Level

30

2.4

27

2.2
●

●

●

●

●

●

2.0

●

●

●

●

24

●

●
●

●
●

●

●

●

21
1.8

18
Shape

200 Year Return Level
35

0.4

30

0.3
●

●
●

●

●
●

0.2
1

●

●

2

3

4

5

●
●

●
●

6

●

7

●

8

●

●

9

10

●

●

25
1

2

3

4

5

6

7

8

r
Fig. 7: Location, scale, and shape parameter estimates, with 95% delta confidence intervals for r = 1 through r = 10
for the Atlantic City precipitation data. Also included are the estimates and 95% profile likelihood confidence intervals
for the 50, 100, and 200 year return levels.

