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 harrell1982). 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
- [HoSM]
NIST/SEMATECH e-Handbook of Statistical Methods, http://www.itl.nist.gov/div898/handbook/, October 6, 2020. https://doi.org/10.18434/M32189 - [Harrell1982]
Harrell, F.E. and Davis, C.E., 1982. A new distribution-free quantile estimator. Biometrika, 69(3), pp.635-640.
https://pdfs.semanticscholar.org/1a48/9bb74293753023c5bb6bff8e41e8fe68060f.pdf