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.
Se define como parametros iniciales la regresión \(Y = 5 + 3*X + \epsilon _i\) donde \(\epsilon_1 \sim N(\mu= 0,\sigma^2 = 4)\)
x <- runif(30,1,30) ## Valores iniciales arbitrarios para X
y <- 5+(3*x) + rnorm(n = 30,mean = 0,sd = 2 ) ##Valores de Y de acuerdo al modelo
beta0_est= seq(3,7,0.01) # rango de valores al rededor de beta0 = 5
beta1_est=seq(1,5,0.01) # rango de valores al rededor de beta1 = 3
betas=expand.grid(beta0_est,beta1_est) # Combinacion de Betas1 y Beta0
names(betas)=c("beta0_est","beta1_est") # Dar nombre a las columnas creadasCon los parametros definidos ateriormente, se procesde a calcular la SCE para cada uno de los Betas definidos.
SCE=array(NA,dim(betas)[1])
for(i in 1:dim(betas)[1]){
y_est=betas$beta0_est[i] + (betas$beta1_est[i]*x)
SCE[i]=sum((y_est-y)^2)
}
resultados=data.frame(betas,SCE)
head(resultados)## beta0_est beta1_est SCE
## 1 3.00 1 49281.23
## 2 3.01 1 49259.49
## 3 3.02 1 49237.76
## 4 3.03 1 49216.03
## 5 3.04 1 49194.31
## 6 3.05 1 49172.59Lo que nos genera una gráfica con una parabola concava hacia arriba. Lo que nos dice que esta estimación nos genera un minimo, al identificarlo obtendríamos la estimación más optima por este medio.
library(plotly)
plot_ly(x=resultados$beta0_est,
y=resultados$beta1_est,
z=resultados$SCE,
size=.5)La estimación de los betas seria:
resultados[resultados$SCE==min(resultados$SCE),]## beta0_est beta1_est SCE
## 82705 3.98 3.06 97.66275La función de verosimilitud se determina con la distribución conjunta de las observaciones.En este caso se tendrán tres parametros desconocidos \(\beta_0\) , \(\beta_1\) ,y \(\sigma^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)) =\sum_{i=1}^n \frac{-1}{2}log(2\pi\sigma^2)-\frac{1}{\sigma^2}(\hat{y}_i-y_i)^2 \\ donde \\ y_i = \beta_0 + \beta_1x_i\]
Con los mismos parametros definidos anteriormente, se procede a calcular el resultado de Maxima Verosimilitud para cada uno de los Betas definidos.
L=array(1,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]= sum((-1/2)*(log(2*pi*sigma)-((1/sigma)*((y_est-y)^2))))
}
resultados_mv=data.frame(betas,L)
head(resultados_mv)## beta0_est beta1_est L
## 1 3.00 1 -123.6297
## 2 3.01 1 -123.6231
## 3 3.02 1 -123.6164
## 4 3.03 1 -123.6098
## 5 3.04 1 -123.6032
## 6 3.05 1 -123.5966Lo que nos genera una gráfica concava hacia abajo, lo que nos dice que esta estimación nos genera un máximo, al identificarlo obtendríamos la estimación más optima por este método
plot_ly(x=resultados_mv$beta0_est,
y=resultados_mv$beta1_est,
z=resultados_mv$L,
size=.5)El resultado de la estimación sería:
resultados_mv[resultados_mv$L==max(resultados_mv$L),]## beta0_est beta1_est L
## 82705 3.98 3.06 -30.273Si comparamos resultados, observamos que estimar \(\beta_1\) y \(\beta_2\) por MCO o por Máxima verosimilitud llegamos a la misma estimación