Capítulo 6 Estimación puntual

El objetivo inmediato de la estimación puntual es el de obtener estimaciones razonables de parámetros de alguna distribución (p.e., la media) desconocida —o parcialmente desconocida— de la que se tiene una muestra, así como las correspondientes estimaciones del error cometido. Un ejemplo paradigmático de la estimación puntual es el de la estimación de la tasa de paro, que es una magnitud desconocida \(\theta\). Para ello, sobre la población española se realiza una encuesta, la EPA, con la que se construye un estimador \(\hat{\theta}\) de \(\theta\).

No obstante, los principios y métodos explicados aquí sirven igualmente para la estimación de parámetros en modelos y son, precisamente, sirven precisamente de base teórica para el ajuste de modelos en ciencia de datos.

En el capítulo se van a estudiar primero dos métodos clásicos de la estimación puntual y se discutirán los correspondientes problemas de la variabilidad de las estimaciones y la magnitud del error cometido y se terminará discutiendo —y tratando de tender puentes con la discusión previa— la técnica ḿas habitual en ciencia de datos para resolver el mismo problema, el de la minimización de determinadas funciones de pérdida.

Queda fuera de este capítulo, y no porque sea irrelevante sino porque se cubrirá en uno posterior específico, la perspectiva bayesiana del problema.

6.1 Estimación puntual clásica

En esta sección se van a presentar dos métodos clásicos de estimación puntual de parámetros: el de los momentos y el de la máxima verosimilitud.

6.1.1 El método de los momentos

El método de los momentos es posiblemente el primero en formalizarse para resolver el problema de la estimación puntual. Se atribuye a Karl Pearson (Pearson 1894) aunque parece improbable que no existan precedentes de su uso en casos simples.

Podemos ilustrarlo con un ejemplo. Podemos suponer que se tiene una muestra de una variable aleatoria con distribución \(X \sim \Gamma(a, b)\) y el problema consiste en estimar los parámetros desconocidos \(a\) y \(b\). Como sabemos que \(E(X) = a/b\) y \(\sigma^2(X) = a / b^2\), podemos calcular la media y la varianza muestral, plantear las ecuaciones correspondientes y despejar nuestras estimaciones \(\hat{a}\) y \(\hat{b}\).

En lugar de resolver este problema analíticamente, podemos plantearlo en R directamente. En primer lugar, extraeremos una muestra de una distribución gamma con parámetros 3 y 4 y calcularemos la media y la varianza de esta muestra:

Efectivamente, podemos comprobar cómo la media y la varianza son aproximadamente iguales a las esperadas:

Despejando de las ecuaciones anteriores, se obtiene \(b = E(X) / \sigma^2(X)\) y \(a = E^2(X)/\sigma^2(X)\), por lo que nuestras estimaciones de los parámetros \(a\) y \(b\) son:

## [1] 3.002717 4.101987

Cabe preguntarse, muy razonablemente, cómo saber si, dada una muestra de números, estos proceden de una distribución gamma. Hay varios procedimientos que se pueden emplear para tratar de dirimir esa cuestión. En la práctica, si uno considera que la distribución gamma es razonable para unos datos es que, efectivamente, hay razones para hacerlo: por ejemplo, que proceden de un sistema en el que históricamente se ha comprobado que los datos generados tienen ese tipo de distribución. Existen métodos informales para, de alguna manera, comprobar que los datos siguen esa distribución, como comparar los histogramas con la función de densidad teórica o usar la función qqplot como se hizo en el capítulo dedicado a la estadística descriptiva. Finalmente, hay métodos de verificación más formales y basados en los test de hipótesis que se mostrarán más adelante.

La siguiente gráfica muestra cómo los valores estimados de \(a\) y \(b\) se aproximan a sus verdaderos valores según aumenta el número de elementos en la muestra:

Los dos objetivos de la estimación puntual de parámetros son que las estimaciones converjan a los valores reales de los parámetros y que, por otro lado, la dispersión —o el error, que corresponde al tamaño de las cajas en el diagrama anterior— alrededor de dicho valor real sea la menor posible. El método de los momentos, en ese sentido, es inferior al de la máxima verosimilitud —que se mostrará más adelante— y por tanto, no suele usarse mucho en la práctica. Sin embargo, por varios motivos, entre ellos el interés por analizar datos muy voluminosos, ha supuesto un renacer parcial del método de los momentos por motivos como los que se explican en el resumen de (Perry 2015):

El artículo propone un procedimiento basado en el método de los momentos […]. El procedimiento sacrifica eficiencia estadística en aras de una mayor eficiencia computacional. […] Al aplicarse a sistemas de recomendación a gran escala y compararse con la estimación por máxima verosimilitud, este método proporciona predicciones comparables pero reduce el tiempo del ajuste de los modelos de horas a minutos.

Para la distribución normal \(N(\mu, \sigma)\), \(E(X) = \mu\) y \(\sigma^2(X) = \sigma^2\), por lo que en tal caso la estimación de los parámetros por el método de los momentos es trivial. En cambio, en otras situaciones, aplicar el método de los momentos conlleva resolver una serie de ecuaciones no lineales numéricamente.

6.1.2 Estimación por máxima verosimilitud

La estimación por máxima verosimilitud está relacionada con el concepto, sumamente intuitivo, de valor más verosímil, que podemos ilustrar con el siguiente ejemplo. Supongamos que hay dos urnas; la primera contiene cinco bolas blancas y un negra; la segunda, cinco bolas negras y una blanca. Alguien extrae de una de ellas (y no nos dice cuál) una bola que resulta ser negra. Parece que la opción más verosímil es la segunda, la que contiene cinco bolas negras.

