The median absolute deviation value of the Gumbel distribution

Andrey Akinshin · 2020-10-06

The Gumbel distribution is not only a useful model in the extreme value theory, but it’s also a nice example of a slightly right-skewed distribution (skewness $\approx 1.14$). Here is its density plot:

In some of my statistical experiments, I like to use the Gumbel distribution as a sample generator for hypothesis checking or unit tests. I also prefer the median absolute deviation (MAD) over the standard deviation as a measure of dispersion because it’s more robust in the case of non-parametric distributions. Numerical hypothesis verification often requires the exact value of the median absolute deviation of the original distribution. I didn’t find this value in the reference tables, so I decided to do another exercise and derive it myself. In this post, you will find a short derivation and the result (spoiler: the exact value is 0.767049251325708 * β). The general approach of the MAD derivation is common for most distributions, so it can be easily reused.

General approach

The median absolute deviation of a sample has a simple definition:

$$ \mathcal{MAD}_0 = \textrm{Median}(|x_i - \textrm{Median}(x_i)|). $$

The exact formula for the distribution looks similar:

$$ \mathcal{MAD}_0 = \textrm{Median}(|X - \textrm{Median}(X)|). $$

It’s convenient to use a scaled version of MAD in many practical cases to make it a consistent estimator for the standard deviation estimation. It’s often denoted as $\mathcal{MAD}$ without any clarification. To avoid misinterpretation and highlight that we are working with a non-scaled version of MAD, we denote it as $\mathcal{MAD}_0$ instead of just $\mathcal{MAD}$.

We also denote the median distribution value as $M = \textrm{Median}(X)$. Thus, the $\mathcal{MAD}_0$ definition can be rewritten as

$$ \mathcal{MAD}_0 = \textrm{Median}(|X - M|). $$

Since $\mathcal{MAD}_0$ is the median of $|X - M|$, we can conclude that the range $[M - \mathcal{MAD}_0; M + \mathcal{MAD}_0]$ contains $50\%$ of the distribution:

This can be expressed using the distribution CDF (let’s call it $F$):

$$ \begin{equation} \label{eq:main}\tag{1} F(M + \mathcal{MAD}_0) - F(M - \mathcal{MAD}_0) = 0.5. \end{equation} $$

Equation ($\ref{eq:main}$) is an important property of the median absolute deviation, which we will use to calculate it’s exact value.

Gumbel distribution

The Gumbel distribution is parametrized by $\mu$ (location) and $\beta$ (scale). For these parameters, the median value and the CDF are well-known (see [HoSM], 1.3.6.6.16):

$$ M = \mu - \beta \cdot \ln(\ln(2)), \quad F(x) = e^{-e^{-(x-\mu)/\beta}}. $$

The $\mathcal{MAD}_0$ value is directly proportional to the scale parameter $\beta$, so let’s introduce an auxiliary variable $p$ for the “descaled version” of $\mathcal{MAD}_0$:

$$ p = \mathcal{MAD}_0 / \beta. $$

Next, let’s express $F(M + \mathcal{MAD}_0)$ via $p$:

$$ \begin{split} F(M + \mathcal{MAD}_0) & = F(\mu - \beta \cdot \ln(\ln(2)) + p \beta) = \\ & = e^{-e^{-(\mu - \beta \cdot \ln(\ln(2)) + p \beta - \mu)/\beta}} = \\ & = e^{-e^{\ln(\ln(2)) - p}} = e^{-\ln(2) e^{-p}} = 0.5^{e^{-p}}, \end{split} $$

In the same way, we can show that

$$ F(M - \mathcal{MAD}_0) = 0.5^{e^{p}}. $$

Now, equation ($\ref{eq:main}$) can be transformed to

$$ 0.5^{e^{-p}} - 0.5^{e^{p}} = 0.5. $$

This equation can be solved numerically. Here is a short R script which calculates the result:

library(rootSolve)
options(digits = 15)
f <- function(p) 0.5^exp(-p) - 0.5^exp(p) - 0.5
uniroot.all(f, c(-10, 10), tol = 1e-15)

The numerical solution is $p = 0.767049251325708 $. Since $p = \mathcal{MAD}_0 / \beta$, the exact solution looks as follows:

$$ \mathcal{MAD}_0 = 0.767049251325708 \beta. $$

Now we can use it in experiments (copy-pastable version: 0.767049251325708 * β).

Numerical simulation

Let’s double-check that our calculations are correct. Below you can see an R script that generates 1000 random samples from the Gumbel distribution with $\mu = 0,\; \beta = 1$ (1000 elements in each sample). Next, it calculates the MAD estimation for each sample using the Harrell-Davis quantile estimator for the median estimations (see A new distribution-free quantile estimator
By Frank E Harrell, C E Davis · 1982harrell1982). After that, it calculates the median of all MAD estimations and prints the result.

library(evd)
library(Hmisc)

# Setting seed for rgumbel to achieve deterministic result
set.seed(38)
# Harrell-Davis-powered median
hdmedian <- function(x) as.numeric(hdquantile(x, 0.5))
# MAD value which is based on the Harrell-Davis-powered median
hdmad <- function(x) hdmedian(abs(x - hdmedian(x)))
# 1000 trails of MAD estimating for a sample from the Gubmel distribution
mads <- sapply(1:1000, function(i) hdmad(rgumbel(1000)))
# Print median of all trials
hdmedian(mads)

This script’s output is 0.7670284 which is close enough to the exact value 0.767049251325708.

Conclusion

The Gumbel distribution can be used as a good model of a slightly right-skewed distribution. If we are playing with an algorithm that involves MAD estimations (e.g., MAD-based outlier detector), we can check the precision of our calculations using the exact MAD value which is $\mathcal{MAD}_0 = 0.767049251325708 \beta$.

The general approach can be reused for any other distributions with known median and CDF. You can get the exact $\mathcal{MAD}_0$ value by solving the equation $F(M + \mathcal{MAD}_0) - F(M - \mathcal{MAD}_0) = 0.5$ which is correct for the distribution of any form.

References