1. Simulación MCO (Mínimos Cuadrados Ordinarios) - Regresión Lineal

1.1. Paso 1

Para la simulación vamos a generar un conjunto de datos basándonos en un modelo con parámetros previamente definidos: \(Y_{i}=5+(3X_{i})\) + \(\varepsilon_{i}\), en donde podemos notar que \(\beta_{0} = 5, \beta_{1} = 3\) y \(\varepsilon_{i} \approx 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)

1.2. Paso 2

Ahora vamos a asumir valores para la estimación de \(\beta_{0} = 4\) y \(\beta_{1} = 2.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= 3)
lines(x1, y1 , col = 2, type = "l")
grid()

1.3. 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=2)
SCE[i]=(y_est-y)^2
}

resultados=data.frame(betas,SCE)

1.4. Paso 4

library(plotly)

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

En la gráfica se puede observar que existe un punto mínimo en la SCE para una combinación de los valores \(\hat{\beta}_{0}\) y \(\hat{\beta}_{1}\), esto en general lo que nos indica es que el método de MCO busca justamente cuales son esos dos valores tal que nos garantizan una mínima suma de cuadrados del error.

Finalmente se realiza la estimación por MCO utilizando la función lm de R, sabiendo que estos valores corresponden a la combinación de valores de los estimadores que minimizan SCE

modelo1=lm(y~x)
summary(modelo1)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -18.8667  -6.8002  -0.2483   3.7490  21.0713 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    9.242      3.120   2.962  0.00617 ** 
## x              3.002      0.182  16.492 5.96e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.36 on 28 degrees of freedom
## Multiple R-squared:  0.9067, Adjusted R-squared:  0.9033 
## F-statistic:   272 on 1 and 28 DF,  p-value: 5.955e-16

2. Máxima Verosimilitud (MV)

2.1. Paso 1

Para la simulación vamos a generar un conjunto de datos basándonos en un modelo con parámetros previamente definidos: \(Y_{i}=5+(3X_{i})\) + \(\varepsilon_{i}\), en donde podemos notar que \(\beta_{0} = 5, \beta_{1} = 3\) y \(\varepsilon_{i} \approx 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="blue")

datos=data.frame(x,y)
desv=10

head(datos,10)

2.2. Paso 2

Ahora vamos a asumir valores para la estimación de \(\beta_{0} = 4\) y \(\beta_{1} = 2.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=3)
lines(x1, y1 , col = 2, type = "l")
grid()

2.3. 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.2) # rango de valores al rededor de beta0 = 5
beta1_est=seq(1,5,0.2)  # rango de valores al rededor de beta1 = 3 

betas=expand.grid(beta0_est,beta1_est)
names(betas)=c("beta0_est","beta1_est")
L=array(NA,dim(betas)[1])
plot(x,y, ylim = c(0,100))
n=30

parte1=(-n/2)*log(2*pi)-(n/2)*log(desv^2)
parte2=(-1/2*desv^2)

for(i in 1:dim(betas)[1]){
y_est=betas$beta0_est[i] + (betas$beta1_est[i]*x)
lines(x,y_est,col=2)
parte3=(y_est-betas$beta0_est[i]-betas$beta1_est[i]*x)^2
parte3=sum(y_est-y)^2
res=(parte1+parte2*parte3)
L[i]=res
} 

resultados=data.frame(betas,L)

2.4. Paso 4

Graficar e interpretar la relación entre los posibles parametros (betas) y la MEL.

library(plotly)

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

En la gráfica se puede observar que existe un punto mínimo 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 MCO busca justamente cuales son esos dos valores tal que nos garantizan una mínima suma de cuadrados del error.

3. Estimación del parámetro \(\lambda\)

Utilizando la función Fitdistr del paquete MASS de R

3.1. Ejemplo

3.1.1. Consideramos la muestra siguiente de tamaño 50 de suponiendo una variable Poisson de parámetro \(\lambda = 5\)

set.seed(98)
muestra.poisson <- rpois(50, lambda = 5)
muestra.poisson
##  [1]  5  4  4  5  3  4  1  4  6  3  7  7  3  5  4  8  4  4  6  3  6  4  6 11  4
## [26]  7  5  2  8  3  5  4  1  5  6  4  7  7  3  4  6 10  5  4  2  9  1  5  2  2

Estimamos el valor del parámetro \(\lambda = 5\) a partir de la muestra anterior

library(MASS)
fitdistr(muestra.poisson, densfun = "poisson")
##     lambda 
##   4.760000 
##  (0.308545)

La función fitdistr nos da el siguiente de \(\lambda = 4.76\), valor se que aproxima al valor real de \(\lambda = 5\), con un error típico de 0.308545

El estimador de máxima verosimilitud de \(\lambda\) es \(\overline{X}\) con error típico \(\frac{\sqrt{\lambda}}{\sqrt{n}}\)

(estimación.lambda <- mean(muestra.poisson))
## [1] 4.76
(estimacion.error.tipico <- sqrt(estimación.lambda/50))
## [1] 0.308545

Compramos que los valores anteriores coinciden con los dados por la función.

3.1.2. Suponiendo que la muestra anterior es normal con media \(\mu\) y la desviación típica \(\sigma\)

fitdistr(muestra.poisson, densfun = "normal")
##      mean         sd    
##   4.7600000   2.1868699 
##  (0.3092701) (0.2186870)

Dichos valores coinciden con la muestra \(\overline{X}\) y la desviación típica “verdadera” de la muestra considerada

sd(muestra.poisson)*sqrt(49/50)
## [1] 2.18687