import numpy as np
import holoviews as hv
hv.extension('bokeh')
from bokeh.plotting import show

8. Modelamiento Paramétrico con MLE#

En esta lección, empezamos a ver la tarea de la estimación de parámetros.

En el caso de inferencia paramétrica asumimos que las observaciones siguen una cierta distribución, i.e. las observaciones son una realización (muestra) de un proceso aleatorio.

Los pasos conceptuales de la inferencia paramétrica para la estimación de parámetros son

  1. Proposición: Se propone un modelo paramétrico en base a un proceso exploratorio sobre los datos

  2. Ajuste: Se encuentran los valores de los parámetros que mejor ajustan (fit) los datos

  3. Evaluación (y modificación): Se estudian los errores del modelo y de ser necesario se proponen modificaciones que acomodan características importantes de los datos que hayan sido omitidas anteriormente. Luego se vuelve al punto 2.

En el paradigma frecuentista el paso 2 se resuelve típicamente utilizando la estimación de máxima verosimilitud (Maximum Likelihood Estimation, MLE). Otro método frecuentista es el método de los momentos. En esta lección nos concentraremos exclusivamente en el MLE.

8.1. La verosimilitud#

En términos simples, la verosimilitud expresa qué tan posible (relativamente) que el conjunto de datos observados haya sido generado por un parámetro (o parámetros).

Nos dice que tan bueno es nuestro modelo (con su parámetro) con respecto a los datos observados.

Veamos ahora una definición matemática

  • Sea el resultado de un experimento el cual modelamos como un conjunto de variables aleatorias (VAs) (discretas o continuas): \(X_1, X_2, \ldots, X_N\)

  • Para cada una tenemos también observaciones (realizaciones): \(\{x_i\}_{i=1...N} = x_1, x_2, \ldots, x_N\)

  • Asumiremos que las VAs siguen una cierta distribución conjunta: \(f(x_1, x_2, \ldots, x_N; \theta)\) con el parámetro \(\theta\) (o el vector de parámetro)

En base a esto la función de verosimilitud se define como

\[ \begin{split} \mathcal{L}(\theta) &= f(x_1, x_2, \ldots, x_N ; \theta) \end{split} \]

¿Por qué usamos “;” en lugar de “|” en la función de densidad?

Dado que nuestra meta es estimar el parámetro, ¿\(\mathcal{L}(\theta)\) es una función de \(\{x_i\}_{i=1...N}\) o del parámetro?

\(\mathcal{L}(\theta)\) es una función del parámetro \(\theta\) (o el vector de parámetro).

¿Cuál es el rango de su valor?

\(\mathcal{L}(\theta)\) tiene valor entre [0, \(\infty\)).

La verosimilitud puede decirnos que tan posible (relativamente) es que \(\{x_i\}_{i=1...N}\) haya sido generado por la distribución \(f\) con parámetros \(\theta\).

En particular, mientras mayor es la verosimilitud, mejor es nuestro modelo, i.e. la maximización de la verosimilitud nos permite encontrar el mejor \(\theta\) para nuestras observaciones.

Advertencia

  • La verosimilitud no es una distribución de probabilidad de \(\theta\), y no integra a uno respecto a \(\theta\). La verosimilitud no puede interpretarse como la probabilidad de \(\theta\)

  • El valor de la verosimilitud en si mismo no tiene mucho significado, pero puede usarse para hacer comparaciones entre distintos valores del parámetro \(\theta\)

Para esta lección, asumiremos que nuestras observaciones son independientes e idénticamente distribuidas (i.i.d.). ¿Cómo podemos simplificar la expresión \(f(x_1, x_2, \ldots, x_N ; \theta)\)?

\[\begin{split} \begin{split} \mathcal{L}(\theta) &= f(x_1; \theta) \cdot f(x_2; \theta) \cdot \ldots \cdot f(x_N; \theta) \\ &= \prod_{i=1}^N f(x_i; \theta) \end{split} \end{split}\]

8.2. Estimación de máxima verosimilitud (MLE)#

