Hace unos años escribí un artículo sobre esto para mis alumnos, así que puedo intentar reescribir esas explicaciones aquí. Veré el IRLS (iteratively reweighted least squares) a través de una serie de ejemplos de complejidad creciente. Para el primer ejemplo necesitamos el concepto de familia localización-escala. Sea $f_0$ sea una función de densidad centrada en cero en algún sentido. Podemos construir una familia de densidades definiendo $$ f(x)= f(x;\mu,\sigma)= \frac{1}{\sigma} f_0\left(\frac{x-\mu}{\sigma}\right) $$ donde $\sigma > 0$ es un parámetro de escala y $\mu$ es un parámetro de localización. En el modelo de error de medición, en el que el término de error se modela habitualmente como una distribución normal, podemos utilizar en lugar de esa distribución normal una familia de escala de localización como la construida anteriormente. Cuando $f_0$ es la distribución normal estándar, la construcción anterior da la $\text{N}(\mu, \sigma)$ familia.
Ahora utilizaremos el IRLS en algunos ejemplos sencillos. Primero encontraremos los estimadores ML (máxima verosimilitud) en el modelo $$ Y_1,Y_2,\ldots,Y_n \hspace{1em} \text{i.i.d} $$ con la densidad $$ f(y)= \frac{1}{\pi} \frac{1}{1+(y-\mu)^2},\hspace{1em} y\in{\mathbb R}, $$ la distribución de Cauchy la familia de localizaciones $\mu$ (por lo que se trata de una familia de ubicación). Pero primero algo de notación. El estimador de mínimos cuadrados ponderados de $\mu$ viene dado por $$ \mu^{\ast} = \frac{\sum_{i=1}^n w_i y_i} {\sum_{i=1}^n w_i}. $$ donde $w_i$ son algunos pesos. Veremos que el estimador ML de $\mu$ puede expresarse de la misma forma, con $w_i$ alguna función de los residuos $$ \epsilon_i = y_i-\hat{\mu}. $$ La función de verosimilitud viene dada por $$ L(y;\mu)= \left(\frac{1}{\pi}\right)^n \prod_{i=1}^n \frac{1}{1+(y_i-\mu)^2} $$ y la función loglikelihood viene dada por $$ l(y)= -n \log(\pi) - \sum_{i=1}^n \log\left(1+(y_i-\mu)^2\right). $$ Su derivada con respecto a $\mu$ es $$ \begin{eqnarray} \frac{\partial l(y)}{\partial \mu}&=& 0-\sum \frac{\partial}{\partial \mu} \log\left(1+(y_i-\mu)^2\right) \nonumber \\ &=& -\sum \frac{2(y_i-\mu)}{1+(y_i-\mu)^2}\cdot (-1) \nonumber \\ &=& \sum \frac{2 \epsilon_i}{1+\epsilon_i^2} \nonumber \end{eqnarray} $$ donde $\epsilon_i=y_i-\mu$ . Escriba a $f_0(\epsilon)= \frac{1}{\pi} \frac{1}{1+\epsilon^2}$ y $f_0'(\epsilon)=\frac{1}{\pi} \frac{-1\cdot 2 \epsilon}{(1+\epsilon^2)^2}$ obtenemos $$ \frac{f_0'(\epsilon)}{f_0(\epsilon)} = \frac{\frac{-1 \cdot2\epsilon}{(1+\epsilon^2)^2}} {\frac{1}{1+\epsilon^2}} = -\frac{2\epsilon}{1+\epsilon^2}. $$ Encontramos $$ \begin{eqnarray} \frac {\partial l(y)} {\partial \mu} & =& -\sum \frac {f_0'(\epsilon_i)} {f_0(\epsilon_i)} \nonumber \\ &=& -\sum \frac {f_0'(\epsilon_i)} {f_0(\epsilon_i)} \cdot \left(-\frac{1}{\epsilon_i}\right) \cdot (-\epsilon_i) \nonumber \\ &=& \sum w_i \epsilon_i \nonumber \end{eqnarray} $$ donde utilizamos la definición $$ w_i= \frac{f_0'(\epsilon_i)} {f_0(\epsilon_i)} \cdot \left(-\frac{1}{\epsilon_i}\right) = \frac{-2 \epsilon_i} {1+\epsilon_i^2} \cdot \left(-\frac{1}{\epsilon_i}\right) = \frac{2}{1+\epsilon_i^2}. $$ Recordando que $\epsilon_i=y_i-\mu$ obtenemos la ecuación $$ \sum w_i y_i = \mu \sum w_i, $$ que es la ecuación de estimación del IRLS. Obsérvese que
- Los pesos $w_i$ son siempre positivos.
- Si el residuo es grande, damos menos peso a la observación correspondiente.
Para calcular el estimador ML en la práctica, necesitamos un valor inicial $\hat{\mu}^{(0)}$ podríamos utilizar la mediana, por ejemplo. Con este valor calculamos los residuos $$ \epsilon_i^{(0)} = y_i - \hat{\mu}^{(0)} $$ y pesos $$ w_i^{(0)} = \frac{2}{1+\epsilon_i^{(0)} }. $$ El nuevo valor de $\hat{\mu}$ viene dado por $$ \hat{\mu}^{(1)} = \frac{\sum w_i^{(0)} y_i} {\sum w_i^{(0)} }. $$ Siguiendo así, definimos $$ \epsilon_i^{(j)} = y_i- \hat{\mu}^{(j)} $$ y $$ w_i^{(j)} = \frac{2}{1+\epsilon_i^{(j)} }. $$ El valor estimado al paso $j+1$ del algoritmo pasa a ser $$ \hat{\mu}^{(j+1)} = \frac{\sum w_i^{(j)} y_i} {\sum w_i^{(j)} }. $$ Continuando hasta la secuencia $$ \hat{\mu}^{(0)}, \hat{\mu}^{(1)}, \ldots, \hat{\mu}^{(j)}, \ldots $$ converge.
Ahora estudiamos este proceso con una familia de ubicación y escala más general, $f(y)= \frac{1}{\sigma} f_0(\frac{y-\mu}{\sigma})$ con menos detalles. Sea $Y_1,Y_2,\ldots,Y_n$ ser independiente con la densidad anterior. Defina también $ \epsilon_i=\frac{y_i-\mu}{\sigma}$ . La función de verosimilitud es $$ l(y)= -\frac{n}{2}\log(\sigma^2) + \sum \log(f_0\left(\frac{y_i-\mu}{\sigma}\right)). $$ Escribir $\nu=\sigma^2$ Obsérvese que $$ \frac{\partial \epsilon_i}{\partial \mu} = -\frac{1}{\sigma} $$ y $$ \frac{\partial \epsilon_i}{\partial \nu} = (y_i-\mu)\left(\frac{1}{\sqrt{\nu}}\right)' = (y_i-\mu)\cdot \frac{-1}{2 \sigma^3}. $$ Cálculo de la derivada de loglikelihood $$ \frac{\partial l(y)}{\partial \mu} = \sum \frac{f_0'(\epsilon_i)}{f_0(\epsilon_i)}\cdot \frac{\partial \epsilon_i}{\partial \mu} = \sum\frac{f_0'(\epsilon_i)}{f_0(\epsilon_i)}\cdot\left(-\frac{1}{\sigma}\right)= -\frac{1}{\sigma}\sum\frac{f_o'(\epsilon_i)}{f_0(\epsilon_i)}\cdot \left(-\frac{1}{\epsilon_i}\right)(-\epsilon_i) = \frac{1}{\sigma}\sum w_i \epsilon_i $$ e igualándola a cero se obtiene la misma ecuación de estimación que en el primer ejemplo. A continuación, la búsqueda de un estimador para $\sigma^2$ : $$ \begin{eqnarray} \frac{\partial l(y)}{\partial \nu} &=& -\frac{n}{2}\frac{1}{\nu} + \sum\frac{f_0'(\epsilon_i)}{f_0(\epsilon_i)}\cdot \frac{\partial \epsilon_i}{\partial\nu} \nonumber \\ &=& -\frac{n}{2}\frac{1}{\nu}+\sum\frac{f_0'(\epsilon_i)}{f_0(\epsilon_i)} \cdot \left(-\frac{(y_i-\mu)}{2\sigma^3}\right) \nonumber \\ &=& -\frac{n}{2}\frac{1}{\nu} - \frac{1}{2}\frac{1}{\sigma^2} \sum\frac{f_0'(\epsilon_i)}{f_0(\epsilon_i)}\cdot \epsilon_i\nonumber \\ &=& -\frac{n}{2}\frac{1}{\nu}-\frac{1}{2}\frac{1}{\nu} \sum\frac{f_0'(\epsilon_i)}{f_0(\epsilon_i)}\cdot \left(-\frac{1}{\epsilon_i}\right) (-\epsilon_i)\cdot\epsilon_i\nonumber \\ &=& -\frac{n}{2}\frac{1}{\nu}+\frac{1}{2}\frac{1}{\nu}\sum w_i \epsilon_i^2 \stackrel{!}{=} 0. \nonumber \end{eqnarray} $$ lo que conduce al estimador $$ \hat{\sigma^2} = \frac{1}{n}\sum w_i (y_i-\hat{\mu})^2. $$ El algoritmo iterativo anterior también puede utilizarse en este caso.
A continuación damos un ejemplo numérico usando R, para el modelo exponencial doble (con escala conocida) y con datos y <- c(-5,-1,0,1,5) . Para estos datos el valor verdadero del estimador ML es 0. El valor inicial será mu <- 0.5 . Una pasada del algoritmo es
iterest <- function(y, mu) {
w <- 1/abs(y - mu)
weighted.mean(y, w)
}
con esta función se puede experimentar haciendo las iteraciones "a mano" Entonces el algoritmo iterativo se puede hacer por
mu_0 <- 0.5
repeat {mu <- iterest(y, mu_0)
if (abs(mu_0 - mu) < 0.000001) break
mu_0 <- mu }
Ejercicio : Si el modelo es un $t_k$ con parámetro de escala $\sigma$ muestran que las iteraciones vienen dadas por el peso $$ w_i = \frac{k + 1}{k + \epsilon_i^2}. $$ Ejercicio : Si la densidad es logística, demuestre que los pesos vienen dados por $$ w(\epsilon) = \frac{ 1-e^\epsilon}{1+e^\epsilon} \cdot - \frac{1}{\epsilon}. $$
De momento lo dejo aquí, continuaré con este post.
