Implementar la puntuación de Fisher para la regresión lineal

Question

Sé que existe una solución analítica para el siguiente problema (OLS). Como trato de aprender y comprender los principios y fundamentos de MLE, implementé el algoritmo de puntuación de Fisher para un modelo de regresión lineal simple.

$$ y = X\beta + \epsilon\\ \epsilon\sim N(0,\sigma^2) $$

La loglikelihood para $\sigma^2$ y $\beta$ está dada por: $$ -\frac{N}{2}\ln(2\pi)-\frac{N}{2}\ln(\sigma^2)-\frac{1}{2\sigma^{2}}(y-X\beta)^{'}(y-X\beta) $$

Para calcular la función de puntuación $S(\theta)$ , donde $\theta$ es el vector de parámetros $(\beta,\sigma^{2})^{'}$ Tomo las primeras derivadas parciales con respecto a $\beta$ y $\sigma^{2}$ : \begin {align} \frac { \partial L}{ \partial \beta } &= \frac {1}{ \sigma ^{2}}(y-X \beta )^{'}X \\ [5pt] \frac { \partial L}{ \partial \sigma ^2} &= - \frac {N}{ \sigma ^{2}}+ \frac {1}{2 \sigma ^{4}}(y-X \beta )^{'}(y-X \beta ) \end {align}

Entonces el algoritmo de puntuación de Fisher se implementa como: $$ \theta_{j+1} = \theta_{j} - (S(\theta_{j})S(\theta_{j})^{'})S(\theta_{j}) $$

Tenga en cuenta que el siguiente código es una implementación muy ingenua (sin regla de parada, etc.)

library(MASS)
x <- matrix(rnorm(1000), ncol = 2)
y <- 2 + x %*% c(1,3) + rnorm(500)

fisher.scoring <- function(y, x, start = runif(ncol(x)+1)){
    n <- nrow(x)
    p <- ncol(x)
    theta <- start
    score <- rep(0, p+1)
    for (i in 1:1e5){
        # betas
        score[1:p] <- (1/theta[p+1]) * t((y - x%*%theta[1:p])) %*% x
        # sigma
        score[p+1] <- -(n/theta[p+1]) + (1/2*theta[p+1]^2) * crossprod(y - x %*% theta[1:p])
        # new
        theta <- theta - MASS::ginv(tcrossprod(score)) %*% score
    }
    return(theta)
}

# Gives the correct result 
lm.fit(cbind(1,x), y)$coefficients
# Does not give the correct result
fisher.scoring(y, cbind(1,x))
# Even if you start with the correct values
fisher.scoring(y, cbind(1,x), start=c(2,1,3,1))

Mi pregunta

¿Qué me he perdido? ¿Dónde está mi error?

Answer 1

He corregido el código según las sugerencias de @Randel. Ahora funciona, excepto $\sigma^2$ tarda mucho tiempo en converger.

Aquí está el código:

fisher.scoring <- function(y, x, start = runif(ncol(x)+1)){
  n <- nrow(x)
  p <- ncol(x)
  theta <- start
  score <- rep(0, p+1)
  for (i in 1:1000){
    # betas
    score[1:p] <- (1/theta[p+1]) * t((y - x%*%theta[1:p])) %*% x
    # sigma
    score[p+1] <- -(n/(2*theta[p+1])) + 
                   (1/(2*theta[p+1]^2)) * 
                   crossprod(y - x %*% theta[1:p])
    # new
    hessMat <- matrix(0,ncol=p+1,nrow = p+1)
    for(j in 1:n)
    {
      # Estimate derivative of likelihood for each observation
      estVec <- c((1/theta[p+1]) * 
                t((y[j] - x[j,]%*%theta[1:p])) %*% x[j,],
                -(n/(2*theta[p+1])) + (1/(2*theta[p+1]^2)) *
                crossprod(y[j] - x[j,] %*% theta[1:p]))
      # Add them up as suggested to get an estimate of the Hessian.
      hessMat <- hessMat + estVec%*%t(estVec)
    }
    theta <- theta + MASS::ginv(hessMat) %*% score
  }
  return(theta)
}

Ahora cuando ejecuto el código obtengo lo siguiente:

> lm.fit(cbind(1,x), y)$coefficients
       x1        x2        x3 
2.0136134 0.9356782 2.9666921 
> fisher.scoring(y, cbind(1,x))
          [,1]
[1,] 2.0136126
[2,] 0.9356782
[3,] 2.9666917
[4,] 0.6534185

Lo he vuelto a ejecutar con 5000 iteraciones en lugar de 1000 y me sale:

> fisher.scoring(y, cbind(1,x))
          [,1]
[1,] 2.0136133
[2,] 0.9356782
[3,] 2.9666920
[4,] 0.8962295

Espero que esto ayude. $\sigma^2$ parece converger muy lentamente, no sé por qué, ¡supongo que eso es material para otra pregunta!

EDITAR: Aquí se puede ver la convergencia de los parámetros. El número de iteraciones está en una escala logarítmica. Los coeficientes de regresión tardan entre 30 y 40 iteraciones en converger, aunque el $\beta_1$ parámetro se sobrepasa y luego vuelve a bajar, (no esperaba ver eso). El $\sigma^2$ converge bastante rápido como los demás, y luego la convergencia se ralentiza mucho. De momento no tengo ni idea de por qué ocurre esto.

EDIT2: Hubo un error en la actualización de $\sigma^2$ . Véase la respuesta de Randel para la corrección.

Answer 2

La respuesta corta es que hay un error en el código de @guðmundur-einarsson.

Para cada observación, la función de puntuación para $\sigma^2$ es $$\frac{\partial L_i}{\partial \sigma^2} = -\frac{1}{2\sigma^{2}}+\frac{1}{2\sigma^{4}}(y_i-X_i\beta)^{'}(y_i-X_i\beta),$$ no $\frac{\partial L_i}{\partial \sigma^2} = -\frac{{\color{red}N}}{2\sigma^{2}}+\frac{1}{2\sigma^{4}}(y_i-X_i\beta)^{'}(y_i-X_i\beta)$ . Así que sólo sustituir n en estVec con 1 . Aunque no es evidente, hay algunas pistas.

• La estimación de $\sigma^2$ extrañamente no cambia mucho a través de las iteraciones, aunque muestra alguna tendencia cuando ylim se ajusta automáticamente para ser pequeño. Este pequeño cambio se debe a que cada movimiento se divide por unos $N$ .
• Como señala @guðmundur-einarsson en la pregunta relacionada:

La convergencia es bastante rápido al principio y luego se estanca y continúa muy lentamente cuando los parámetros de regresión convergen.

Cuando $\beta$ no converge (en la ronda 4 del eje x en las figuras anteriores) , $\sigma^2$ depende de $\beta$ Así que está cambiando más rápido en comparación con el después $\beta$ converge. La escala logarítmica de los tiempos de iteración también parece producir algún efecto artificial.

• La razón por la que $\beta$ aún puede converger con el error puede explicarse porque su solución de forma cerrada no depende de $\sigma^2$ en absoluto.

Las siguientes cifras se basan en 100 iteraciones y la semilla aleatoria se fija en 1.

La lección que aprendí de esto es mantener la depuración, aunque es relativamente más fácil en este pequeño ejemplo.