Esta idea intuitiva se puede formalizar para luego poder ser extendida. En efecto, la probabilidad de que la bola sea negra de haber sido extraída de la primera urna es \(P(N \; | \; U = 1) = 1/6\). La probabilidad correspondiente a la segunda urna es \(P(N \; | \; U = 2) = 5/6\). Por ese motivo nos decantamos por \(U = 2\).

Ya se habían tratado problemas como este en capítulos previos. En particular, al discutir el teorema de Bayes. La diferencia esencial entre aquella aproximación y esta es que entonces, la decisión se basaba en \(P(U = x \; | \; N)\) y no en \(P(N \; | \; U = x)\). Esa sutil diferencia, expresamente manifiesta aquí, es la que, como se verá, separa la estadística frecuentista de la bayesiana.

Hay que hacer observar, en todo caso, que en este problema en particular, la decisión es la misma dado que, en principio, las probabilidades a priori de cada urna son las mismas.

El ejemplo anterior se puede extender a una situación más general donde se tienen unos datos extraídos de una distribución \(P_\theta\), donde \(\theta\) es un parámetro desconocido, ese del que se quiere obtener una estimación. Igual que con las urnas se ha construido \(P(N \; | \; U = x)\), es decir, la probabilidad de los datos (una bola negra, respresentada por \(N\)) en función de la urna (\(U = 1\) o \(U=2\)), se puede construir la función \(P(\text{datos} \; | \; \text{parámetros}) = P(D \; | \; \theta)\) y plantear, por analogía, que de tener que optar por algún valor de \(\theta\) es razonable hacerlo por aquel que maximiza esa expresión. Dicha expresión se conoce como función de verosimilitud y se representa así:

\[L(\theta) = L(\theta, X) = P(X \; | \; \theta).\]

En el siguiente ejemplo, se va a usar la estimación por máxima verosimilitud para estimar el parámetro desconocido en un problema de lanzamiento de monedas: se tira una 100 veces y se obtienen 60 caras. En este caso, interesa estimar el parámetro \(\theta\) de una distribución binomial, la probabilidad de cara. La función de verosimilitud es

\[L(\theta) = P(X = 100 \; | \; \theta) = {{100}\choose{60}} \theta^{60} (1-\theta)^{40}\]

que es la que hay que maximizar. Es evidente que maximizar esa expresión equivale a maximizar \(\theta^{60} (1-\theta)^{40}\) porque el factor faltante no depende de \(\theta\). También lo es que maximizar esa expresión equivale a maximizar su logaritmo, \(60 \log \theta + 40 \log (1-\theta)\). Ahora es sencillo derivar con respecto a \(\theta\), igualar a cero, despejar y obtener la estimación \(\hat{\theta} = 0.6\).

De hecho, es habitual que la verosimilitud tenga forma de producto de términos (por ejemplo, cuando las observaciones, como en este caso, son independendientes) y es muy habitual operar con el logaritmo de la verosimilitud en lugar de con esa función directamente.

Existe otro motivo más pedestre por el que tomar logaritmos: el de la precisión numérica cuando las operaciones se realizan con un ordenador. En el ejemplo anterior, \(\theta^{60}\), cuando \(\theta = 0.5\), vale 8.673617e-19.

También es posible resolver analíticamente otro caso importante: \(X\) sigue una distribución normal \(N(\mu, \sigma)\) y se obtiene una muestra \(x_1, \dots, x_n\) de \(X\), entonces se puede probar que la estimación por máxima verosimilitud de la media, \(\hat{\mu}\) no es otra cosa que \(\frac{1}{n} \sum_i x_i\), la media muestral de las observaciones. En efecto —por una vez merece la pena explicitar las operaciones necesarias para obtener este resultado—, partiendo de la verosimilitud

\[L(\mu, \sigma) = P(x_1, \dots, x_n \; | \; \mu, \sigma) = \prod_i \frac{1}{\sqrt{2 \pi \sigma^2}} \exp\left(-\frac{1}{2} \left(\frac{x_i - \mu}{\sigma}\right)^2\right),\]

podemos tomar logaritmos para obtener la expresión

\[-\frac{n}{2} \log(2 \pi \sigma^2) - \frac{1}{2} \sum_i \left(\frac{(x_i - \mu)^2}{\sigma^2} \right),\]

cuya maximización (respecto a \(\mu\)) equivale a la mimimización de

\[\sum_i (x_i - \mu)^2.\]

Esa última expresión es una suma cuadrática y se puede probar que su mínimo es precisamente \(\sum_i x_i / n\), la media de las observaciones. De hecho, en este caso, las estimaciones por máxima verosilitud, con una pérdida cuadrática y por el método de los momentos, coinciden.

Hay que tener en cuenta de que el método de la máxima verosimilitud se aplica no solo en casos simples como los discutidos hasta este punto sino también en situaciones más complejas. De hecho, la práctica totalidad de los modelos estadísticos que se discuten en los capítulos siguientes se ajustan —es decir, se realizan estimaciones puntuales de los parámetros y coeficientes que intervienen en ellos— mediante el procedimiento de la máxima verosimilitud. La maximización de estas versomimilitudes, en ocasiones funciones muy complejas de muchos parámetros, se realiza necesariamente usando métodos computacionales muy intensivos.

Por ejemplo, en ciencia de datos, dada una secuencia de pares de valores \((y_i, x_i)\) (donde \(x_i\) puede ser un vector) se postula más o menos implícitamente que \(y_i \sim N(f_\theta(x_i), \sigma)\), donde \(f_\theta\) es una función más o menos compleja —p.e., una red neuronal— y utilizando una variante de la demostración de más arriba, se concluye que sus parámetros \(\hat{\theta}\) óptimos son los que minimizan

\[\sum_i \left( y_i - f_\theta(x_i)\right)^2,\]

