Using Gumbel distribution to approximate distribution of sample maximum --- formulae for the parameters?

Question

Suppose you have an observable sample $X_1,...,X_n \sim \text{IID } F_X$ which has a right-tail that decreases sufficiently rapidly to apply the extreme-value theorem (e.g., a normal distribution) to give the approximation:

$$\max \{ X_1,...,X_n \} \overset{\text{approx}}{\sim} \text{Gumbel}(\psi_n, \theta_n),$$

where $\psi_n$ and $\theta_n$ are approriate parameters to approximate the true distribution of the maximum. As far as I am aware, there are a number of methods used to determine the appropriate parameters for the approximation (see e.g., here). However, I am not sure which formulae are the most accurate and/or the most simple.

---

My question: What are some of the methods used to obtain the parameters of the Gumbel approximation to the sample maximum? (One method per answer please.) Does the method result in an explicit formula for the parameters, or does it require numerical evaluation of the parameters? How well do these approximations perform?

Answer 1

One method of approximation is to equate the true density to the approximation at two distinct quantiles, and then use these equations to obtain the two parameters of the Gumbel distribution. The advantage of this method is that it results in an explicit formula for the parameters, which means that it is simple to implement. In this answer I will derive the relevant formulae for the method and give an example to show the accuracy of the method for a sample of normal random variables.

---

Parameter formulae via quantile correspondence: To derive the formulae for this method, let's start by looking at what happens when we equate the two distributions at some point $x \in \mathbb{R}$ which gives the quantile $0<p<1$. This gives the equations:

$$p = F_X(x)^n = \exp \Big( -\exp \Big( - \frac{x-\psi}{\theta} \Big) \Big).$$

Let $Q(p) \equiv |\ln|\ln (p)||$ denote the quantile function of the standard Gumbel distribution, and let $Q_X = F_X^{-1}$ be the quantile function for the distribution of the observable values. Then we can write these equations as:

$$Q_X (p^{1/n}) = x \quad \quad \quad Q(p) = \frac{x-\psi}{\theta}.$$

Substituting the values of $x$ are re-arranging gives the parameter equation:

$$\psi = Q_X (p^{1/n}) - \theta \cdot Q(p).$$

This gives us a single equation for the parameters of the Gumbel distribution when we equate the true distribution to the Gumbel distribution at the quantile $p$. Now, since there are two parameters, we need to set the distributions equal at two distinct quantiles $p_0$ and $p_1$. Solving the resulting simultaneous equations gives the parameter estimators:

$$\hat{\psi}_n = \frac{Q_X(p_1^{1/n}) Q(p_0) - Q_X(p_0^{1/n}) Q(p_1)}{Q(p_0)-Q(p_1)} \quad \quad \quad \hat{\theta}_n = \frac{Q_X(p_0^{1/n}) - Q_X(p_1^{1/n})}{Q(p_0)-Q(p_1)}.$$

This gives explicit formulae for the parameters for any two distinct quantiles. Choice of different quantiles leads to different approximations, so the accuracy of the method depends on this choice. Generally we will want to use quantiles that are spaced apart by a reasonable amount so that the approximating distribution is close to the true distribution over a wide range of quantiles. A simple choice is to take $p_0 = \tfrac{1}{3}$ and $p_1 = \tfrac{2}{3}$ which ensures that the approximating distribution is equal to the true distribution at four equidistant quantiles (those two quantiles, plus the quantiles $p=0$ and $p=1$).

---

Implementation: We can implement this distributional approximation in R as follows. We create a function dmaxapprox which gives the density of the approximating distribution. This function takes inputs for the sample size n the quantile function qx and the quantiles p0 and p1 and uses these to form the Gumbel approximation to the true distribution of the maximum sample value. The function gives the density at a vector of points x and also allows the user to set the logical variable log to stipulate whether to return the density or log-density.

dmaxapprox <- function(x, n, qx, p0, p1) {

  #Set the parameters of the approximating distribution
  Q0    <- -log(-log(p0));
  Q1    <- -log(-log(p1));
  T0    <- qx(p0^(1/n));
  T1    <- qx(p1^(1/n));
  THETA <- (T0-T1)/(Q0-Q1);
  PSI   <- (T1*Q0-T0*Q1)/(Q0-Q1);

  #Output the approximating density
  extraDistr::dgumbel(x, mu = PSI, sigma = THETA, log = FALSE); }

We will test this function by comparing it to the kernel density estimator from $M = 10^6$ simulations of the maximum of a sample of $n=1000$ standard normal random variables, which is a large enough number of simulations that we will take this as the true density. We will generate the approximating density using the quantiles $p_0 = \tfrac{1}{3}$ and $p_1 = \tfrac{2}{3}$ as recommended above. We plot both densities overlayed on the same plot; the true density is shown by the unbroken line and the approximating density is shown by the dotted line.

#Set parameters
M <- 10^6;
n <- 1000;

#Simulate maximum values from standard normal distribution
set.seed(1);
RAND <- matrix(rnorm(n*M), nrow = M, ncol = n);
MAX  <- rep(NA, M);
for (i in 1:M) { MAX[i] <- max(RAND[i,]); }
DENS_SIM <- density(MAX);

#Generate approximating distribution
qx <- qnorm;
p0 <- 1/3;
p1 <- 2/3;
xx <- DENS_SIM$x;
yy <- dmaxapprox(xx, n, qx, p0, p1);
DENS_APPROX <- list(x = DENS_SIM$x, y = yy, n = NULL, bw = NULL,
                    data.name = NULL, has.na = FALSE);
class(DENS_APPROX) <- 'density';

#Plot the densities
plot(DENS_SIM, ylim = c(0, 1.5), lty = 1,
     main = 'True Density vs Approximating Density');
lines(DENS_APPROX, lty = 2, add = TRUE);

As can be seen from the figure, this particular method gives a reasonable approximation to the true density, by setting the parameters so that an equivalence in the quantiles holds at two interior points of the distribution. The method has the advantage that it leads to an explicit formula that does not involve any approximation --- the quantiles are exactly equivalent at the chosen points.