Simulación Metodo de Máxima Verosimilidad

Paso 1: Generación del modelo

Tratando de generar una simulación vamos a generar un conjunto de datos que tengan unos parametros adjudicados como con la siguinte función \[Y_{i}=5+(3X_{i})+\varepsilon_{i}\], con \(\beta_{0}\) y \(\beta_{1}\) y \(\varepsilon_{i}\) ~ \(Normal(0,5)\).

X=1:50  ## generamos unos valores para x entre 1 y 50
Y = 5+(3*X) + rnorm(n = 50,mean = 0,sd = 3 ) ## 50 datos con supuestos estadisticos de media 0 y desviación estandar 3

plot(X,Y)

Una vez tenemos los datos generados e posible adicionar el conjunto de datos a un dataframe.

datos=data.frame(X,Y)
head(datos,10)
##     X         Y
## 1   1  5.634616
## 2   2 10.970110
## 3   3 14.354572
## 4   4 17.396457
## 5   5 21.999962
## 6   6 27.486941
## 7   7 26.495180
## 8   8 32.454943
## 9   9 31.804372
## 10 10 35.499775

En una segunda instancia se puede suponer algunos valores para \(\beta_{0}\) y \(\beta_{1}\) por ejemplo \(1\) y \(2\). Podemos graficar la linea recta que terminará ajustando los datos bajo los supuestos valores de \(\beta_{0}=1\) y \(\beta_{1}=2\) junto con los datos originales.

beta0_est=1
beta1_est=2
Y_est=beta0_est + (beta1_est*X)

plot(X,Y)
lines(X,Y_est,col="green")

Ahora podemos establecer muchas mas funciones lineales que tengan diferentes valores de \(\beta_{0}\) y \(\beta_{1}\) usando la generación de secuencias de números.

beta0_est=seq(0,8,0.5)      #Generación desde 0 hasta 8 de 0.5 en 0.5
beta1_est=seq(0,8,0.5)
betas=expand.grid(beta0_est,beta1_est)    #Datos apartir de todas las posibles combinaciones posibles de beta0_est y beta1_est 

names(betas)=c("beta0_est","beta1_est") #Asignar nombres a las variables.

SCE=array(NA,dim(betas)[1])   #Definición de un arreglo para la suma de cuadratica de los errores  
ln_L=array(NA,dim(betas)[1])  #Definición de un arreglo para guardar las funciones de Máxima Verosimilitud
plot(X,Y)

#Con este ciclo se construyen los valores que se le pueden asignar a la función Y. 

for(i in 1:dim(betas)[1]){ # Desde 1 hasta la dimensión del arreglo betas
Y_est=betas$beta0_est[i] + (betas$beta1_est[i]*X) # Cálculo de los valores de cada dato de Y
lines(X,Y_est,col="green")                  # Gráfica de todas las funciones lineales conseguidas.
SCE[i]=sum(Y_est-Y)^2                       # Cálculo de cada una de las diferencias cuadraticas.
ln_L[i]=((-30/2)*log(2*3.1416))-((30/2)*log(5^2))-((1/(2*5^2))*SCE[i])  #Cálculo de la función de verosimilitud.
}

resultados=data.frame(betas,SCE,ln_L)   # Se alamcenan todos los resultados en un arreglo.

Por ultimo se grafica e interpreta la relación entre los posibles Parametros (betas) y la SCE. Mediante un grafico en tres dimensiones.

require(plotly)
plot_ly(x=resultados$beta0_est,y=resultados$beta1_est,z=resultados$ln_L)

En la Grafica se puede observar que existe un punto máximo en la SCE para una combinación de los valores \(\beta_0\) y \(\beta_1\), esto en general lo que nos indica es que el método de máxima verosimilitud busca justamente cuales son esos dos valores tal que nos garantizan una máxima función de verosimilitud.