una expresión que se suele postular sin mayor justificación en muchos textos de la materia.

Pero la minimización de una suma cuadrática es consecuencia de asumir que los errores, es decir, las diferencias \(y_i - f_\theta(x_i)\) tienen una distribución normal con varianza \(\sigma\) constante, una hipótesis que no siempre es sostenible.

6.2 Variabilidad de las estimaciones

La cuestión que se plantea en esta sección es cómo variarían las estimaciones de repetirse el proceso tal vez con otros datos pero producidos por el mismo fenómeno aleatorio.

Un experimento mental: el Instituto Nacional de Estadística español realiza una encuesta trimestral para medir la tasa de paro. Cada tres meses, aplica un determinado procedimiento de muestreo, realiza una serie de entrevistas personales, postprocesa los datos y publica en su portal un número (p.e., 15.4%). Esa se considera la tasa oficial de paro en España. El experimento mental es el siguiente: imaginad que 10 equipos independientes del INE aplicasen exactamente los mismos métodos, salvo que a la hora de seleccionar los hogares entrevistados, a cada uno de ellos les correspondiese por azar hogares distintos. A la hora de publicar los resultados, ¿coincidirían todos ellos? ¿Cómo serían de parecidos? ¿Qué, en definitiva, podríamos aprender en el caso de que fuesen muy similares (o diferentes)?

En esta sección se van a describir dos procedimientos habituales para estimar la variabilidad de las estimaciones \(\hat{\theta}\) y para representarla, particularmente a través de los llamados intervalos de confianza.

De nuevo, los intervalos de confianza son resúmenes de ancho de banda estrecho para representar la variabilidad de los estimadores. Casi siempre hay opciones mejores para representarla, incluidas las gráficas.

En secciones posteriores se justificará cómo se relaciona la variabilidad de las estimaciones con, entre otras cosas, el error cometido con ellas.

6.2.1 El método plug-in

El método del plug-in puede aplicarse solo cuando las distribuciones de partida de los datos son conocidas. Aquí, conocidas no significa exactamente que sean una de las que aparecen con nombre en los libros: pueden ser producto de algún tipo de proceso aleatorio conocido y parametrizado. Se supondrá, en particular aunque sin pérdida de generalidad, que los datos proceden de una distribución con función probabilidad \(F_{\theta}\).

A partir de dicha distribución se ha obtenido una muestra \(x_1, \dots, x_n\) y, a partir de ella, un estimador \(\hat{\theta}\). Representemos el proceso completo por medio del operador \(G\) tal que

\[G(\theta) \longmapsto \hat{\theta},\]

es decir, el que nos ha proporcionado la estimación. Para obtener una muestra de ellas y analizar su distribución, podríamos ponerlo en marcha otra vez (pasando, en cada iteración, por una muestra distinta del mismo tamaño \(n\) que la original). Sin embargo, eso no es posible dado que no conocemos \(\theta\).

El método del plug-in consiste en utilizar \(G\) operando sobre \(\hat{\theta}\) en lugar de sobre \(\theta\).

Por ejemplo, si se tiene una muestra de 100 observaciones \(x_i\) de una distribución normal \(N(\mu, \sigma)\),

podemos estimar su media como

y para analizar su variabilidad, generar otras muestras usando como reemplazos de mu y sigma, desconocidos, sus estimaciones muestrales est_mu y

Es decir, se puede generar una muestra de 10000 valores posibles de est_mu haciendo

La gráfica que sigue muestra la distribución de esos valores junto con el valor estimado, marcado en rojo:

A partir de esa distribución de valores se pueden calcular estadíticos relevantes suyos como su varianza, sus cuantiles, etc. Si bien es cierto que en ocasiones es posible evitar el paso por el remuestreo para estimarlos. Efectivamente, en ciertos casos de libro, es posible calcularlos directamente. En este caso en particular, sabemos que la distribución de las muestras es normal con parámetros est_mu y est_sigma, por lo que se pueden calcular directamente:

  • La desviación estándar, que es est_sigma.
  • Los cuantiles q haciendo qnorm(q, est_mu, est_sigma).
  • Los intervalos de confianza, que son parejas de cuantiles (típicamente, del .025 al .975).

Otro ejemplo: si se tira una moneda al aire 100 veces y salen 60 caras, el intervalo de confianza teórico al 95% para la media sería

## [1] 0.61 0.79

donde p es, precisamente, el valor desconocido que se quiere estimar. Por eso, como \(\hat{\theta}_n = 0.6\) (que es un valor obtenido, por ejemplo, por máxima verosimilitud), en virtud del principio del plug-in, el intervalo de confianza típico al 95% sería

## [1] 0.50 0.69

Cabe, lógicamente, preguntarse por qué funciona el método del plug-in en lugar de, sencillamente, no hacerlo. Como se verá más adelante, cabe esperar que \(\hat{\theta}\) esté próximo al valor real \(\theta\), por lo que las distribuciones \(F_\theta\) y \(F_\hat{\theta}\) deberían ser parecidas (y más, al aumentar el tamaño de la muestra). Esto se puede ilustrar comparando la función de probabilidad \(F_\theta\) con diversas \(F_\hat{\theta}\) construidas a partir de muestras cada vez más grandes, como en el siguiente gráfico, donde se muestra cómo se aproximan las funciones de probabilidad plug-in (trazo gris claro) a la original (trazo rojo), una normal estándar, conforme aumenta el número de observaciones en la muestra:

6.2.2 Intervalos de confianza mediante remuestreos