En modelamiento paramétrico estamos interesados en encontrar el valor de \(\theta\) que mejor ajusta a nuestras observaciones. El método que utilizaremos es MLE:

  • 1 Seleccionar una distribución (modelo) para las observaciones y formular una verosimilitud \(\mathcal{L}(\theta)\)

  • 2 Buscar el valor de \(\theta\) que maximiza \(\mathcal{L}(\theta)\) dado los datos, i.e.

\[ \hat \theta = \text{arg} \max_\theta \mathcal{L}(\theta), \]

donde el estimador puntual \(\hat \theta\) se llama estimador de máxima verosimilitud de \(\theta\) (en esta lección, también usamos MLE para referir este tipo de estimador).

Luego de esto podemos

  • 3 Determinar el intervalo de confianza de \(\hat \theta\) ya sea de forma analítica o numérica (bootstrap, etc)

  • 4 Evaluar el modelo y modificarlo si es necesario.

Importante

Un supuesto equivocado en el paso 1 invalidará nuestra inferencia. ¿Cómo podemos proponer/encontrar un buen modelo?

  • Proponemos un modelo (o una distribución) en base a un proceso exploratorio sobre los datos

  • Evalua que el método sea apropiado con algunas métricas o gráficos

  • Compara con otros modelos

  • Sugiere mejoras incrementales

8.2.1. Log-verosimilitud#

En muchos casos es más práctico utilizar logaritmo de la verosimilitud (log-verosimilitud).

¿Cómo se escribe log-verosimilitud y cómo simplificarla en el caso de i.i.d.?

\[\begin{split} \begin{split} log(\mathcal{L}(\theta)) &= log(f(x_1, x_2, \ldots, x_N ; \theta))\\ &= log(f(x_1; \theta) \cdot f(x_2; \theta) \cdot \ldots \cdot f(x_N; \theta)) \\ &= log(\prod_{i=1}^N f(x_i; \theta)) \\ &= \sum_{i=1}^N log f(x_i; \theta) \end{split} \end{split}\]

En contexto de verosimilitud, \(log\) normalmente se refiere a \(log_e\) o \(ln\).

¿Por qué usamos log-verosimilitud en lugar de verosimilitud en práctica?

  • El logaritmo es una función monotónica.

  • Es más fácil manejar la suma que la multiplicación.

  • Mejor estabilidad numérica: al realizar una transformación logarítmica se pueden convertir decimales pequeños en valores más grandes que un ordenador con precisión finita puede manejar mejor.

  • El logaritmo es particularmente conveniente para las distribuciones de familia exponencial.

¿Cuál es el rango de la log-verosimilitud?

8.2.2. Ejemplo: MLE para la media de una distribución Gaussiana#

Consideremos un conjunto de N mediciones \(\{x_i\}_{i=1,\ldots, N}\) obtenidas de un sensor y que todas fueron obtenidas bajo las mismas condiciones y son independientes (i.i.d).

Basado en un análisis exploratorio, suponemos que el sensor tiene un error que sigue una distribución Gaussiana con varianza conocida \(\sigma^2\).

¿Cómo podemos modelar los datos medidos por el sensor?

Las mediciones pueden verse como muestras ruidosas en torno a un valor verdadero (intrínseco) que denominaremos \(\mu\):

\[ x_i = \mu + \epsilon_i, \quad \epsilon_i \sim \mathcal{N}(0, \sigma^2), \]

y por lo tanto

\[ x_i \sim \mathcal{N}(\mu,\sigma^2) \]
\[ f(x_i) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp \left( -\frac{(x_i-\mu)^2}{2\sigma^2}\right) \]

¿Cuál es la verosimilitud de \(\mu\) dadas las mediciones \(\{x_i\}_{i=1...N}\) y la varianza \(\sigma^2\)?

\[ \mathcal{L}(\mu) = f(\{x_i\}_{i=1...N}; \mu, \sigma^2) = \prod_{i=1}^N f(x_i; \mu, \sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}^N} \prod_{i=1}^N \exp \left( -\frac{(x_i-\mu)^2}{2\sigma^2} \right) \]

Podemos encontrar el valor de \(\mu\) maximizando la expresión anterior. Pero es mejor maximizar la log-verosimilitud por las razones expuestas anteriormente.

