Simulación Máxima Verosimilitud

MAESTRÍA EN CIENCIA DE DATOS UNIVERSIDAD JAVERIANA DE CALI

MÉTODOS ESTADÍSTICOS PARA LA TOMA DE DECISIONES

Carolina Galindres Bernal, Adrian Rodriguez Amaya

Para el desarrollo de la presente simulación utilizaremos un modelo definido con la siguiente ecuación:

\[Y_i = 5+(3X_i)+ε_i\]

Donde:

\[ε_i \textit{ ~ } Normal(0,5)\]

Iniciaremos creando un conjunto de datos arbitrarios y graficaremos un diagrama de dispersión con los datos generados de la siguiente manera:

x = runif(30, 1, 30)
y = 5 + (3*x) + rnorm(n = 30, mean = 0, sd = 5)

ggplot(mapping = aes(x = x, y = y)) + geom_point(color='darkblue', size = 2)

data = data.frame(x, y)
head(data,10)

##            x        y
## 1  18.900760 63.86069
## 2  25.539165 72.14943
## 3  25.203739 82.14224
## 4  21.929459 66.31325
## 5  21.185844 69.23378
## 6  13.282503 51.43282
## 7  15.301204 53.80480
## 8   5.008064 11.71808
## 9   4.386582 14.92313
## 10  4.578089 20.68375

Ahora, asumiremos aleatorios para cada uno de los coeficientes y graficaremos el modelo con base a la ecuación definida anteriormente:

beta_0 = 4
beta_1 = 2.5

x1 = 1:30
y1 = beta_0 + beta_1 * x1

ggplot(mapping = aes(x = x, y = y)) + geom_point(color='darkblue', size = 2) + 
  geom_line(mapping = aes(x = x1, y = y1), color = 'darkorange2')

Luego, construiremos vectores de valores para los betas alrededor de beta0 y beta1 definidos anteriormente. Con estos valores procederemos a cacular la suma de cuadrados del error SCE, así:

beta0_est = seq(3, 7, 0.5)
beta1_est = seq(1, 5, 0.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=3)
  SCE[i] = (y_est - y)^2
}

Con los valores obtenidos de la sumatoria de los cuadrados del error construimos un nuevo data frame de resultados para construir la gráfica. También, obtener los valores mínimos para el SEC:

resultados=data.frame(betas,SCE)
kable(resultados[resultados$SCE == min(resultados$SCE), ], col.names = c('β0', 'β1', 'SCE'), align = 'c', booktabs=TRUE, caption = "Estimadores MCO") %>%
  kable_styling(latex_options = c("striped", "scale_down"), full_width = F, font_size = 15) 

Estimadores MCO

	ß0	ß1	SCE
69	5.5	3.1	0.053668

plot_ly(x=resultados$beta0_est,
        y=resultados$beta1_est,
        z=resultados$SCE,
        size=.5)

Una vez obtenidos los valores para el SCE continuamos con el inicio de los cálculos para la máxima verosimilitud. Así:

MVS=array(NA,dim(betas)[1])

for(i in 1:dim(betas)[1]){
  MVS[i] = -((30/2)*log(10^2))-((30/2)*log(2*3.1416))-((1/2)*SCE[i]/(10^2))
}

resultados2=data.frame(betas,MVS)
kable(resultados2[resultados2$MVS == max(resultados2$MVS), ], col.names = c('β0', 'β1', 'MVS'), align = 'c', booktabs=TRUE, caption = "Estimadores MVS") %>%
  kable_styling(latex_options = c("striped", "scale_down"), full_width = F, font_size = 15) 

Estimadores MVS

	ß0	ß1	MVS
69	5.5	3.1	-96.64601

plot_ly(x=resultados2$beta0_est,
        y=resultados2$beta1_est,
        z=resultados2$MVS,
        size=.5)

De acuerdo con los cálulos anteriores podemos observar que los coeficientes que máximizan la ecuación del modelo conrresponden a los indicados a continuación. Además podemos observar que el resultado de los coeficientes tanto para la simulación de MCO y para la simulación de máxima verosimilitud son iguales, esto se cumple para los casos de regresión simple y bajo el supuesto de normalidad del error:

\[β_0 = 5 \textit{ , } β_1=2.8\]