La técnica de los remuestreos (o bootstrap) apareció a finales de los años 70 y es más intensiva computacionalmente que las anteriores. Está relacionada con la técnica del plug-in en el sentido de que si este considera \(F_\hat{\theta}\) como una a proximación a \(F_\theta\) para construir intervalos de confianza, el remuestreo, de alguna manera, reemplaza \(F_\theta\) por su aproximación \(\hat{F}_\theta\).

En concreto, la técnica del remuestreo simula otras muestras posibles de la distribución subyacente realizando remuestreos (con reemplazamiento) de la muestra original. Por ejemplo, si la muestra original está contenida en la variable x de R y el estimador buscado \(\hat{\theta}\) está dado por la media muestral mean(x), entonces se puede obtener una distribución aproximada de este estimador repitiendo la operación

tantas veces como sea preciso. Finalmente, a partir de esas muestras se pueden calcular los correspondientes cuantiles y a partir de estos, los intervalos de confianza.

Al aplicar el método del plug-in, todavía era necesario conocer la distribución de partida para insertar en ella el parámetro estimado puntualmente. Ahora ya no es necesario: solo se utiliza la muestra, independientemente de la naturaleza de la distribución de partida.

Este procedimiento está justificado en la siguiente observación (que, de hecho, es un teorema, el de Glivenko-Cantelli): que la función de distribución empírica \(F_n\) construida a partir de muestras \(x_1, \dots, x_n\) de una distribución \(F\) converge a \(F\) conforme aumenta el tamaño muestral. Gráficamente, si \(F\) es la distribución lognormal,

En el gráfico anterior se aprecia cómo la función de probabilidad empírica con n = 1000 en este caso es prácticamente indistinguible de la teórica. Es decir, muestrear una u otra debería ser prácticamente equivalente a efectos prácticos.

El siguiente gráfico compara la distribución teórica de la media de una muestra de 100, 500 y 1000 observaciones de una normal estándar con las estimadas tanto por el método del plug-in, como mediante remuestreos mostrando cuatro resultados simulados por tamaño muestral:

La ventaja de los remuestreos es —como se ha indicado más arriba— que permiten estimar la variabilidad de parámetros incluso cuando estos no tienen una distribución conocida o evidente. Un ejemplo extraído de (Efron and Hastie 2016) tiene que ver con la regresión local, lowess, introducida en el capítulo dedicado a la estadística descriptiva. El siguiente código —basado en un conjunto de datos, kidney, disponible en la página web del libro anterior— y el gráfico que generan ilustran esta técnica:

6.3 Convergencia y error de las estimaciones

En esta sección se estudia una cuestión muy importante relacionada con la de las anteriores: si esas estimaciones son fiables. En particular, si están próximas a los valores que intentan estimar y en estudiar su variabilidad. El planteamiento del problema está muy bien ilustrado, que prácticamente se ha convertido en un meme obligatorio al discutir estas cuestiones:

6.3.1 Convergencia de las estimaciones

En las secciones anteriores se han planteado y justificado dos procedimientos para obtener estimaciones puntuales de parámetros de distribuciones. Pero surge una pregunta natural: ¿funcionan verdaderamente? Más concretamente, si hemos estimado \(\hat{\theta}\), ¿tenemos garantías de que esa estimación esté próxima al valor verdadero \(\theta\)?

Existe una serie de teoremas, los llamados de convergencia, que estudian qué sucede si construimos el estimador \(\hat{\theta}_n\) a partir de muestras crecientes en tamaño, \(n\), de la distribución subyacente. Una de las preguntas que intentan responder es si como se espera (o en qué casos), \(\hat{\theta}_n\) tiende a \(\theta\), i.e., si al crecer el tamaño de la muestra, el valor estimado se aproxima progresivamente al real.

Estos teoremas se preocupan también del sentido que tiene esa convergencia. En concreto, si lo es en probabilidad o casi seguro. Esta es una discusión muy técnica que excede el alcance de este libro.

En realidad, en algunos de los casos vistos anteriormente, el estimador \(\hat{\theta}\) no es otra cosa que una proporción o, más en general, una media. En esos casos (y bajo ciertas condiciones no particularmente estrictas), la ley de los grandes números garantiza la convergencia. Los teoremas de convergencia arriba mencionados permiten extender estos resultados a (prácticamente) cualquier estimador \(\hat{\theta}\) de los considerados más arriba.

El lector interesado, de todos modos, podrá encontrar demostraciones accesibles de la convergencia del estimador de máxima verosimilitud —en el caso unidimensional y en condiciones sumamente generales— en, por ejemplo, la sección 4.2 de (Efron and Hastie 2016).

6.3.2 Error de las estimaciones e intervalos de confianza

En la práctica, no solo interesa tener garantías de que \(\hat{\theta}_n\) tiende a \(\theta\): también se busca alguna indicación de la distancia entre ambos valores, es decir, del error cometido.

Por ejemplo, en la ficha técnica de las encuestas electorales se encuentran expresiones del tipo:

Error de la muestra: \(\pm 3.5\)% para un nivel de confianza del 95%.

Eso significa que alrededor de la estimación de la intención de voto de un partido —que no deja de ser una estimación puntual— existe un rango de \(\pm 3.5\) puntos porcentuales dentro de los cuales se espera esté contenida la proporción de apoyo real a ese partido con una confianza —piénsese: probabilidad— del 95%. Las encuestas con mayor número de encuestados —y, por lo tanto, más caras— tienen un error muestral menor, es decir, un rango de variabilidad más estrecho; dicho de otra manera, aspiran a una mayor precisión.

Lo que proporciona la teoría es esencialmente una distribución —usualmente no conocida o solo conocida aproximadamente— de las diferencias entre el estimador y el valor real. En la práctica, estas distribuciones se representan en rangos de valores —como en el ejemplo anterior— alredededor del estimador: los intervalos de confianza.