¿Cuál es la log-verosimilitud?

\[\begin{split} \begin{split} \log \mathcal{L} (\mu) &= \log \prod_{i=1}^N f(x_i ; \mu, \sigma^2) \\ &= \sum_{i=1}^N \log f(x_i; \mu, \sigma^2) \\ &= - \frac{1}{2} \sum_{i=1}^N \log 2\pi\sigma^2 - \frac{1}{2} \sum_{i=1}^N \frac{(x_i-\mu)^2}{\sigma^2} \\ &= - \frac{N}{2} \log 2\pi\sigma^2 - \frac{1}{2\sigma^{2}} \sum_{i=1}^N (x_i-\mu)^2 \end{split} \end{split}\]

¿Cómo encontrar el valor de \(\mu\) dado la log-verosimilitud?

  • Podemos maximizar la log-verosimilitud tomando su derivada e igualando a cero. Por qué?

Puedes derivar el estimador?

\[ \frac{d \log \mathcal{L} (\mu)}{d\mu} = \frac{1}{\sigma^{2}} \sum_{i=1}^N (x_i-\mu) =0 \]
\[ \hat \mu = \frac{1}{N} \sum_{i=1}^N x_i, \quad \sigma >0, \]

Esto corresponde al estimador de máxima verosimilitud de \(\mu\).

Nota

  • La media muestral es el MLE de la media para una distribución Gaussiana. También es el MLE de la media para varias distribuciones.

  • En la práctica, las librerías suelen minimizar la log-verosimilitud negativa (como función de pérdida/coste).

8.2.3. Ejemplo: MLE para la varianza de una distribución Gaussiana#

Supongamos ahora que no se conoce la varianza del sensor en el ejemplo anterior. Cuál es el MLE de la varianza?

\[ \log \mathcal{L} (\mu, \sigma^2) = - \frac{N}{2} \log 2\pi\sigma^2 - \frac{1}{2\sigma^{2}} \sum_{i=1}^N (x_i-\mu)^2 \]

El MLE de la varianza puede obtenerse utilizando el mismo procedimiento. En el proceso abajo, podemos crear una nueva variable \(\beta=\sigma^2\) para calcular la derivada, y al final reemplazamos \(\beta\) por \(\sigma^2\). Podría ser más fácil llegar al término \(\sigma^{4}\) de esta manera.

\[ \frac{d \log \mathcal{L} (\mu, \sigma^2)}{d\sigma^2} = - \frac{N}{2} \frac{1}{\sigma^2} + \frac{1}{2\sigma^{4}}\sum_{i=1}^N (x_i-\mu)^2 =0 \]
\[ \hat \sigma^2 = \frac{1}{N} \sum_{i=1}^N (x_i- \mu)^2 \]

Generalmente, \(\mu\) no se conoce, entonces se usa media muestral \(\hat \mu\):

\[ \hat \sigma^2 = \frac{1}{N} \sum_{i=1}^N (x_i- \hat \mu)^2 \]

8.2.4. Ejemplo: MLE para una mezcla de Gaussianas#

Imaginemos ahora que nuestros datos iid provienen de una mezcla finita de \(K\) distribuciones Gaussianas:

\[ f(x_i; \pi,\mu,\sigma^2) = \sum_{k=1}^K \pi_k \mathcal{N}(x_i; \mu_k, \sigma_k^2), \]

donde \(\sum_{k=1}^K \pi_k = 1\) and \(\pi_k \in [0, 1] ~~ \forall k\)

La verosimilitud es

\[ \log \mathcal{L}(\pi,\mu,\sigma^2) = \sum_{i=1}^N \log \sum_{k=1}^K \pi_k \mathcal{N}(x_i; \mu_k, \sigma_k^2) \]

Error

No se puede obtener una expresión analítica de los parámetros. En casos como este se utilizan métodos aproximados y/o iterativas, e.g. gradiente descendente o expectation maximization. Veremos más sobre esto en el futuro.

8.3. Propiedades de los estimadores#

Ahora, veamos las propieades de un estimador!

