Block maxima and the GEV distribution

extreme value theory

statistics

ecology tutorial

Fit the generalised extreme value distribution to annual maxima in R from scratch, read the shape parameter, and see why the tail is hard to pin down.

Author

Tidy Ecology

Published

2026-08-17

Most of the statistics on this blog is about the middle of a distribution: the mean response, the typical slope, the average abundance. But a great deal of ecology is decided at the edges. A population is not killed by an average summer; it is killed by the hottest day in a decade. A shore community is not reshaped by a typical wave; it is reshaped by the largest storm in fifty years. Gaines and Denny (1993) made this point plainly: biostatistics trains ecologists to think about central tendency, yet many of the questions that matter concern the largest, smallest, highest and lowest values a system ever sees.

Extreme value theory gives those questions a proper footing. This post covers the block maxima approach: split a long record into blocks (usually years), keep the maximum of each block, and fit the one distribution that block maxima are guaranteed to converge to. We build the likelihood by hand so nothing is hidden, and we end on the honest catch that the whole cluster returns to: the parameter that governs the tail is the hardest one to estimate, and a few decades of maxima barely constrain it.

The one distribution block maxima can converge to

Suppose you take the maximum of a large number of independent observations. As the block size grows, the (suitably rescaled) maximum can only settle down to one family of distributions, whatever the underlying data look like. This is the extremal types theorem of Fisher and Tippett (1928), later given its general form by Gnedenko (1943), and it is the extreme value analogue of the central limit theorem. The limiting family is the generalised extreme value (GEV) distribution, with distribution function

\[G(z) = \exp\left\{-\left[1 + \xi\,\frac{z-\mu}{\sigma}\right]^{-1/\xi}\right\}, \qquad 1 + \xi\,\frac{z-\mu}{\sigma} > 0,\]

for location $\mu$, scale $\sigma > 0$ and shape $\xi$. The shape parameter is the whole story. It sorts the GEV into three tail behaviours:

$\xi = 0$ (the Gumbel limit): a light, exponential-like tail. Maxima of Gaussian, exponential and lognormal data head here.
$\xi > 0$ (the Frechet limit): a heavy, polynomial tail with no upper bound. Maxima of heavy-tailed data head here.
$\xi < 0$ (the Weibull limit): a bounded tail with a finite upper endpoint at $\mu - \sigma/\xi$.

Base R has no GEV functions, so we write the density, quantile and random generator ourselves. There is nothing to install.

library(ggplot2)

dgev <- function(x, mu, sigma, xi) {
  z <- (x - mu) / sigma
  if (abs(xi) < 1e-8) {
    exp(-z) * exp(-exp(-z)) / sigma
  } else {
    t <- 1 + xi * z
    d <- numeric(length(x))
    ok <- t > 0
    d[ok] <- (1 / sigma) * t[ok]^(-1 / xi - 1) * exp(-t[ok]^(-1 / xi))
    d
  }
}
qgev <- function(p, mu, sigma, xi) {
  if (abs(xi) < 1e-8) mu - sigma * log(-log(p))
  else mu - (sigma / xi) * (1 - (-log(p))^(-xi))
}
rgev <- function(n, mu, sigma, xi) qgev(runif(n), mu, sigma, xi)

# a Weibull-type GEV has a hard ceiling; a Frechet-type one does not
weibull_endpoint <- 31 - 1.6 / (-0.25)

A Weibull-type fit with $\mu = 31$, $\sigma = 1.6$ and $\xi = -0.25$ places a physical ceiling on the variable at 37.4: no matter how long you wait, the block maximum cannot exceed it. That is a strong biological claim, and whether it is warranted comes down to the sign of $\xi$, which is exactly the quantity the data struggle to resolve.

Fitting the GEV to annual maxima

Take a concrete ecological setting: the annual maximum daily temperature at a site, recorded for sixty years. High-temperature extremes drive thermal stress and mortality, so the quantity of interest is not the average summer but the hottest day each year. We simulate sixty annual maxima from a GEV with a mildly heavy tail, then recover the parameters as if we did not know them.

