Contenido

Ajuste de modelos paramétricos con MLE

import holoviews as hv
hv.extension('bokeh')

import numpy as np
import scipy.signal

7. Ajuste de modelos paramétricos con MLE

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 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 los datos

  3. 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). Otros métodos frecuentistas son el método de los momentos y el estimador-M. En esta lección nos concentraremos exclusivamente en el MLE

7.1. La verosimilitud

La función de verosimilitud es una descripción cuantitativa de nuestro experimento o proceso de medición. Además es el punto de partida del modelamiento paramétrico en ambos paradigmas (F y B). En términos simples:

la verosimilitud nos dice que tan bueno es nuestro modelo 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 (VA): \(X_1, X_2, \ldots, X_N\)

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

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

En base a esto la verosimilitud se define como

\[\begin{split} \begin{split} \mathcal{L}(\theta) &= P(X_1=x_1, X_2=x_2, \ldots, X_N=x_n) \\ &= f(x_1, x_2, \ldots, x_N | \theta), \end{split} \end{split}\]

que es una función de los parámetros \(\theta\)

Para los ejemplos de esta clase asumiremos adicionalmente que nuestras observaciones son independientes e idénticamente distribuidas (i.i.d.). Con lo anterior podemos simplificar la expresión como

\[\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}\]

Nota

El valor de la verosimilitud en si mismo no tiene mucho significado, pero puede usarse para hacer comparaciones entre distintos valores del vector de 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

Verosimilitud no es lo mismo que probabilidad

  • La verosimilitud de una VA no integra a uno, i.e. en general la verosimilitud no es una densidad de probabilidad válida

  • La verosimilitud en su misma no puede interpretarse como la probabilidad de \(\theta\)

  • La verosimilitud sólo puede decirnos que tan posible es que \(\{x_i\}\) haya sido generado por la distribución \(f\) con parámetros \(\theta\)

7.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\)

Luego de esto podemos

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

  • 4 Extraer conclusiones sobre nuestro modelo (test de hipótesis)

Importante

Un supuesto equivocado en el paso 1 invalidará nuestra inferencia. Siempre evalua que el método sea apropiado, compara con otros modelos y sugiere mejoras incrementales

7.2.1. 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. El sensor que se utilizó para medir los datos tiene un error que sigue una distribución Gaussiana con varianza conocida \(\sigma^2\).

Si las condiciones son las mismas entonces 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 podemos escribir la distribución de \(x_i\) como

\[ f(x_i) = \mathcal{N}(x_i |\mu,\sigma^2) \quad \forall i. \]

Luego la verosimilitud de \(\mu\) dadas las mediciones y la varianza es

\[ \mathcal{L}(\mu) = f(\{x_i\}| \mu, \sigma^2) = \prod_{i=1}^N f(x_i| \mu, \sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}} \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

Truco

En muchos casos es más práctico encontrar el máximo del logaritmo de la verosimilitud. El logaritmo es una función monotónica y su máximo es idéntico al de su argumento

Lo anterior es particularmente conveniente en las distribuciones de familia exponencial. En el caso anterior la log-verosimilitud es

\[\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}\]

Que podemos maximizar tomando su derivada e igualando a cero:

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

Despejando se tiene que:

\[ \hat \mu = \frac{1}{N} \sum_{i=1}^N x_i, \quad \sigma >0, \]

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

Nota

La media muestral es el MLE de la media para una verosimilitud Gaussiana

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

Supongamos ahora que no se conoce la varianza del sensor en el ejemplo anterior. El MLE de la varianza puede obtenerse utilizando el mismo procedimiento:

\[ \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 \]
\[ \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- \hat\mu)^2 \]

Nota

Si el valor real de \(\mu\) no se conoce entonces este es un estimador sesgado de la varianza real. En general, el procedimiento MLE no garantiza que los estimadores sean insesgados

El siguiente ejemplo muestra la evolución del MLE de estos parámetros 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 s2')*hv.HLine(1).opts(color='k')
(mu_plot + s2_plot).opts(hv.opts.Curve(logx=True, width=350), hv.opts.HLine(line_dash='dashed'))

7.2.3. 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|\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|\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 o y/o iterativas, e.g. gradiente descendente o expectation maximization. Veremos más sobre esto en el futuro

7.3. Una nota sobre estimadores sesgados e insesgados

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

\[ \mathbb{E}[\hat \theta] = \theta, \]

entonces \(\hat \theta\) es un estimador insesgado de \(\theta\)

¿Es el MLE de \(\mu\) insesgado?

\[\begin{split} \begin{split} \mathbb{E}[\hat \mu] &= \mathbb{E} \left[ \frac{1}{N} \sum_{i=1}^N x_i \right] \\ &= \frac{1}{N} \sum_{i=1}^N \mathbb{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\) insesgado?

\[\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} \mathbb{E}[\hat \sigma^2] &= \frac{N-1}{N^2} \sum_{i=1}^N \mathbb{E} [x_i^2] - \frac{1}{N^2} \sum_{i=1}^N \sum_{j \neq i} \mathbb{E} [x_i] \mathbb{E} [x_j] \\ &= \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 \]

7.4. Condiciones de optimalidad e incertidumbre de los MLEs

Asumiendo que los datos vienen del modelo especificado entonces el MLE es

Consistente: El estimador converge al parámetro real a medida que los datos aumentan

\[ \lim_{N\to \infty} \hat \theta = \theta \]

Asintóticamente normal: La distribución del estimador se aproxima a una Gaussiana centrada en el valor real del parámetro

\[ \lim_{N\to \infty} p(\hat \theta) = \mathcal{N}(\hat \theta | \theta, \sigma_\theta^2), \]