En el ejemplo del sensor, ya hemos derivado las expresiones de estimadors MLE para la media y varianza de la districución Gaussiana. Veamos cómo estos MLEs cambiar a medida que se observan más datos!

# data from the sensor
np.random.seed(1234)
x = 80 + np.random.randn(10000)
#x = 80 + 2*np.random.rand(1000)  # What happens if the data is not normal

# Computing the MLE
Ns = np.round(np.logspace(0, 4, num=16)).astype('int')
hat_mu = np.array([np.mean(x[:N]) for N in Ns])
hat_s2 = np.array([np.mean((x[:N]-mu)**2) for mu, N in zip(hat_mu, Ns)])

mu_plot = hv.Curve((Ns, hat_mu), 'Number of samples', 'hat mu')*hv.HLine(80).opts(color='k')
s2_plot = hv.Curve((Ns, hat_s2), 'Number of samples', 'hat sigma2')*hv.HLine(1).opts(color='k')
plot = (mu_plot + s2_plot).opts(hv.opts.Curve(logx=True, width=350), hv.opts.HLine(line_dash='dashed'))

show(hv.render(plot))

Lo que hemos observado es sobre la consistencia de un estimador.

8.3.1. Consistencia#

Para un parámetro \(\theta\) y un estimador \(\hat \theta\), si

\[\lim_{N \to \infty} P\left(| \hat \theta - \theta| < \epsilon \right) = 1 \]

para cualquier \(\epsilon >0\) (\(N\): tamaño de la muestra). O se escribe como

\[\begin{split}\begin{array}{c} \hat \theta \overset{p}{\to} \theta \\ (\hat \theta \text{ converge en probabilidad a } \theta \text{ cuando } N \to \infty) \end{array}\end{split}\]

entonces \(\hat \theta\) es un estimador consistente. Es decir, un estimador consistente aproxima el valor del parámetro cuanto mayor es \(N\) (tamaño de la muestra).

Hay otras propiedades para evaluar un estimador!

8.3.2. Insesgadez#

Para un parámetro \(\theta\) y un estimador \(\hat \theta\), si

\[E[\hat \theta] = \theta\]

entonces \(\hat \theta\) es un estimador insesgado de \(\theta\). Es decir, la distribución de un estimador insesgado se centra en el parámetro verdadero.

Nota

La distribución de un estimador es una distribución de muestreo. Un estimador es un estadístico cuyo valor se puede determinar a partir de los datos muestrales. Cuando un estadístico se usa para estimar un parámetro de la población, el estadístico se llama un estimator.

Consistencia vs. Insesgadez

  • “Consistente” es diferente que “insesgado”. Un estimador puede ser consistente pero sesgado (ej el MLE estimador de varianza arriba), o inconsistente pero insesgado (ej \(x_1\) de \(\{x_i\}\) para estimar la media).

  • Ambos se refieren a la distribución del estimador, pero “consistente” se define considerando cuando el tamaño de la muestra se vuelve grande, y “sesgado” se define considerando la esperanza del estimador dado cualquier tamaño de la muestra.

  • La propiedad insesgado solía recibir mucha atención, pero en estos días se considera menos importante; muchos de los estimadores que usaremos son sesgados. Un requisito razonable para un estimador es que debe converger al valor del parámetro verdadero a medida que recopilamos más y más datos, es decir, debe ser consistente.

8.3.3. Eficiencia#

Un estimador es eficiente cuando posee varianza mínima, o alcanza la cota inferior de Cramér-Rao.

Aquí consideramos la cota inferior de Cramér-Rao para los estimadores insesgados. El teorema de la cota de Cramér-Rao determina que la varianza de un estimador insesgado \(\hat\theta\) de un parámetro \(\theta\) tiene una cota inferior como el inverso de la información de Fisher:

\[ Var(\hat\theta) = \sigma_{\hat\theta}^2 \geq \frac{1}{{\cal I}_N(\theta)} \]

  • \({\cal I}_N(\theta)\) es la información de Fisher de una muestra con el tamaño \(N\), y se puede calcular basado en la negativa esperanza de la segunda derivada de la log-verosimilitud de cada observación \(x\) (o \(x_i\)) con respecto a \(\theta\) (bajo ciertas condiciones):