En secciones previas de este capítulo se han introducido los intervalos de confianza como rangos de variabilidad esperada de los estimadores. El objetivo de esta sección consiste, esencialmente, en ilustrar cómo ambos conceptos vienen a ser uno mismo.

6.3.2.1 Error de las estimaciones

Igual que existen resultados teóricos que extienden la ley de los grandes números y que garantizan la convergencia de \(\hat{\theta}_n\) a \(\theta_0\), las extensiones correspondientes del teorema central del límite concluyen que los ejemplos de la sección anterior son casos particulares de una propiedad bastante habitual: es frecuente que \(\hat{\theta}_n\ - \theta\), al menos para valores de \(n\) bastante altos, tenga una una distribución normal con desviación estándar del orden de \(1 / \sqrt{n}\). O, escrito de otra forma, que

\[\sqrt{n}(\hat{\theta}_n - \theta) \longrightarrow N(0,\sigma),\]

es decir, que las diferencias entre los estimadores y el valor real tienen una distribución aproximadamente normal (para valores altos de \(n\)) con una varianza que depende de un parámetro en principio desconocido \(\sigma\) pero que decrece como \(1/\sqrt{n}\).

Estos resultados no son universales y dependen del problema, de las distribuciones implicadas, etc. Pero es relativamente seguro darlo por bueno en general. Por supuesto, la enumeración y demostración de los teoremas relevantes queda fuera del alcance de estas páginas y, casi seguro, del interés más inmediato de sus lectores.

Eso significa que la distribución de prácticamente cualquier estimador razonable \(\hat{\theta}\), es aproximadamente \(N(\theta, \sigma / \sqrt{n})\). La siguiente gráfica ilustra este resultado. En ella se muestran las distribuciones esperadas de las estimaciones de tres parámetros distintos —la desviación estándar de una normal, la intensidad de una Poisson y la correlación de dos vectores normales— haciendo variar el número de muestras:

Antes de ello conviene subrayar, porque es un hecho casi universal, que distancia entre el estimador y el valor verdadero del parámetro es una variable aleatoria con desviación estándar inversamente proporcional a \(\sqrt{n}\). Es decir, en general, para conseguir el doble de precisión, hace falta una muestra 4 veces mayor. Incrementar el tamaño de las muestras está sujeto a un proceso de incrementos notablemente decrecientes.

Para terminar, merece la pena subrayar un resultado muy importante, la llamada cota de Cramér-Rao. Este resultado asegura que el estimador de máxima verosimilitud el mejor en un sentido muy concreto: en la fórmula anterior, es el que tiene el error, es decir, la \(\sigma\), más pequeña de todos. Es decir, garantiza que los intervalos de confianza son los más estrechos.

Al menos, es el mejor y en ese sentido entre todos los estimadores centrados, sin sesgo. Pueden existir y, de hecho, existen estimadores sesgados que tienen un mejor comportamiento.

Eso quiere decir, por ejemplo, que con datos con una distribución normal, la media —que es el estimador por máxima verosimilitud de \(\mu\)— es superior (en el sentido de que tenderá a estar más próximo a \(\mu\)) que, por ejemplo, la mediana. La salvedad es que eso ocurre cuando —y el autor está tentado de añadir la mentira piadosa y solo cuando— los datos proceden de una distribución normal. En otros casos, la mayoría, la ventaja revierte en favor de la mediana.

6.3.2.2 Intervalos de confianza para la distribución normal

Antes de discutir los intervalos de confianza en general, se van a considerar los relacionados con problemas en los que interviene la distribución normal. Se hace por dos motivos: porque es un caso notable y porque ilustra perfectamente, además, el procedimiento que se aplicará en el caso general.

Si \(X\) es una variable aleatoria normal de media \(\mu\) y desviación estándar \(\sigma\) y se extrae de ella una muestra \(x_i\) de tamaño \(n\), entonces, la media \(\bar{X}\) de la muestra, típicamente, se encuentra cerca de \(\mu\). Esto es así porque ya sabemos que

\[\frac{\bar{X} - \mu}{\sigma / \sqrt{n}}\] tiene distribución normal estándar. Es decir, el 95% de las veces estará comprendida entre \(q_{0.025}\) y \(q_{0.975}\) de dicha distribución, es decir, los valores

## [1] -1.959964  1.959964

Estos valores se redondean habitualmente a $1.96% —como en lo que sigue— o, incluso, a \(\pm 2\).

Lo anterior equivale a decir que, muy frecuentemente, \(\bar{X}\) está dentro del rango

\[\left[\mu - 1.96 \sigma / \sqrt{n}\;, \;\; \mu + 1.96 \sigma / \sqrt{n}\right].\]

O, desplazando el intervalo de confianza convenientemente, que muy frecuentemente, el 95% de las veces aproximadamente (por el redondeo), \(\mu\) estará comprendido en el intervalo

\[\left[\bar{X} - \frac{1.96 \sigma}{\sqrt{n}}\;, \;\; \bar{X} + \frac{1.96 \sigma}{\sqrt{n}}\right].\]

Eso es lo que en este caso se conoce como intervalo de confianza (al 95%) para la media. Las afirmaciones anteriores son comprobables computacionalmente: podemos construir muchos intervalos de confianza como se indica más arriba y estimar la proporción de los casos que contienen el valor \(\mu\). En efecto,

## [1] 0.9517

En este punto, muchos estadísticos observan y son muy celosos de la siguiente sutileza. Así como podemos decir que la probabilidad de que \(\bar{X}\) esté contenido en el intervalo

\[\left[\mu - \frac{1.96 \sigma}{\sqrt{n}}\;, \;\; \mu + \frac{1.96 \sigma}{\sqrt{n}}\right]\]

