Cómo calcular la mediana sin tener que ordenar los datos

Se explica aquí cómo aprovechar que la mediana es, en cierto sentido, un valor extremo, para calcularla sin tener que ordenar los datos. También se comenta por qué esta observación tiene implicaciones relevantes en lo que se conoce como regresión de cuantiles.

La mediana resuelve un problema de programación lineal

Dice Aristóteles en la Ética a Nicómaco que “la virtud es un término medio, pero con respecto a lo mejor y al bien, es un extremo”, e igualmente en estadística las medidas de la posición central de un conjunto de datos se pueden considerar extremos, en el sentido de ser las que minimizan alguna medida de la magnitud de las desviaciones de los datos a una posición central dada. Por ejemplo, la mediana es el valor que minimiza la suma de desviaciones absolutas. Si definimos la función objetivo $$f(t)=\sum _{i=1}^n|x_i-t|,$$

entonces las medianas son los mínimos globales de $f$. Se puede consultar esta entrada del blog para una demostración sencilla de esta propiedad.

Desde este punto de vista, en lugar de ordenar los datos, el cálculo de la mediana requiere resolver un problema de optimización, que por un lado es fácil porque es convexo (no hay mínimos locales) pero por otra parte en la función objetivo hay un valor absoluto, lo que impide el uso de métodos que requieran derivadas.

La observación clave es que el problema de optimización que hay que resolver para calcular la mediana se puede escribir como un problema de programación lineal, por lo que basta usar cualquier método estándar para este tipo de problemas (por ejemplo, el famoso algoritmo del simplex).

Veamos cómo escribir el problema de forma apropiada:

Para $i=1,\dots ,n$ definimos las variables auxiliares $$x_i^{+}=max\{x_i-t,0\}\ge 0,\text{y}{x}_i^{-}=max\{-(x_i-t),0\}\ge 0.$$ Es muy fácil ver que $|x_i-t|=x_i^{+}+x_i^{-}$ y que $x_i-t=x_i^{+}-x_i^{-}$. En función de las nuevas variables lo que queremos es minimizar en $x_i^{+}$, $x_i^{-}$ y $t$ el objetivo $$g(x_1^{+},\dots ,x_n^{+},x_1^{-},\dots ,x_n^{-},t)=\sum _{i=1}^n{x}_i^{+}+\sum _{i=1}^n{x}_i^{-},$$ sujeto a

$$\begin{array}{rl}{x}_i^{+}-x_i^{-}+t& =x_i,i=1,\dots ,n,\\ x_i^{+}& \ge 0,i=1,\dots ,n,\\ x_i^{-}& \ge 0,i=1,\dots ,n.\end{array}$$

El valor de la mediana es el valor óptimo $t^{\ast }$ del problema anterior, que es un problema de optimización lineal (tanto la función objetivo como las restricciones están definidas a partir de funciones lineales de las variables).

En forma estándar

Recordemos que un problema de programación lineal en forma estándar es aquel en el que queremos encontrar el vector $y$ que resuelve $$\begin{array}{}\text{(1)}& min{c}^{\prime }y\text{sujeto a}Ay=b,y\ge 0,\end{array}$$ para una matriz $A$ y unos vectores $c$ y $b$ de las dimensiones adecuadas. He usado la notación habitual en los libros que tratan de estos problemas (en particular, $y\ge 0$ significa que todas las coordenadas de $y$ son no negativas).

En función de las nuevas variables $x_i^{+}$, $x_i^{-}$ y $t$ ya tenemos un problema de optimización lineal, pero aún no está en forma estándar porque la variable $t$ puede ser negativa. Pero esta dificultad se resuelve fácilmente mediante uno de los trucos habituales que vienen en los libros de programación lineal: escribir una variable cuyo signo es arbitrario como diferencia de dos variables no negativas.

Obsérvese que $t=t^{+}-t^{-}$, donde $t^{+}=max\{t,0\}\ge 0$ y $t^{-}=max\{-t,0\}\ge 0$. Por lo tanto, nuestro problema equivale a minimizar en $x_i^{+}$, $x_i^{-}$, $t^{+}$ y $t^{-}$ el objetivo $$g(x_1^{+},\dots ,x_n^{+},x_1^{-},\dots ,x_n^{-},t^{+},t^{-})=\sum _{i=1}^n{x}_i^{+}+\sum _{i=1}^n{x}_i^{-},$$ sujeto a