\[ {\cal I}_N(\theta) = - N \cdot E_{X;\theta}\left[ \frac{d^2 \log \mathcal{L} (\theta)}{d\theta^2} \right] = - N \int \frac{d^2 \log f(x;\theta)}{d\theta^2} f(x;\theta) dx \]

  • \(E_{X;\theta}\) denota la esperanza con respecto a la densidad \(f(x; \theta)\). Si no se indica, en lo que sigue, la esperanza se toma con respecto a \(X\).

  • La información de Fisher puede verse como la curvatura del gráfico de log-verosimilitud.

  • En algunos casos, no existe un estimador insesgado que alcance la cota inferio.

  • En algunos casos, un estimador sesgado puede tener menor varianza que la cota inferior de Cramér-Rao definida arriba definida para los estimadores insesgados.

Para un estimador sesgado, se aplica otra formula más general de la cota inferior de Cramér-Rao: General scalar case.

8.3.4. Normalidad asintótica#

Para un parámetro \(\theta\) y un estimador \(\hat \theta\), si

\[\begin{split}\begin{array}{c} \hat \theta \overset{d}{\to} \cal{N}(\theta, \frac{\sigma_\theta^2}{\textit{N}})\\ \\ (\hat \theta \text{ converge en distribución a una distribución Normal } \cal{N}(\theta, \frac{\sigma_\theta^2}{\textit{N}}) \text{ cuando } \textit{N} \to \infty) \end{array}\end{split}\]

o equivalente

\[\begin{array}{c} \sqrt{n} (\hat \theta - \theta) \overset{d}{\to} \cal{N}(0, \sigma_\theta^2) \end{array}\]

entonces \(\hat \theta\) es asintóticamente normal.

Ahora, veamos si los MLEs alcanzan todas estas propiedades!

8.4. La propiedad sobre insesgadez de los MLEs#

Veamos la propiedad sobre insesgadez de los MLEs primero. Continuamos el ejemplo de la distribución Gaussiana arriba. Preguntamos:

¿Es el MLE de \(\mu\) para la distribución Gaussiana insesgado?

\[\begin{split} \begin{split} E[\hat \mu] &= E \left[ \frac{1}{N} \sum_{i=1}^N x_i \right] \\ &= \frac{1}{N} \sum_{i=1}^N E[x_i] = \frac{1}{N} \sum_{i=1}^N \mu = \mu \end{split} \end{split}\]

La respuesta es SI.

¿Es el MLE de \(\sigma^2\) para la distribución Gaussiana insesgado?

Tenemos varias maneras de derivaciones. Aquí les muestra una manera. Puedes dar una derivación más simple?

\[\begin{split} \begin{split} \hat \sigma^2 &= \frac{1}{N} \sum_{i=1}^N \left(x_i- \frac{1}{N}\sum_{j=1}^N x_j \right)^2 \\ &= \frac{1}{N} \sum_{i=1}^N x_i^2 - \frac{1}{N^2} \sum_{i=1}^N \sum_{j=1}^N x_i x_j \\ &= \frac{1}{N} \sum_{i=1}^N x_i^2 - \frac{1}{N^2} \sum_{i=1}^N \sum_{j\neq i} x_i x_j - \frac{1}{N^2} \sum_{i=1}^N x_i^2 \\ &= \frac{N-1}{N^2} \sum_{i=1}^N x_i^2 - \frac{1}{N^2} \sum_{i=1}^N \sum_{j \neq i} x_i x_j \end{split} \end{split}\]

Luego aplicando el valor esperado tenemos

\[\begin{split} \begin{split} E[\hat \sigma^2] &= \frac{N-1}{N^2} E[\sum_{i=1}^N x_i^2] - \frac{1}{N^2} E[\sum_{i=1}^N \sum_{j \neq i} x_i x_j] \\ &= \frac{N-1}{N^2} \sum_{i=1}^N E [x_i^2] - \frac{1}{N^2} \sum_{i=1}^N \sum_{j \neq i} E [x_i] E [x_j] \text{ (Por qué?)} \\ &= \frac{N-1}{N} (\sigma^2 + \mu^2) - \frac{N-1}{N} \mu^2 \\ &= \frac{N-1}{N} \sigma^2 \neq \sigma^2 \end{split} \end{split}\]

