Simulacion MCO vs Máxima Verosimilitud Regresión Lineal
Katherin Quiñones
Agosto 2022
El método de mínimos cuadrados ordinarios suele ser comunmente utilizado para estimar los parametros en un modelo de regresion lineal independiente de la forma de distribución de los errores \(\epsilon\).
Para lograr hacer inferencias a partir de los modelos construidos, utilizar mco no es suficiente. Entre los procedimientos adicionales que se pueden requerir calcular se encuentran por ejemplo la prueba de hipótesis o construccion de intervalos de confianza, en donde los errores se distribuyen normalmente. Es aquí donde el método por máxima verosimilitud coge relevancia.
Por consiguiente, se decide hacer una comparación de los resultados entre mco y máxima verosimilitus, para un grupo de valores aleatorios.
Simulación MCO
Paso 1:
Para la simulación vamos a generar un conjunto de datos basandonos en un modelo con parametros previamente definidos: \(Y_i=5+(3X_i)+ε_i\), en donde podemos notar que \(β_0=5\), \(β_1=3\) y \(εi \sim Normal(0,10)\).
x=runif(30,1,30) ## Valores iniciales arbitrarios para X
y = 5+(3*x) + rnorm(n = 30,mean = 0,sd = 10 ) ##Valores de Y de acuerdo al modelo
plot(x,y, xlim = c(0,30), ylim = c(0,100), las=1, pch=19, col="#b0394a")
datos=data.frame(x,y)
head(datos,10)## x y
## 1 5.254957 13.51626
## 2 22.677414 69.65715
## 3 29.035909 95.49465
## 4 13.338048 38.96966
## 5 28.396207 86.52463
## 6 27.753702 78.74299
## 7 11.195036 47.13603
## 8 6.910198 12.98350
## 9 2.312347 11.88881
## 10 16.409397 43.55568Paso 2:
Ahora vamos a asumir valores para la estimación de \(β_0=4\) y \(β_1=3.5\), posteriormente graficamos la linea y obtenemos el valor de SCE.
b0=4 # intercepto estimado
b1=2.5 # pendiente estimada
x1=1:30
y1 = b0 + b1 * x1
plot(x,y, xlim = c(0,30), ylim = c(0,100), las=1,pch=19 ,col='blue')
lines(x1, y1 , col = 2, type = "l")
grid()
Paso 3:
Se realiza el Paso 2 para otros posibles valores de los parametros y se guarda el valor de SCE (Suma de Cuadrados del Error)
beta0_est= seq(3,7,0.5) # rango de valores al rededor de beta0 = 5
beta1_est=seq(1,5,0.3) # rango de valores al rededor de beta1 = 3
betas=expand.grid(beta0_est,beta1_est)
names(betas)=c("beta0_est","beta1_est")
SCE=array(NA,dim(betas)[1])
plot(x,y, ylim = c(0,100))
for(i in 1:dim(betas)[1]){
y_est=betas$beta0_est[i] + (betas$beta1_est[i]*x)
lines(x,y_est,col='red')
SCE[i]=(y_est-y)^2
}
resultados=data.frame(betas,SCE)Paso 4:
Graficar e interpretar la relación entre los posibles Parametros (betas) y la SCE.
library(plotly)
plot_ly(x=resultados$beta0_est,
y=resultados$beta1_est,
z=resultados$SCE,
size=.5)Acontinuación se hace la construcción del modelo de regresión lineal el cual mostrará los betas que se encontraron anteriormente por MCO.
modelo1=lm(y~x)
summary(modelo1)##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -23.4120 -5.8771 -0.6516 6.2420 19.4608
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.6615 3.4509 0.771 0.447
## x 2.9793 0.1932 15.417 3.28e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.365 on 28 degrees of freedom
## Multiple R-squared: 0.8946, Adjusted R-squared: 0.8909
## F-statistic: 237.7 on 1 and 28 DF, p-value: 3.28e-15Simulación Máxima Verosimilitud
La función de verosimilitud se determina con la distribución conjunta de las observaciones.En este caso se tendrán tres parametros desconocidos $ _0 , _1 ,y , ^2$. Después de los cálculos en el cual se hace el logaritmo natura de la funcion L, obtenemos el siguiente resultado:
\(ln \, L(y_i , x_i , \beta_0 , \beta_1 , \sigma^2 )= -(\frac{n}{2})* \ln 2\pi- (\frac{n}{2}) * \ln \sigma^2 -(\frac{1}{2\sigma^2}) \sum_{i=1}^{n} (y_i - \beta_0 - \beta_1 x_i )^2\)
L=array(1,dim(betas)[1])
n=dim(betas)[1]
for(i in 1:dim(betas)[1]){
y_est=betas$beta0_est[i] + (betas$beta1_est[i]*x)
sigma = SCE[i]/length(x)
L[i]= -(n/2)*log(2*pi)- (n/2)*log(sigma) - (1/(2*sigma))*(y-y_est)^2
}
resultados_mv=data.frame(betas,L)
head(resultados_mv)## beta0_est beta1_est L
## 1 3.0 1 -125.71852
## 2 3.5 1 -113.13649
## 3 4.0 1 -99.15725
## 4 4.5 1 -83.43124
## 5 5.0 1 -65.45877
## 6 5.5 1 -44.48932La gráfica que se encuentra al visualizar los betas y el valor de L, es una gráfica concava hacia abajo, lo cual no indica los mejores betas que hacen que encontremos el minimo error a partir de la máxima verosimilitud.
plot_ly(x=resultados_mv$beta0_est,
y=resultados_mv$beta1_est,
z=resultados_mv$L,
size=.5)min_ml=resultados_mv[resultados_mv$L==max(resultados_mv$L),]
min_ml## beta0_est beta1_est L
## 29 3.5 1.9 517.8128min_mco=resultados[resultados$SCE==min(resultados$SCE),]
min_mco## beta0_est beta1_est SCE
## 29 3.5 1.9 0.001013824Conclusión
Los estimadores de máxima verosimilitud de la ordenada en el origen y la pendiente, \(\beta_0 \, y \,\beta_1\) son idénticos a los obtenidos con mínimos cuadrados.
Casi siempre los estimadores por máxima verosimilitud tienen mejores propiedades estadísticas, que los estimadores por MCO. Los estimadores por ML, son insesgados, y tienen varianza mínima.