es del 95%, no se puede decir lo mismo de la probabilidad de que \(\mu\) lo esté en

\[\left[\bar{X} - \frac{1.96 \sigma}{\sqrt{n}}\;, \;\; \bar{X} + \frac{1.96 \sigma}{\sqrt{n}}\right].\]

En este segundo caso, sostienen, \(\mu\) estará dentro o no. Además, el intervalo depende de la muestra obtenida. Lo único que se puede decir, como ilustra la simulación anterior, es que de repetirse el experimento muchas veces, en el 95% de los casos el intervalo contendría el valor verdadero. Pero el intervalo cambia (porque cambia \(\bar{X}\)) en cada simulación.

No obstante, en la práctica, siendo \(\mu\) desconocido, la anterior sutileza carece prácticamente de relevancia. Además, si yo pensase —y no hiciese público— un valor de \(\mu\), obtuviese a partir de él una muestra y calculase un intervalo de confianza al 95%, ¿hasta cuánto estarías dispuesto a apostar por ganar un euro si el intervalo efectivamente contiene el valor que he pensado previamente?

Si no se conoce \(\sigma\), que es lo más habitual, se puede recurrir al principio del plug-in y reemplazar en la fórmula anterior por su estimación, la desviación estándar de la muestra, i.e.,

\[s^2 = \frac{1}{n-1} \sum_i (x_i - \bar{X})^2\]

y entonces se sabe que la expresión

\[\frac{\bar{X} - \mu}{s / \sqrt{n}}\]

tiene una distribución t de Student con \(n-1\) grados de libertad. Por lo que ya sabemos que para valores de \(n\) no particularmente grandes, es normal de nuevo. Por eso los intervalos de confianza usados más arriba pueden aplicarse en este caso con dos salvedades:

  • El valor \(\sigma\) desconocido se reemplaza por su aproximación \(s\).
  • Se pueden usar los cuantiles de la distribución normal solo si \(n\) es grande; en otro caso habría que usar qt(c(0.025, 0.975), df = n-1). (Por fijar ideas, con 20 observaciones, el intervalo de confianza usando la distribución t de Student es apenas un 6% más ancho que el obtenido con la distribución normal; con 100 observaciones, un 1.2%).

Los dos ejemplos planteados más arriba son sumamente particulares. Pero haciendo abstracción, hemos considerado en ambos un estimador \(\hat{\theta}_n\) de un parámetro de interés \(\theta_0\) y hemos podido deducir la distribución de \(\hat{\theta}_n - \theta_0\), que resulta ser \(N(0, \sigma / \sqrt{n})\). Luego, a partir de esta distribución, hemos podido construir intervalos de confianza. En el resto de la sección vamos a extender estos resultados a otros casos más generales.

El intervalo \(\bar{X} \pm 1.96 s / \sqrt{n}\) recibe el nombre de intevalo (de confianza) estándar para la media. Así, en general, todos los intervalos de confianza construidos de esa manera —se verá a continuación que no es la única— reciben también el adjetivo de estándar.

6.3.2.3 Intervalos de confianza

Es conveniente hacer un resumen de las dos secciones anteriores. Por un lado, al estudiar el caso de la media de valores de una distribución normal, se ha visto cómo dicho estimador tiene distribución \(N(\mu, \sigma / \sqrt{n})\) y se ha llegado a la conclusión de que un intervalo de confianza al 95% para dicho estimador es \(\bar{X} \pm 1.96 s / \sqrt{n}\) donde, recuérdese, \(s / \sqrt{n}\) es una estimación de la desviación estándar de \(\bar{X}\).

En el caso general, se sabe que la distribución de un estimador razonable \(\hat{\theta}\) es aproximadamente \(N(\theta, \sigma / \sqrt{n})\), luego generalizando el razonamiento aplicado a la distribución normal, se puede concluir que un intervalo de confianza al 95% para \(\hat{\theta}\) estimador es \(\hat{\theta} \pm 1.96 \sigma_\hat{\theta}\), donde \(\sigma_\hat{\theta}\) es una estimación de la desviación estándar de \(\hat{\theta}\).

De todos modos, el intervalo \(\hat{\theta} \pm 1.96 \sigma_\hat{\theta}\) está construido a partir de una aproximación a la distribución de \(\hat{\theta}\), que se ha construido de dos maneras distintas:

  1. Teóricamente, según la cual es normal (y de donde procede dicho intervalo).
  2. De una manera ḿas explícita a través de la técnica del plugin o de los remuestreos.

Eso unifica las dos definiciones de intervalo de confianza empleadas a lo largo de del capítulo, es decir, aquella según la cual el intervalo de confianza contiene al valor verdadero del parámetro con cierto nivel de garantías por un lado y aquella según la cual es el rango de variabilidad esperado del estimador al repetir el experimento cierto número de veces. Que ambas definiciones coincidan no es obvio; lo hacen en tanto que es posible garantizar la normalidad de los errores y aplicar el razonamiento anterior, que a decir de Fisher (véase (Efron and Hastie 2016)) es obviamente correcto.

De hecho, como la normalidad de la distribución del error solo puede garantizarse para valores de \(n\) muy grandes, en ocasiones, el intervalo de confianza estándar, es decir, \(\hat{\theta} \pm 1.96 \sigma_\hat{\theta}\), es incorrecto (o muy mejorable). Existe teoría que queda fuera del alcance de este libro acerca de la construcción de intervalos de confianza más finos en esos casos.

No obstante, no son necesarios valores de \(n\) muy grandes siempre para que los intervalos de confianza basados en remuestreos (véase el capítulo 11.2 de (Efron and Hastie 2016)) proporcionen intervalos de confianza razonables. Existen otros argumentos distintos del planteado hoy aquí que lo justifican en otros casos. Por ejemplo, que es posible construir transformaciones ad hoc de los datos que los normalicen de manera que pueda volver a aplicarse el ranzonamiento obviamente correcto de Fisher.

