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β+ϵϵ∼N(0,σ2)
La loglikelihood para σ2 y β está dada por: −N2ln(2π)−N2ln(σ2)−12σ2(y−Xβ)′(y−Xβ)
Para calcular la función de puntuación S(θ) , donde θ es el vector de parámetros (β,σ2)′ Tomo las primeras derivadas parciales con respecto a β y σ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: θj+1=θj−(S(θj)S(θj)′)S(θ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?
1 votos
Un primer vistazo muestra ∂L∂σ2=−N2σ2+12σ4(y−Xβ)′(y−Xβ) . Y θj+1=θj+(S(θj)S(θj)′)S(θj) Aunque la gente usaría una suma de primeras derivadas al cuadrado para cada observación, no una suma al cuadrado. Un ejemplo es aquí .