nll_gev <- function(par, x) {
  mu <- par[1]; sigma <- exp(par[2]); xi <- par[3]   # log-scale keeps sigma > 0
  z <- (x - mu) / sigma
  if (abs(xi) < 1e-6) {
    ll <- -log(sigma) - z - exp(-z)
  } else {
    t <- 1 + xi * z
    if (any(t <= 0)) return(1e10)
    ll <- -log(sigma) - (1 / xi + 1) * log(t) - t^(-1 / xi)
  }
  -sum(ll)
}
fit_gev <- function(x) {
  s0 <- sd(x) * sqrt(6) / pi                 # moment starts (Gumbel)
  mu0 <- mean(x) - 0.5772 * s0
  o <- optim(c(mu0, log(s0), 0.1), nll_gev, x = x,
             method = "BFGS", hessian = TRUE)
  se <- sqrt(diag(solve(o$hessian)))
  list(mu = o$par[1], sigma = exp(o$par[2]), xi = o$par[3],
       se_mu = se[1], se_xi = se[3])
}

set.seed(4246)
mu_true <- 31.0; sigma_true <- 1.6; xi_true <- 0.13
ymax <- rgev(60, mu_true, sigma_true, xi_true)   # sixty annual maxima
fit <- fit_gev(ymax)
xi_lo <- fit$xi - 1.96 * fit$se_xi
xi_hi <- fit$xi + 1.96 * fit$se_xi

The maximum likelihood fit gives a location of 30.82 (standard error 0.24) and a scale of 1.68, both close to the true 31.0 and 1.6. The location and scale are pinned down well: sixty maxima are plenty to say where the distribution sits and how spread out it is.

The shape is another matter. The estimate is 0.060 with a standard error of 0.081, so a Wald interval runs from -0.10 to 0.22. That interval covers the true value of 0.13, but it also comfortably covers zero and a good stretch of negative values. From sixty years of maxima we cannot confidently say whether the tail is heavy, light or bounded, even though we simulated it as heavy. This is not a failing of the code; it is intrinsic. The shape parameter is estimated from the handful of largest maxima, and a few decades simply do not contain many of those.

Left panel shows three GEV density curves for Gumbel, Frechet and Weibull shapes, with the Weibull curve cut off at a vertical dotted bound. Right panel shows a histogram of sixty annual maxima with a smooth generalised extreme value density curve fitted on top. — Figure 1: The generalised extreme value distribution. Left: the three tail types set by the shape parameter, with the Weibull type ending at a finite upper bound (dotted line). Right: the fitted GEV density (dark curve) over the sixty simulated annual maxima; location and scale are well estimated, the tail shape is not.

Return levels: reading the fit as a design number

The reason to fit a GEV is rarely the density itself. It is the return level: the value expected to be exceeded once, on average, every $m$ blocks. For annual maxima that is the $m$-year event. The $m$-year return level is the $(1 - 1/m)$ quantile of the GEV,

\[z_m = \mu - \frac{\sigma}{\xi}\left[1 - \{-\log(1 - 1/m)\}^{-\xi}\right],\]

which is just our quantile function evaluated at $p = 1 - 1/m$.

return_level <- function(fit, m) qgev(1 - 1 / m, fit$mu, fit$sigma, fit$xi)
rl50  <- return_level(fit, 50)
rl100 <- return_level(fit, 100)
rl50_true  <- qgev(1 - 1 / 50,  mu_true, sigma_true, xi_true)
rl100_true <- qgev(1 - 1 / 100, mu_true, sigma_true, xi_true)

The fitted 50-year return level is 38.2 and the 100-year level is 39.7. The true values, which we happen to know, are 39.1 and 41.1. The fit underestimates both, and the reason is the shape: our $\hat\xi$ came out below the truth, and the return level for a long return period is dominated by $\xi$. A small change in the shape, well within its standard error, moves the 100-year temperature by more than a degree. The return level inherits all of the shape parameter’s uncertainty and then some, because we are extrapolating past the largest value we ever observed. The uncertainty in that extrapolation is the subject of a later tutorial.