$$\begin{array}{rl}{x}_i^{+}-x_i^{-}+t& =x_i,i=1,\dots ,n,\\ x_i^{+}& \ge 0,i=1,\dots ,n,\\ x_i^{-}& \ge 0,i=1,\dots ,n,\\ t^{+}& \ge 0,\\ t^{-}& \ge 0.\end{array}$$

Cómo resolverlo con R

Hace unos años escribí una guía sobre cómo resolver problemas de optimización lineal con R en la que se explicaban las principales características del paquete lpSolve. Aquí lo voy a usar para calcular la mediana. Para entender mejor la función hay que escribir nuestro problema en la forma (1). Para ello, definimos $y=(x_1^{+},\dots ,x_n^{+},x_1^{-},\dots ,x_n^{-},t^{+},t^{-}{)}^{\prime }\in {\mathbb{R}}^{2n+2}$, $c=(1,\dots ,1,0,0)\in {\mathbb{R}}^{2n+2}$, $b=(x_1,\dots ,x_n{)}^{\prime }\in {\mathbb{R}}^n$ y $A=(\mathbb{I},-\mathbb{I},u,-u)\in {\mathbb{R}}^{n\times 2n}$, donde $\mathbb{I}$ es la identidad de orden $n$ y $u=(1,\dots ,1{)}^{\prime }\in {\mathbb{R}}^n$. Es decir, nuestro problema al final tiene $2n+2$ variables y $n$ restricciones, donde $n$ es el tamaño muestral, además de la restricción de no negatividad de todas las variables.

library(lpSolve)

mediana.lp <- function(y){
  # Calcula el valor t que minimiza sum_i |y_i-t|
  n <- length(y)
  cvec <-c(rep(1,2*n), 0, 0)
  b <- y
  unos <- rep(1, n)
  A <- cbind(diag(unos), diag(-unos), unos, -unos)
  dir <- rep('=', n)
  
  optimo <- lp('min', cvec, A, dir, b)$solution
  return(optimo[2*n+1] - optimo[2*n+2])
}

# Ejemplos
y <- c(2, 3, 4, 5, 6)
mediana.lp(y)
#> [1] 4

set.seed(123)
y <- rnorm(500)
mediana.lp(y)
#> [1] 0.02045071

Regresión de cuantiles

Este ejemplo lo usé como ejercicio de programación lineal cuando impartí la asignatura Investigación Operativa en la UAM, y puede parecer que es una forma muy rebuscada de calcular la mediana, pero su interés es mayor de lo que parece. El modelo de regresión de cuantiles lineal consiste en suponer que la mediana (u otros cuantiles) de la distribución de una variable respuesta $y$ se pueden poner, salvo por un término de error aleatorio, como función lineal de un vector $x$ de variables regresoras: $$\text{Mediana}(y|x)={\beta }^{\prime }x+\epsilon .$$ Una forma de estimar el vector $\beta$ de coeficientes a partir de una muestra $(y_1,x_1),\dots ,(y_n,x_n)$ consiste en encontrar el vector $\beta$ que minimiza $$f(\beta )=\sum _{i=1}^n|y_i-{\beta }^{\prime }{x}_i|.$$ En este modelo el criterio de mínimos cuadrados habitual se sustituye por el mínimo de los valores absolutos ya que lo que queremos es estimar la mediana de $y$ condicionada a $x$ en lugar de la media.

Por un razonamiento similar al que hemos descrito en esta entrada, el vector $\beta$ de parámetros se puede estimar resolviendo un problema de programación lineal. Además, si añadimos unos pesos $w_1,\dots ,w_n$ adecuados y optimizamos la función $$f(\beta )=\sum _{i=1}^n{w}_i|y_i-{\beta }^{\prime }{x}_i|$$ la técnica descrita no solo es válida para estimar la mediana condicionada sino cualquier cuantil de interés. Así, resolviendo problemas de programación lineal se puede obtener una descripción bastante completa de la distribución de $y$ condicionada a $x$.