2.2. Máxima Verosimilitud
luismor
apuntes para estudiantes
Modelo de máxima verosimilitud
El método de MCO suele ser el más usado para hacer estimaciones lineales, aunque existen otras técnicas también apropiadas. Entre ellas está el método de Máxima Verosimilitud (MV), que supone algunas propiedades teóricas más fuertes que el MCO. Sin adentrarnos en una explicación teórica demasiado profunda, diremos que suponiendo un modelo tal que \(Y_i\) sean independientes y distribuidas de forma:
\[Y_i ∼ N(β_0 + β_1X_i, σ^2)\] Formalmente, esta función puede denotarse como:
\[f(Y_1, Y_2,..., Y_n | β_0 + β_1X_i , σ^2)=\\ f(Y_1 | β_0 + β_1X_i, σ^2)\ f (Y_2|β_0 + β_1X_i, σ^2)... f(Y_n | β_0 + β_1X_i, σ^2)=\\ \frac{1}{σ^n\sqrt{2π}}exp\lbrace-\frac{1}{2}\sum\frac{(Yi − β_0 − β_1X_i)^2}{σ^2}\rbrace\]
siendo la última parte la función de densidad de una función normalmente distribuida con media y varianza dadas, y que se conoce como función de verosimilutd cuando las \(β_i\) son desconocidas, que se escribe como \(FV(β_0,β_1,σ^2)\).
El método de máxima verosimilitud consiste en estimar los parámetros betas de manera que la probabilidad de observar el valor dado \(Y\) sea el mayor posible. Pero al resolver la ecuación anterior, llegaríamos a un resultado tal que los estimadores \(β_i\) de MV son los mismos que los estimadores del modelo MCO. O lo que es lo mismo: la maximización de la probabilidad de un valor Y es equivalente a la minimización de los errores, que es justo el enfoque de los mínimos cuadardos ordinarios.
¿Dónde está la diferencia entonces? La diferencia la encontrarmos en la esperanza de la desviación estándar, que resulta ser insesgada en el caso del modelo MCO, y sesgada para el modelo de MV. De forma concreta:
\[E(\hat σ^2)_{MCO} = [1/(n − 2)]\sum\hat u^2_i\\ E(\hat σ^2)_{MV} = [1/n]\sum\hat u^2_i=σ^2-\frac{2}{n}σ^2\]
Como vemos, el sesgo es hacia abajo, es decir, la esperanza de la varianza en el modelo MV tiende a subestimar el verdadero valor de la varianza. Pero según se aumenta el valor de n, este efecto desaparece, pudiendo decir que hay una convergencia tal que \(E(\hatσ^2)_{MV}=σ^2\) a medida que \(n→∞\).
Diferencias entre MV y MCO
#install.packages("maxLik")
library(maxLik)## Loading required package: miscTools##
## Please cite the 'maxLik' package as:
## Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
##
## If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
## https://r-forge.r-project.org/projects/maxlik/Vamos a generar una muestra aleatoria con un volumen de datos igual a 100 valores. Nuestra variable independiente sige una distribución normal, con media \(\mu = 30\) y una desviación típica de 2. La dependiente se obtiene a partir de la propia independiente:
#Creamos una muestra aleatoria de variables x, y
set.seed(123)
x <- rnorm(100, mean = 30, sd = 2)
u <- runif(100, min = 1, max = 5)
y <- 40+(2*x)+u
plot(x,y)
Para obtener nuestros estimadores MV, debemos crear la función de Máxima Verosimilitud, que se deriva de la función de densidad explicada previamente, aplicando logaritmos:
\[l(θ) = −0.5 N log(2π) − N log σ − 0.5\sum^N_{i=1}\frac{y-β_0 − β_1x_i}{σ^2}\]
mv <- function(param) {
b0 <- param[1]
b1 <- param[2]
sd0 <- param[3]
N <- length(x)
ll <- -0.5*N*log(2*pi) - N*log(sd0) - sum(0.5*(y - b0 - b1*x)^2/sd0^2)
return(ll)
}Dado que el modelo de Máxima Verosimilitud buscará entre todas las opciones la que mejor probabilidad de y obtenga, indicamos un punto de inicio: \(\beta_i = 0\).
start0 = c(0,0,1)
reg <- maxLik(mv, start=start0)## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNs
## Warning in log(sd0): Se han producido NaNssummary(reg)## --------------------------------------------
## Maximum Likelihood estimation
## Newton-Raphson maximisation, 21 iterations
## Return code 1: gradient close to zero (gradtol)
## Log-Likelihood: -157.5611
## 3 free parameters
## Estimates:
## Estimate Std. error t value Pr(> t)
## [1,] 41.86111 1.85247 22.60 <2e-16 ***
## [2,] 2.03589 0.06128 33.22 <2e-16 ***
## [3,] 1.16961 0.08270 14.14 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## --------------------------------------------Para comprobar que los valores obtenidos, es decir \(\beta_0 = 41.86\), \(\beta_1 = 2.03\), son los óptimos, realizamos también el método de MCO. Los resultados deben coincidir.
summary(lm(y~x))##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.92382 -1.05274 -0.08004 0.99649 1.97320
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 41.86111 1.96659 21.29 <2e-16 ***
## x 2.03589 0.06504 31.30 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.181 on 98 degrees of freedom
## Multiple R-squared: 0.9091, Adjusted R-squared: 0.9081
## F-statistic: 979.7 on 1 and 98 DF, p-value: < 2.2e-16Como apuntábamos en la teoría, ambos resultados son los mismos.