MEL - intento de simulación-Eliana Poveda

Simulación MCO

#Para la simulación vamos a generar un conjunto de datos basandonos en un modelo con parametros previamente definidos: Yi=5+(3Xi)+εi, en donde podemos notar que β0=5, β1=3 y εi∼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)

##Paso 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
}

#Guardamos como dataframe

resultados=data.frame(betas,SCE)

##Paso 4. #Se hace el grafico para betas y SCE

library(plotly)

## Loading required package: ggplot2

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

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

## No trace type specified:
##   Based on info supplied, a 'scatter3d' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter3d

## No scatter3d mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

#Después de hallar los betas por MCO podemos hacer la somulación por Máxima Verosimilitud

##Paso1 #Se utiliza formula para hallar MLE. se guarda el valor de la probabilidad de la ecuación L de verosimilitud

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)

##Paso.2 #Se gráfican los resultados

plot_ly(x=resultados_mv$beta0_est,
        y=resultados_mv$beta1_est,
        z=resultados_mv$L,
        size=.5)

##Conclusión ##Mientras que el MCO busca saber cuales son los valores que nos garantizan una mínima suma de cuadrados del error, el MLE busca maximizar la función. Tambièn es posible ver que las graficas varian, la de MCO muestra una U, mientras que la de máxima verosimilitud se ve una gráfica concava hacia abajo, lo cual no indica los mejores betas .