La respuesta es NO.

¿Es posible corregir este sesgo?

Si. Si multiplicamos por \(\frac{N}{N-1}\) se obtiene el conocido estimador insesgado de la varianza

\[ \hat \sigma_{u}^2 = \frac{N}{N-1} \hat \sigma^2 = \frac{1}{N-1} \sum_{i=1}^N (x_i- \hat\mu)^2 \]

El estimador \(\hat \sigma^2 = \frac{1}{N} \sum_{i=1}^N \limits (x_i- \bar{x})^2\) también se llama varianza muestral sesgado (biased sample varianze).

El estimador \(\hat \sigma_{u}^2 = \frac{1}{N-1} \sum_{i=1}^N \limits (x_i- \bar{x})^2 = S^2\) también se llama varianza muestral insesgado (unbiased sample varianze), o varianza muestral.

Nota

En general, los MLEs no son estimadores insesgados. Pero, son asintóticamente insesgados, según su distribución asintótica que presentaremos próximamente.

8.5. Propiedadeds de los MLEs#

MLE tiene otras propiedades! Asumiendo que los datos vienen de un modelo especificado (Gaussiana o no), entonces el MLE usualmente tiene los siguientes propiedades:

  • Asintóticamente insesgado

  • Consistente

  • Asintóticamente eficiente, y alcanca la minima varianza asintótica de cualquier estimador insesgado

  • Asintóticamente Normal

En particular, veamos la última propiedad, asintóticamente Normal, que proviene del teorema de la distribución asintótica del MLE:

Para un parámetro \(\theta\) y un estimador MLE \(\hat \theta\) basado en una muestra de N observaciones i.i.d.,

\[\begin{split}\begin{array}{c} \hat \theta \overset{d}{\to} \cal{N}(\theta, \frac{1}{\cal{I}_\textit{N}(\theta)})\\ \\ (\hat \theta \text{ converge en distribución a una distribución Normal } \cal{N}(\theta, \frac{1}{\cal{I}_\textit{N}(\theta)})\text{ cuando } \textit{N} \to \infty) \end{array}\end{split}\]

La prueba se puede encontra aquí (y se aplica el Teorema del Límite Central en algún paso).

Nota

  • El teorema arriba ayuda a deducir las propiedades de asintóticamente insesgado y asintóticamente eficiente.

  • Las propiedades arriba solo se mantienen bajo ciertas condiciones regulares. Estas son esencialmente condiciones de suavidad en \(f(x;\theta)\). Usualmente en las situaciones que encuentres, estas condiciones se mantienen.

Truco

Según las propiedades, el MLE es un buen estimador, por lo menos en muestra grande. Aunque no sabemos qué tan grande la muestra debe ser, en general es bueno usar MLE si no sabemos mejores estimadores.

8.5.1. Ejemplo: Desviación estándar del MLE de \(\mu\) para una distribución Gaussiana#

Cuando hacemos la estimación de parámetros desde una muestra, es importante saber el incertidumbre o presición de un estimador. Para cuantificar la incertidumbre o presición de un estimador, ¿qué podemos calcular?

Desviación estándar del estimador!

Aquí, mostramos un ejemplo dell estimador MLE \(\hat \mu\) de \(\mu\) para una distribución Gaussiana del ejemplo anterior. Tenemos dos pasos.

(1) En el primer paso

Calculamos la varianza del estimador. Podemos obtener la varianza de a través de dos maneras!

En la primera manera, ¿cuál es el valor de \(Var(\hat \mu)\) (o \(\sigma^2_{\hat \mu}\)), dado que \(\hat \mu = \bar x = \frac{1}{N}\sum_{i=1}^N \limits x_i\) (la media muestral)?

Para la varianza de \(\hat \mu = \frac{1}{N} \sum_{i=1}^N \limits x_i\),