Así que en general, es bastante seguro construir intervalos de confianza —y ahora sí, en el sentido de que contienen el valor vedadero del parámetro con un nivel de certeza dado— a partir de cuantiles de la distribución de \(\hat{\theta}\) obtenidos por remuestreos.

6.4 Funciones de pérdida

Como se ha visto previamente en este capítulo, el principio de la máxima verosimilitud conlleva maximizar una determinada función. También cómo en algunos casos, la aplicación de este principio lleva a la mimimización de expresiones del tipo \(\sum_i \left(y_i - f_\theta(x_i)\right)^2\). En ciencia de datos es frecuente estimar los parámetros de una distribución —donde esos parámetros son típicamente coeficientes de un modelo— recurriendo al argumento de la minimización de algún tipo de función de pérdida similar a la anterior. En esta sección se discute esta aproximación a la estimación puntual de parámetros.

Es perfectamente posible —sin entrar a valorar si eso es positivo o no— avanzar significativamente en la carrera de científico de datos sin haber oído hablar de la función de verosimilitud. Pero es totalmente imposible sin haber hecho algún tipo de estimación puntual. La ciencia de datos parece en ocasiones una mera colección de algoritmos y es la estadística la que nos guía y nos indica cuál cabe aplicar en un momento dado.

De nuevo, ilustraremos cómo estimar parámetros de una distribución mediante un ejemplo. Si \(P_\mu\) es una distribución con media \(\mu\), entonces podemos definir la función del parámetro \(\theta\) como

\[l(\theta) = \int (x - \theta)^2 dP_\mu(x).\]

Se puede probar que esta función tiene su mínimo en \(E(X) = \mu\).

La función se llama \(l\) se llama función de pérdida y su nombre viene del inglés, loss. La integral lo es sobre la expresión \((x - \theta)^2\) que tiene valores grandes cuando \(x\) está muy alejado de \(\theta\). Miminizar la pérdida significa encontrar un punto \(\theta\) razonablemente próximo a todos los \(x\).

Así que dada una muestra \(x_i\) de esa distribución, como

\[l(\theta) \approx \frac{1}{N} \sum_i (x_i - \theta)^2\]

y el mínimo de la expresión de la derecha —donde frecuentemente se obvia el factor \(1/N\)— es la media muestral, \(\frac{1}{n} \sum_i x_i\), podemos usarla como aproximación a \(\mu\).

La discusión anterior permite redefinir la media de una manera inhabitual: aquél valor de una distribución que miminiza el error cuando este se mide cuadráticamente.

En general, si existe una función \(l(x, \theta)\) —vamos a usar, abusando del lenguaje, el mismo nombre, \(l\), tanto para la pérdida puntual como para la global— tal que

\[l(\theta) = \int f(x, \theta) \; dP_\mu(x)\]

tiene un mínimo en el parámetro de interés, podemos obtener una aproximación a él minimizando la expresión

\[\sum_i l(x_i, \theta)\]

En el caso anterior, \(l(x, \theta) = (x - \theta)^2\) es la llamada pérdida cuadrática. Este tipo de funciones de pérdida, que se usa muy frecuentemente, está relacionado como se ha indicado más arriba, con la distribución normal. Pero que es solo una de las posibles. Como alternativa se puede usar \(l(x, \theta) = |x - \theta|\), que tiene su mínimo en la mediana de la distribución y que es menos sensible a outliers.

Estas dos funciones de pérdida son las más populares. Sin embargo, es posible usar otras distintas con propiedades deseables. Por ejemplo, que sean aún más robustas a outliers. Existe, de hecho, una subdisciplina entera de la estadística, la de los M-estimadores, que estudia estas funciones alternativas y sus propiedades. Los M-estimadores son el fundamento de los llamados estimadores (y modelos) robustos, que son altenativas a los tradicionales cuando las observaciones presentan patrones de ruido excesivo.

En el gráfico anterior se muestra el perfil de tres funciones de pérdida: la cuadrática, la del valor absoluto y la de Tukey. Se observa cómo el impacto de una observación muy alejada en el último caso tiene el mismo peso que una razonablemente alejada; no ocurre lo mismo con las otras dos. El cuadrado de la cuadrática, además, da un peso superior a estas observaciones que el que les asigna la pérdida basada en el valor absoluto.

En concreto, con pérdidas cuadráticas, un outlier puede tener una influencia enorme sobre la estimación. El uso de otro tipo de funciones de pérdida reduce el impacto de los outliers y por lo tanto, proporcionar estimaciones más robustas.

El uso de funciones de pérdida es muy habitual cuando se ajustan modelos, algo que ya se ha comentado arriba y a lo que se volverá más adelante, pero sobre que merece la pena abundar en este capítulo. Supongamos que dos variables aleatorias \(X\) e \(Y\) dependen la una de la otra de manera que, p.e., \(Y | X \sim N(a + b X, \sigma)\). Eso significa que \(Y\) depende de \(X\) y de parámetros \(a\) y \(b\): para cada valor de \(X\), \(Y\) tiene una distribución normal y su media (supuesto \(b > 0\)) es tanto mayor cuando mayor sea \(X\).

Dada una muestra \((x_i, y_i)\) de la distribución conjunta de \(X\) e \(Y\), podemos estimar \(a\) y \(b\) minimizando la función de pérdida

\[l(a,b) = \sum_i (y_i - a - b x_i)^2.\]