The block size is a bias against variance trade

The GEV is a limit as the block size grows. For a finite block it is an approximation, and the quality of the approximation depends on the underlying data. Gaussian data are the notorious case: the maximum of a Gaussian sample converges to the Gumbel form ($\xi = 0$), but it does so at a logarithmic rate, which is painfully slow. Fit a GEV to maxima taken over short blocks and $\hat\xi$ is pulled well below zero, wrongly suggesting a bounded tail.

We can see the systematic bias directly by taking maxima over blocks of increasing length and looking at the average $\hat\xi$, rather than relying on any single noisy dataset.

set.seed(202)
block_sizes <- c(10, 30, 90, 365, 2000)
bias_tab <- data.frame(block = block_sizes, xi = NA_real_, sd = NA_real_)
for (i in seq_along(block_sizes)) {
  xis <- replicate(30, {
    M <- replicate(400, max(rnorm(block_sizes[i])))
    fit_gev(M)$xi
  })
  bias_tab$xi[i] <- mean(xis)
  bias_tab$sd[i] <- sd(xis)
}
xi_small <- bias_tab$xi[bias_tab$block == 10]
xi_big   <- bias_tab$xi[bias_tab$block == 2000]

Blocks of ten give an average shape of -0.163, blocks of two thousand give -0.067, and the true limit is exactly zero. Even blocks of two thousand have not arrived. This is the bias side of the trade. The variance side is the mirror image: for a fixed length of record, longer blocks mean fewer maxima to fit, so the estimates get noisier. There is no block size that is right in general. Annual blocks are the convention because a year is a natural unit and it usually keeps the approximation acceptable, but the choice is a modelling decision, not a fact about the data.

Left panel plots return level against return period on a logarithmic axis, a rising curve with observed points scattered around it. Right panel plots estimated shape parameter against block size, rising from about minus 0.16 towards zero. — Figure 2: Left: the return level plot for the fitted GEV, with empirical maxima placed at their plotting-position return periods. Right: the shape parameter estimated from Gaussian maxima creeps towards its true value of zero only slowly as the block size grows.

What to take away

Block maxima give a principled way to model the largest value a system produces in each period, and the GEV is the one distribution they are entitled to follow. Fitting it by hand is a short piece of code. The honest parts are the ones worth remembering. The location and scale are easy; the shape parameter, which decides whether the tail is heavy, light or bounded, and which drives every long return level, is estimated from only the largest few maxima and carries wide uncertainty. The block size trades bias for variance with no objective optimum. And a return level is an extrapolation beyond the observed data, so it is only as trustworthy as the shape parameter that steers it. The next tutorial keeps every value above a high threshold instead of one per block, which uses the data more efficiently and shifts the difficulty from the block size to the threshold.

References

Fisher, R.A. and Tippett, L.H.C. 1928. Proceedings of the Cambridge Philosophical Society 24(2):180-190 (10.1017/S0305004100015681).

Gnedenko, B. 1943. Annals of Mathematics 44(3):423-453 (10.2307/1968974).

Gaines, S.D. and Denny, M.W. 1993. Ecology 74(6):1677-1692 (10.2307/1939926).

Katz, R.W., Brush, G.S. and Parlange, M.B. 2005. Ecology 86(5):1124-1134 (10.1890/04-0606).

Coles, S. 2001. An Introduction to Statistical Modeling of Extreme Values. Springer (ISBN 978-1-85233-459-8).