lo cual es consecuencia del teorema central del límite

Para \(\{X_i\}, i=1,\ldots,N\) i.i.d., con \(\mathbb{E}[X] < \infty\) y \(\text{Var}[X] < \infty\) entonces

\[ \lim_{N\to\infty} \sqrt{N} (\bar X - \mathbb{E}[X]) = \mathcal{N}(0, \sigma^2) \]

Nota

Dado que el MLE tiene distribución asintotica Gaussiana entonces el cociente entre log verosimilitudes tiene distribución asintótica chi-cuadrado

Mínima varianza: El estimador alcanza la varianza mínima teórica dada por la cota de Cramer-Rao:

\[ \sigma_{nm}^2 = \left (- \frac{d^2 \log \mathcal{L} (\theta)}{d\theta_n \theta_m} \bigg\rvert_{\theta = \hat\theta}\right)^{-1} \]

La cota corresponde al inverso de la “información de Fisher”, i.e la segunda derivada de \(- \log L\) con respecto a \(\theta\)

Nota

  • \(\sigma_{nm}^2\) es la varianza mínima que puede alcanzar un estimador insesgado

  • \(\sigma_{nn}^2\) entrega las barras de error marginales

  • If \(\sigma_{nm} \neq 0 ~~ n\neq m\), entonces los errores son correlacionados, i.e algunas combinaciones de los parámetros puede ser mejor determinadas que otras

7.4.1. Ejemplo: Cota de Cramer-Rao para el MLE de \(\mu\)

Considerando la verosimilitud Gaussiana del ejemplo anterior

\[ \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 \]

¿Cuál es la incertidumbre del MLE de \(\mu\)?

En este caso la cota de Cramer-rao es

\[\begin{split} \begin{split} \sigma_{\hat\mu}^2 &= \left (- \frac{d^2 \log \mathcal{L}(\mu, \sigma^2)}{d\mu^2} \bigg\rvert_{\mu=\hat\mu}\right)^{-1} \\ &= \left (- \frac{1}{\sigma^2} \frac{d}{d\mu} \sum_{i=1}^N (x-\mu) \bigg\rvert_{\mu=\hat\mu}\right)^{-1} \\ &= \left ( \frac{N}{\sigma^2} \bigg\rvert_{\mu=\hat\mu}\right)^{-1} = \frac{\sigma^2}{N} \end{split} \end{split}\]

una expresión conocida como el error estándar de la media

Luego tenemos

\[ p(\hat \mu) \to \mathcal{N} \left(\mu, \frac{\sigma^2}{N} \right) \]

Veamos como cambia el MLE y su incertidumbre 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)
(mu_plot*se_plot*mu_real_plot).opts(hv.opts.Curve(logx=True, width=500))

7.5. Tests de hipótesis basados en la verosimilitud

Considerando la distribución asintótica qeu se mostró antes podemos formular un test de hipótesis apra el MLE de \(\theta\). A continuación se presentan dos tests clásicos.

7.5.1. Test de Wald

Supongamos que deseamos probar si

\[ \mathcal{H}_0: \theta = \theta_0 \]
\[ \mathcal{H}_A: \theta \neq \theta_0 \]

Bajo la hipótesis nula podemos escribir

\[ W = \frac{(\hat \theta - \theta_0)^2}{\left (- \frac{d^2 \log \mathcal{L} (\theta)}{d\theta^2} \bigg\rvert_{\theta = \hat\theta}\right)^{-1}} = (\hat \theta - \theta_0)^2 \sigma_{\hat \theta}^2 \to \chi^2_1 \]

El estadística de test tiene una distribución \(\chi^2\) con un grado de libertad

Si \(W\) es mayor que el percentil \((1-\alpha)100\%\) de \(\chi^2_1\) entonces rechazamos la hipótesis nula

7.5.2. Test de Wilks

Supongamos que deseamos probar si

\[ \mathcal{H}_0: \theta = \theta_0 \]
\[ \mathcal{H}_A: \theta =\theta_1 \]

Podemos escribir el cociente entre las verosimilitudes:

\[ \lambda(\mathcal{D}) = \frac{\mathcal{L}(\theta_0|\mathcal{D})}{\mathcal{L}(\theta_1|\mathcal{D})} \]

Que bajo la hipótesis nula tiene distribución asintótica

\[ -2 \log \lambda(\mathcal{D}) \to \chi^2_1 \]

Luego si \(-2 \log \lambda(\mathcal{D})\) es mayor que el percentil \((1-\alpha)100\%\) de \(\chi^2_1\) entonces rechazamos la hipótesis nula

7.6. Criterios para comparación de modelos

¿Cómo comparamos modelos que tienen distinto número de parámetros?

En general mientras más parámetros tenga un modelo mayor es su capacidad de ajuste (y de sobreajuste). La verosimilitud no toma en consideración la complejidad de un modelo, es decir su cantidad de parámetros libres.

Una opción para comparar modelos que toma en cuenta la complejidad es el criterio de información de Akaike (AIC). Para un modelo con \(k\) parámetros que se ajuste con \(N\) datos el AIC se define como

\[ \text{AIC} = -2 \log \mathcal{L}(\hat \theta) + 2k + \frac{2k(k+1)}{N-k-1}, \]

el cual buscamos minimizar. El AIC combina la verosimilitud y la complejidad del modelo.

Nota

Esta es una idea que está relacionada al concepto de regularización, el cual revisaremos en detalle más adelante

anterior

6. Inferencia con modelos paramétricos y no paramétricos

siguiente

8. Modelamiento con métodos no-paramétricos