De hecho, ese procedimiento es muy habitual, prácticamente el procedimiento por defecto, en ciencia de datos. Pero hay que señalar que:

  1. Esa función de pérdida, \(l(a, b)\), es precisamente a la que se llegaría aplicando el principio de la maximización de la verosimilitud; la demostración es sencilla, una extensión de la realizada más arriba para probar cómo la media muestral es el estimador por máxima verosimilitud de la media de una distribución normal.
  2. Que el quee la périda cuadrática es consecuencia de considerar errores normales. En ciencia de datos es habitual recurrir a la pérdida cuadrática simplemente por tradición o facilidad de uso, obviando la relación que existe con la normalidad en la distribución de los errores.

El uso de funciones de pérdida tiene una relación muy estrecha con la toma de decisiones. Las funciones de pérdida consideradas más arriba (p.e., la cuadrática, \((x - \theta)^2\)) son tan comunes como arbitrarias. En situaciones reales pueden existir funciones de pérdida propiamente dichas que cuantifiquen el error cometido en, por ejemplo, euros. En tales circunstancias, es perfectamente razonable utilizarlas directamente en lugar de cualquiera de las de libro. Por ejemplo, en un hospital pueden querer estimar la capacidad reservada para el servicio de urgencias. En tal caso, una sobreestimación (que dilapida recursos), aunque indeseable, es mucho menos grave que una infraestimación (que dejaría a pacientes con una antención inadecuada). En particular, la pérdida no es simétrica y la función de pérdida debería penalizar más severamente el error por defecto que el por exceso. Ese tipo de funciones de pérdida puede llevar a elegir como parámetro no ya una medida de centralidad de la distribución sino un cuantil de la misma.

El siguiente ejercicio numérico desarrolla el ejemplo anterior. Se van a simular datos de pacientes que llegan a un servicio de urgencias usando una distribución de Poisson (obviando cuestiones como la tendencia y la estacionalidad) y se van a asignar unos costes fijos \(a * \theta\), donde \(\theta\) es la capacidad y \(a\) es el coste fijo por plaza, y otro variable \(b * (n - \theta)\), cuando el número de pacientes que lleguen en un día concreto, \(n\), exceda la capacidad reservada.

Se puede calcula el coste anual asociado a cada decisión \(\theta\), así:

El coste anual según el valor \(\theta\) elegido es

Que tiene un valor mínimo en 39 —aunque puede variar según la simulación de los datos—. En este caso, se ha estimado un parámetro de la distribución de Poisson que no es su intensidad, \(\lambda\), para lo que habría que haber utilizado la función de péridida cuadrática, sino otro, \(\theta\), asociado a nuestra particular y circunstancial definción de la función de pérdida \(l\).

Es posible —y queda propuesto como ejercicio al final del capítulo— estudiar la distribución del coste anual asociado a cada valor de \(\theta\), así como la posible variabilidad de la misma \(\theta\).

Este ejemplo ilustra la estrecha vinculación entre el análisis de los datos y el proceso de toma de decisiones. Desafortunadamente, en demasiadas situaciones existe una importante brecha de comunicación entre quienes analizan los datos y los responsables de tomar las últimas decisiones. Este es el motivo, de hecho, de muchas de las decisiones incorrectas que toman instituciones púbicas y privadas.

6.5 Bibliografía razonada

Existen objeciones a la atribución a Pearson (Pearson 1894) del método de los momentos. Parece haber un precedente en Chebyshev (Encyclopedia of Mathematics 2021), aunque puede argumentarse que Pearson usó el procedimiento de la manera más próxima a cómo se emplea habitualmente en estadística.

El material de esta sección está desarrollado con mayor profundidad y rigor en otras fuentes, como (Efron and Hastie 2016) o (Sköld 2006). En particular, (Efron and Hastie 2016) discute aspectos sobre la creación de intervalos de confianza cuando las aproximaciones asintóticas no tienen justificación.

En este capítulo se ha dado importancia a la estimación por máxima verosimilitud y a las pérdidas cuadráticas que se deducen de ella. Entre otras cosas, porque producen estimadores insesgados de varianza mínima (cota de Cramér-Rao); sin embargo, es posible conseguir estimadores sesgados con una varianza menor. Algunos de ellos, como los que estudia la regresión ridge (Hastie 2020) son muy habituales en ciencia de datos y están basados en una modificación ad hoc (aunque justificada) de la función de pérdida cuadrática habitual.

6.6 Ejercicios

Ejercicio 6.1 En el ejemplo del hospital, crear un intervalo de confianza para el coste anual asociado al valor \(\theta = 39\) mediante remuestreos.
Ejercicio 6.2 En el ejemplo del hospital, crear un intervalo de confianza \(\theta\) mediante remuestreos.

Referencias

Efron, B., and T. Hastie. 2016. Computer Age Statistical Inference: Algorithms, Evidence, and Data Science. https://doi.org/10.1017/CBO9781316576533.

Encyclopedia of Mathematics. 2021. “Moments, Method of (in Probability Theory).” http://encyclopediaofmath.org/index.php?title=Moments,_method_of_(in_probability_theory)&oldid=47882.

Hastie, Trevor. 2020. “Ridge Regularizaton: An Essential Concept in Data Science.” http://arxiv.org/abs/2006.00371.

Pearson, Karl. 1894. “Contributions to the Mathematical Theory of Evolution.” Philosophical Transactions of the Royal Society of London 3. https://doi.org/10.1098/rsta.1894.0003.

Perry, Patrick. 2015. “Fast Moment-Based Estimation for Hierarchical Models.” Journal of the Royal Statistical Society: Series B (Statistical Methodology), April. https://doi.org/10.1111/rssb.12165.

Sköld, M. 2006. Computer Intensive Statistical Methods. https://www.math.kth.se/matstat/gru/sf2955/Kompendium/kompendium_sf2955.pdf.