--- title: "Block maxima and the GEV distribution" description: "Fit the generalised extreme value distribution to annual maxima in R from scratch, read the shape parameter, and see why the tail is hard to pin down." date: "2026-08-17 09:00" categories: [R, extreme value theory, statistics, ecology tutorial] image: thumbnail.png image-alt: "Three generalised extreme value densities and a fitted curve over annual temperature maxima" --- Most of the statistics on this blog is about the middle of a distribution: the mean response, the typical slope, the average abundance. But a great deal of ecology is decided at the edges. A population is not killed by an average summer; it is killed by the hottest day in a decade. A shore community is not reshaped by a typical wave; it is reshaped by the largest storm in fifty years. Gaines and Denny (1993) made this point plainly: biostatistics trains ecologists to think about central tendency, yet many of the questions that matter concern the largest, smallest, highest and lowest values a system ever sees. Extreme value theory gives those questions a proper footing. This post covers the block maxima approach: split a long record into blocks (usually years), keep the maximum of each block, and fit the one distribution that block maxima are guaranteed to converge to. We build the likelihood by hand so nothing is hidden, and we end on the honest catch that the whole cluster returns to: the parameter that governs the tail is the hardest one to estimate, and a few decades of maxima barely constrain it. ## The one distribution block maxima can converge to Suppose you take the maximum of a large number of independent observations. As the block size grows, the (suitably rescaled) maximum can only settle down to one family of distributions, whatever the underlying data look like. This is the extremal types theorem of Fisher and Tippett (1928), later given its general form by Gnedenko (1943), and it is the extreme value analogue of the central limit theorem. The limiting family is the generalised extreme value (GEV) distribution, with distribution function $$G(z) = \exp\left\{-\left[1 + \xi\,\frac{z-\mu}{\sigma}\right]^{-1/\xi}\right\}, \qquad 1 + \xi\,\frac{z-\mu}{\sigma} > 0,$$ for location $\mu$, scale $\sigma > 0$ and shape $\xi$. The shape parameter is the whole story. It sorts the GEV into three tail behaviours: - $\xi = 0$ (the Gumbel limit): a light, exponential-like tail. Maxima of Gaussian, exponential and lognormal data head here. - $\xi > 0$ (the Frechet limit): a heavy, polynomial tail with no upper bound. Maxima of heavy-tailed data head here. - $\xi < 0$ (the Weibull limit): a bounded tail with a finite upper endpoint at $\mu - \sigma/\xi$. Base R has no GEV functions, so we write the density, quantile and random generator ourselves. There is nothing to install. ```{r} #| label: gev-helpers #| message: false #| warning: false library(ggplot2) dgev <- function(x, mu, sigma, xi) { z <- (x - mu) / sigma if (abs(xi) < 1e-8) { exp(-z) * exp(-exp(-z)) / sigma } else { t <- 1 + xi * z d <- numeric(length(x)) ok <- t > 0 d[ok] <- (1 / sigma) * t[ok]^(-1 / xi - 1) * exp(-t[ok]^(-1 / xi)) d } } qgev <- function(p, mu, sigma, xi) { if (abs(xi) < 1e-8) mu - sigma * log(-log(p)) else mu - (sigma / xi) * (1 - (-log(p))^(-xi)) } rgev <- function(n, mu, sigma, xi) qgev(runif(n), mu, sigma, xi) # a Weibull-type GEV has a hard ceiling; a Frechet-type one does not weibull_endpoint <- 31 - 1.6 / (-0.25) ``` A Weibull-type fit with $\mu = 31$, $\sigma = 1.6$ and $\xi = -0.25$ places a physical ceiling on the variable at `r sprintf("%.1f", weibull_endpoint)`: no matter how long you wait, the block maximum cannot exceed it. That is a strong biological claim, and whether it is warranted comes down to the sign of $\xi$, which is exactly the quantity the data struggle to resolve. ```{r} #| label: setup-theme #| echo: false te_forest <- "#275139"; te_rust <- "#9c5b34"; te_gold <- "#b8860b" theme_te <- function() { theme_minimal(base_size = 12) + theme(panel.grid.minor = element_blank(), panel.grid.major = element_line(colour = "grey90", linewidth = 0.3), plot.title = element_text(face = "bold", size = 11.5), axis.title = element_text(size = 10.5), legend.position = "top", legend.title = element_blank(), plot.margin = margin(6, 10, 6, 6)) } ``` ## Fitting the GEV to annual maxima Take a concrete ecological setting: the annual maximum daily temperature at a site, recorded for sixty years. High-temperature extremes drive thermal stress and mortality, so the quantity of interest is not the average summer but the hottest day each year. We simulate sixty annual maxima from a GEV with a mildly heavy tail, then recover the parameters as if we did not know them. ```{r} #| label: sim-fit nll_gev <- function(par, x) { mu <- par[1]; sigma <- exp(par[2]); xi <- par[3] # log-scale keeps sigma > 0 z <- (x - mu) / sigma if (abs(xi) < 1e-6) { ll <- -log(sigma) - z - exp(-z) } else { t <- 1 + xi * z if (any(t <= 0)) return(1e10) ll <- -log(sigma) - (1 / xi + 1) * log(t) - t^(-1 / xi) } -sum(ll) } fit_gev <- function(x) { s0 <- sd(x) * sqrt(6) / pi # moment starts (Gumbel) mu0 <- mean(x) - 0.5772 * s0 o <- optim(c(mu0, log(s0), 0.1), nll_gev, x = x, method = "BFGS", hessian = TRUE) se <- sqrt(diag(solve(o$hessian))) list(mu = o$par[1], sigma = exp(o$par[2]), xi = o$par[3], se_mu = se[1], se_xi = se[3]) } set.seed(4246) mu_true <- 31.0; sigma_true <- 1.6; xi_true <- 0.13 ymax <- rgev(60, mu_true, sigma_true, xi_true) # sixty annual maxima fit <- fit_gev(ymax) xi_lo <- fit$xi - 1.96 * fit$se_xi xi_hi <- fit$xi + 1.96 * fit$se_xi ``` The maximum likelihood fit gives a location of `r sprintf("%.2f", fit$mu)` (standard error `r sprintf("%.2f", fit$se_mu)`) and a scale of `r sprintf("%.2f", fit$sigma)`, both close to the true `r sprintf("%.1f", mu_true)` and `r sprintf("%.1f", sigma_true)`. The location and scale are pinned down well: sixty maxima are plenty to say where the distribution sits and how spread out it is. The shape is another matter. The estimate is `r sprintf("%.3f", fit$xi)` with a standard error of `r sprintf("%.3f", fit$se_xi)`, so a Wald interval runs from `r sprintf("%.2f", xi_lo)` to `r sprintf("%.2f", xi_hi)`. That interval covers the true value of `r sprintf("%.2f", xi_true)`, but it also comfortably covers zero and a good stretch of negative values. From sixty years of maxima we cannot confidently say whether the tail is heavy, light or bounded, even though we simulated it as heavy. This is not a failing of the code; it is intrinsic. The shape parameter is estimated from the handful of largest maxima, and a few decades simply do not contain many of those. ```{r} #| label: fig-gev #| echo: false #| fig-width: 8.6 #| fig-height: 3.7 #| fig-cap: "The generalised extreme value distribution. Left: the three tail types set by the shape parameter, with the Weibull type ending at a finite upper bound (dotted line). Right: the fitted GEV density (dark curve) over the sixty simulated annual maxima; location and scale are well estimated, the tail shape is not." #| fig-alt: "Left panel shows three GEV density curves for Gumbel, Frechet and Weibull shapes, with the Weibull curve cut off at a vertical dotted bound. Right panel shows a histogram of sixty annual maxima with a smooth generalised extreme value density curve fitted on top." xx <- seq(26, 42, length.out = 400) fam <- rbind( data.frame(x = xx, d = dgev(xx, 31, 1.6, 0), type = "Gumbel (xi = 0)"), data.frame(x = xx, d = dgev(xx, 31, 1.6, 0.35), type = "Frechet (xi > 0)"), data.frame(x = xx, d = dgev(xx, 31, 1.6, -0.25), type = "Weibull (xi < 0)") ) pA <- ggplot(fam, aes(x, d, colour = type)) + geom_line(linewidth = 0.9) + geom_vline(xintercept = weibull_endpoint, linetype = "dotted", colour = te_rust) + scale_colour_manual(values = c("Gumbel (xi = 0)" = "grey35", "Frechet (xi > 0)" = te_gold, "Weibull (xi < 0)" = te_forest)) + labs(title = "Three tail types, one shape parameter", x = "block maximum", y = "density") + theme_te() grid_x <- seq(min(ymax) - 1, max(ymax) + 3, length.out = 400) fit_df <- data.frame(x = grid_x, d = dgev(grid_x, fit$mu, fit$sigma, fit$xi)) pB <- ggplot(data.frame(y = ymax), aes(y)) + geom_histogram(aes(y = after_stat(density)), bins = 16, fill = "grey85", colour = "white") + geom_line(data = fit_df, aes(x, d), colour = te_forest, linewidth = 1) + labs(title = "Fitted GEV over sixty annual maxima", x = "annual maximum temperature", y = "density") + theme_te() library(grid) grid.newpage() pushViewport(viewport(layout = grid.layout(1, 2))) print(pA, vp = viewport(layout.pos.row = 1, layout.pos.col = 1)) print(pB, vp = viewport(layout.pos.row = 1, layout.pos.col = 2)) ``` ## Return levels: reading the fit as a design number The reason to fit a GEV is rarely the density itself. It is the return level: the value expected to be exceeded once, on average, every $m$ blocks. For annual maxima that is the $m$-year event. The $m$-year return level is the $(1 - 1/m)$ quantile of the GEV, $$z_m = \mu - \frac{\sigma}{\xi}\left[1 - \{-\log(1 - 1/m)\}^{-\xi}\right],$$ which is just our quantile function evaluated at $p = 1 - 1/m$. ```{r} #| label: return-levels return_level <- function(fit, m) qgev(1 - 1 / m, fit$mu, fit$sigma, fit$xi) rl50 <- return_level(fit, 50) rl100 <- return_level(fit, 100) rl50_true <- qgev(1 - 1 / 50, mu_true, sigma_true, xi_true) rl100_true <- qgev(1 - 1 / 100, mu_true, sigma_true, xi_true) ``` The fitted 50-year return level is `r sprintf("%.1f", rl50)` and the 100-year level is `r sprintf("%.1f", rl100)`. The true values, which we happen to know, are `r sprintf("%.1f", rl50_true)` and `r sprintf("%.1f", rl100_true)`. The fit underestimates both, and the reason is the shape: our $\hat\xi$ came out below the truth, and the return level for a long return period is dominated by $\xi$. A small change in the shape, well within its standard error, moves the 100-year temperature by more than a degree. The return level inherits all of the shape parameter's uncertainty and then some, because we are extrapolating past the largest value we ever observed. The uncertainty in that extrapolation is the subject of a [later tutorial](../return-levels-and-uncertainty/). ## The block size is a bias against variance trade The GEV is a limit as the block size grows. For a finite block it is an approximation, and the quality of the approximation depends on the underlying data. Gaussian data are the notorious case: the maximum of a Gaussian sample converges to the Gumbel form ($\xi = 0$), but it does so at a logarithmic rate, which is painfully slow. Fit a GEV to maxima taken over short blocks and $\hat\xi$ is pulled well below zero, wrongly suggesting a bounded tail. We can see the systematic bias directly by taking maxima over blocks of increasing length and looking at the average $\hat\xi$, rather than relying on any single noisy dataset. ```{r} #| label: block-bias set.seed(202) block_sizes <- c(10, 30, 90, 365, 2000) bias_tab <- data.frame(block = block_sizes, xi = NA_real_, sd = NA_real_) for (i in seq_along(block_sizes)) { xis <- replicate(30, { M <- replicate(400, max(rnorm(block_sizes[i]))) fit_gev(M)$xi }) bias_tab$xi[i] <- mean(xis) bias_tab$sd[i] <- sd(xis) } xi_small <- bias_tab$xi[bias_tab$block == 10] xi_big <- bias_tab$xi[bias_tab$block == 2000] ``` Blocks of ten give an average shape of `r sprintf("%+.3f", xi_small)`, blocks of two thousand give `r sprintf("%+.3f", xi_big)`, and the true limit is exactly zero. Even blocks of two thousand have not arrived. This is the bias side of the trade. The variance side is the mirror image: for a fixed length of record, longer blocks mean fewer maxima to fit, so the estimates get noisier. There is no block size that is right in general. Annual blocks are the convention because a year is a natural unit and it usually keeps the approximation acceptable, but the choice is a modelling decision, not a fact about the data. ```{r} #| label: fig-returnbias #| echo: false #| fig-width: 8.6 #| fig-height: 3.7 #| fig-cap: "Left: the return level plot for the fitted GEV, with empirical maxima placed at their plotting-position return periods. Right: the shape parameter estimated from Gaussian maxima creeps towards its true value of zero only slowly as the block size grows." #| fig-alt: "Left panel plots return level against return period on a logarithmic axis, a rising curve with observed points scattered around it. Right panel plots estimated shape parameter against block size, rising from about minus 0.16 towards zero." mgrid <- exp(seq(log(2), log(200), length.out = 120)) rl_df <- data.frame(m = mgrid, z = sapply(mgrid, function(m) return_level(fit, m))) n <- length(ymax) emp <- data.frame(z = sort(ymax), p = (seq_len(n) - 0.44) / (n + 0.12)) emp$m <- 1 / (1 - emp$p) pC <- ggplot(rl_df, aes(m, z)) + geom_line(colour = te_forest, linewidth = 1) + geom_point(data = emp, aes(m, z), colour = te_rust, size = 1.4, alpha = 0.8) + scale_x_log10() + labs(title = "Return level plot", x = "return period (years, log scale)", y = "return level") + theme_te() pD <- ggplot(bias_tab, aes(block, xi)) + geom_hline(yintercept = 0, linetype = "dashed", colour = "grey55") + geom_errorbar(aes(ymin = xi - sd, ymax = xi + sd), width = 0.06, colour = "grey70") + geom_line(colour = te_gold, linewidth = 0.9) + geom_point(colour = te_gold, size = 2) + scale_x_log10() + labs(title = "Shape bias shrinks slowly with block size", x = "block size (log scale)", y = "estimated shape") + theme_te() grid.newpage() pushViewport(viewport(layout = grid.layout(1, 2))) print(pC, vp = viewport(layout.pos.row = 1, layout.pos.col = 1)) print(pD, vp = viewport(layout.pos.row = 1, layout.pos.col = 2)) ``` ## What to take away Block maxima give a principled way to model the largest value a system produces in each period, and the GEV is the one distribution they are entitled to follow. Fitting it by hand is a short piece of code. The honest parts are the ones worth remembering. The location and scale are easy; the shape parameter, which decides whether the tail is heavy, light or bounded, and which drives every long return level, is estimated from only the largest few maxima and carries wide uncertainty. The block size trades bias for variance with no objective optimum. And a return level is an extrapolation beyond the observed data, so it is only as trustworthy as the shape parameter that steers it. The [next tutorial](../peaks-over-threshold-and-the-gpd/) keeps every value above a high threshold instead of one per block, which uses the data more efficiently and shifts the difficulty from the block size to the threshold. ## References Fisher, R.A. and Tippett, L.H.C. 1928. Proceedings of the Cambridge Philosophical Society 24(2):180-190 (10.1017/S0305004100015681). Gnedenko, B. 1943. Annals of Mathematics 44(3):423-453 (10.2307/1968974). Gaines, S.D. and Denny, M.W. 1993. Ecology 74(6):1677-1692 (10.2307/1939926). Katz, R.W., Brush, G.S. and Parlange, M.B. 2005. Ecology 86(5):1124-1134 (10.1890/04-0606). Coles, S. 2001. An Introduction to Statistical Modeling of Extreme Values. Springer (ISBN 978-1-85233-459-8). ## Related tutorials - [Peaks over threshold and the GPD](../peaks-over-threshold-and-the-gpd/) - [Return levels and uncertainty](../return-levels-and-uncertainty/) - [Checking an extreme value model](../checking-an-extreme-value-model/) - [Quantile regression in ecology](../quantile-regression-in-ecology/)