\[\begin{split}\begin{split} Var(\hat \mu) = Var(\frac{1}{N}\sum_{i=1}^N \limits x_i) &= \frac{1}{N^2} ( \sum_{i=1}^N \limits Var(x_i)) \text{ (Por qué?)} \\ &= \frac{1}{N^2} ( N \sigma^2 ) \\ &= \frac{\sigma^2}{N} \end{split}\end{split}\]

Nota

Esta derivación se mantiene en general para \(\{x_i\}_{i=1...N}\) i.i.d., incluso si \(X_i\) no sigue la distribución normal (Gaussiana).

Por otro lado, podemos calcular la cota inferior de Crámer-Rao. Ya sabemos que \(\hat \mu\) es un estimador insesgado, y la verosimilitud es:

\[ \log \mathcal{L} (\mu) = - \frac{N}{2} \log 2\pi\sigma^2 - \frac{1}{2\sigma^{2}} \sum_{i=1}^N (x_i-\mu)^2 \]

La cota inferior de Crámer-Rao (o la varianza mínima) de los estimadores insesgados es:

\[\begin{split} \begin{split} \frac{1}{\cal{I}_N(\mu)} &= \frac{1}{- E\left [ \frac{d^2 \log \mathcal{L}(\mu)}{d\mu^2}\right]} \\ &= \frac{1}{- E \left [ \frac{d}{d\mu} \left( \frac{1}{\sigma^2}\sum_{i=1}^N \limits (x_i-\mu) \right)\right]} \\ &= \frac{1}{E\left [ \frac{N}{\sigma^2} \right]}\\ &= \frac{\sigma^2}{N} \end{split} \end{split}\]

Se puede ver que en este caso, la varianza del estimador MLE \(\hat \mu\) alcanza la cota inferior de Crámer-Rao (o la varianza mínima) de los estimadores insesgados en cualquier tamaño de muestra dado. Es decir, el estimador MLE \(\hat \mu\) para una distribución Gaussiana es eficiente.

Podemos obtener la varianza de \(\hat \mu\) a través de dos maneras muestradas arriba!

(2) En el segundo paso

Calculamos la raíz cuadrada de la varianze del estimador, que es el error estándar, la desviación estándar de la distribución muestral de un estadístico muestral o estimador en nuestro contexto.

Cuando el estimador (\(\hat \mu\)) es de la media, este error estándar (que es la desviación estándar de la media muestral) también se llama el error estándar de la media \(se(\hat \mu)\). Lo usamos para cuantificar la incertidumbre de \(\hat \mu\) (en el futuro, lo usaremos para estimar el intervalo de confianza).

En general, supongamos que una muestra de N observaciones i.i.d. se toma de una distribución (Gaussiana o no) con una desviación estándar \(\sigma\) o varianza \(\sigma^2\), la media muestral \(\hat \mu\) tendrá un error estándar de la media dado por:

\[ se(\hat \mu) = \sqrt{Var(\hat \mu)} = \sqrt{\frac{\sigma^2}{N}} = \frac{\sigma}{\sqrt{N}} \]

Ahora, veamos como cambia el MLE \(\hat \mu\) y su incertidumbre \(se(\hat \mu)\) para una distribución Gaussiana a medida que observamos más y más datos:

# Generate data
np.random.seed(12345)
mu_real, s_real = 2.23142, 1.124123
x = mu_real + s_real*np.random.randn(10000)
# MLE and its standard error
hat_mu = np.array([np.mean(x[:n]) for n in range(1, len(x))])
standard_error = s_real/np.sqrt(np.arange(1, len(x)))

mu_plot = hv.Curve((range(1, len(x)), hat_mu), 
                   kdims='Number of samples', vdims='mu', label='Estimated mu')
se_plot = hv.Spread((range(1, len(x)), hat_mu, standard_error), 
                    kdims='Number of samples', label='Standard error')
mu_real_plot = hv.HLine(mu_real).opts(line_dash='dashed', color='k', alpha=0.5)
plot = (mu_plot*se_plot*mu_real_plot).opts(hv.opts.Curve(logx=True, width=500))

show(hv.render(plot))