¿Cuál es la varianza de la varianza?

Para introducir el bootstrap sugiero aquí el problema de estimar la precisión de la varianza muestral. Este ejemplo tiene alguna ventaja respecto al de estimar la precisión de la media muestral ya que las fórmulas teóricas para la varianza no son tan conocidas para los estudiantes como las de la media, lo que creo que hace que valoren mejor el interés del bootstrap. Concluyo con una breve reflexión sobre la docencia en asignaturas de introducción a la estadística.

Tres varianzas distintas

La varianza \(\sigma^2\) de una distribución de probabilidad \(F\) (la primera de las varianzas) mide la dispersión de los valores aleatorios generados a partir de esa distribución (o población). Es un parámetro importante porque da una medida de lo heterogénea que es la población e interviene de manera decisiva en cualquier inferencia que hagamos sobre, por ejemplo, la media poblacional.

Si por ejemplo generamos una muestra de 200 datos normales con media 0 y varianza \(\sigma^2 = (1.4)^2 = 1.96\), podemos estimar \(\sigma^2\) calculando la varianza muestral (la segunda de las varianzas) \[S^2 = \frac{1}{n-1}\sum_{i=1}^n (X_i - \bar{X})^2.\] En la siguiente simulación vemos que el valor de \(S^2\) obtenido es parecido a \(\sigma^2\), pero no igual:

set.seed(123)  # para reproducir los resultados
n <- 200
sigma <- 1.4
sigma^2
#> [1] 1.96
datos <- rnorm(n, sd = sigma)  # genera la muestra
var(datos)
#> [1] 1.743519

De hecho, si repetimos el experimento anterior muchas veces obtenemos muchos valores distintos de \(S^2\) que tienen su propia dispersión, determinada por el valor de \(\mbox{Var}(S^2)\) (la tercera de las varianzas):

R <- 10000 # número de réplicas
varianzas <- replicate(R, var(rnorm(n, sd = sigma)))
var(varianzas)
#> [1] 0.03785888

Esta tercera de las varianzas es importante porque cuantifica hasta qué punto podemos estar seguros de que \(S^2\approx\sigma^2\). De hecho, la desigualdad de Chebychev (dado que \(S^2\) es insesgado) implica que para cualquier margen de error \(\epsilon>0\), \[\mbox{P}(|S^2-\sigma^2|\geq \epsilon) \leq\frac{\mbox{Var}(S^2)}{\epsilon^2}.\]

Luego si \(\mbox{Var}(S^2)\approx 0\) la probabilidad de que \(S^2\) difiera mucho de \(\sigma^2\) será muy baja. ¿Qué es lo que explica el valor 0.038 que hemos obtenido en la simulación? ¿Cómo podemos aproximarlo a partir de la muestra sin conocer la distribución de la que proceden los datos? En esta entrada se muestra cómo valorar la precisión de \(S^2\) como estimador de \(\sigma^2\) sin recurrir a fórmulas matemáticas o teoría estadística.

Si asumimos que la distribución de los datos es normal

Si podemos suponer que los datos vienen de una distribución normal (con la media y la varianza desconocidas), algo muy sencillo que podemos hacer es repetir la simulación anterior pero sustituyendo los valores de \(\mu\) y \(\sigma\) por los estimadores \(\bar{x}\) y \(s\) calculados a partir de los datos. Este procedimiento se denomina bootstrap paramétrico:

R <- 10000 # número de réplicas
varianzas <- replicate(R, var(rnorm(n, mean = mean(datos), sd = sd(datos))))
var(varianzas)
#> [1] 0.03011815

Este método evita tener que conocer la teoría clásica estadística, de acuerdo con la cual si los datos son normales la distribución de \(S^2\) es una \(\chi^2_{n-1}\) (salvo constantes): \[\frac{(n-1)S^2}{\sigma^2} \equiv \chi^2_{n-1}.\] A partir del resultado anterior se deduce la relación teórica entre las tres varianzas del apartado anterior: \[\mbox{Var}(S^2) = \frac{2\sigma^4}{n-1}. \tag{1}\] La ecuación (1) explica el valor que obtuvimos en la simulación inicial:

2*sigma^4/(n-1)
#> [1] 0.03860905

De hecho, a partir de (1), el estimador más natural de la varianza de \(S^2\) es \(2S^4/(n-1)\)

2*var(datos)^2/(n-1)
#> [1] 0.03055134

Este valor es prácticamente idéntico al que se obtiene con bootstrap paramétrico. No es extraño ya que, si uno se para a pensarlo, el valor obtenido con lo que he llamado bootstrap paramétrico es en realidad una aproximación de Monte Carlo de \(2S^4/(n-1)\), que en realidad es el estimador bootstrap paramétrico ideal, el que se obtiene al sustituir la verdadera distribución \(\mbox{N}(\mu,\sigma^2)\) por la distribución estimada \(\mbox{N}(\bar{x},s^2)\). La ventaja de la aproximación de Monte Carlo es que no requiere saber nada sobre la \(\chi^2\) ni conocer la fórmula (1), y esto hace que la idea del bootstrap combinada con el procedimiento de Monte Carlo sea tan potente.

Si no sabemos nada sobre la distribución de los datos

Si no tenemos ni idea de la distribución \(F\) de la que proceden los datos, podemos aplicar una versión más radical del bootstrap llamada bootstrap no paramétrico. Se repite la simulación pero ahora los datos se extraen de la distribución empírica \(F_n\) (la que asigna probabilidad \(1/n\) a cada uno de los datos originales \(X_1,\ldots,X_n\)). Este método es razonable al menos para tamaños muestrales grandes ya que en este caso el teorema de Glivenko-Cantelli garantiza que \(F_n\approx F\).

Veamos qué pasa:

R <- 10000 # número de réplicas
varianzas <- replicate(R, var(sample(datos, replace = TRUE)))
var(varianzas)
#> [1] 0.03349739

Mediante el procedimiento anterior, extremadamente simple y sin suponer nada sobre \(F\) ni conocer ninguna teoría, hemos podido valorar la precisión de \(S^2\). El valor obtenido no está lejos del correcto, 0.0386.

La relación general entre las tres varianzas

Para los/as curiosos/as, una bonita fórmula que viene en algunos libros es la equivalente a (1), pero para datos que no son normales. La relación general entre las tres varianzas es: \[\mbox{Var}(S^2) = \frac{\sigma^4}{n}\left(\kappa -1 + \frac{2}{n-1}\right). \tag{2}\] En (2), \(\kappa\) representa el coeficiente de curtosis. En el caso de la distribución normal \(\kappa=3\), lo que muestra que (1) es un caso particular de (2). Si prescindimos de términos de orden \(n^{-2}\), tenemos la aproximación \[\mbox{Var}(S^2) \approx \frac{\sigma^4(\kappa-1)}{n}\] válida para cualquier distribución tal que \(\kappa<\infty\). Esta aproximación es la misma que se deduce de la aplicación del teorema central del límite a \(S^2\) para obtener su distribución asintótica.

Una breve reflexión sobre la docencia

Una comprensión del bootstrap como método para valorar la precisión de un estimador es mucho más formativa conceptualmente y se acerca más a la esencia de la inferencia estadística que conocer las propiedades de la distribución \(\chi^2\) o de la \(t\) de Student. Además, su rango de aplicaciones es mucho más amplio. Creo que se debería explicar en los cursos introductorios y se le debería dar la importancia que merece de acuerdo con su impacto en la práctica estadística de las últimas décadas. El problema es la inercia. Por una parte en los docentes, ya que suele ser difícil cambiar los contenidos de asignaturas que a veces son impartidas por muchas personas diferentes. Pero también la inercia en los investigadores aplicados a quienes a menudo ahuyenta cualquier desviación de los métodos que usan tradicionalmente.