Rincón de Práctica R

MLE usando R

---

Ejemplo 1.

Sea \(\epsilon_i\sim i.i.d.\, N(\mu,\sigma^2)\).

1. Planear la función log-verosimilitud.

2. Realizar en el computador simulación de 100 datos con \(\mathbb{N}(\mu_0,\sigma_0^2)=\mathbb{N}(1.7,3^2)\) y estimación ML.

Respuesta propuesta. Como vimos, la log-likelihood function estará dada por

\[\ell(\boldsymbol{\theta})=-\frac{n\cdot log(2\pi)}{2}-\frac{n\cdot log(\sigma^2)}{2}-\sum_{i=1}^N{\frac{(\epsilon_i-\mu)^2}{2\cdot \sigma^2}}\]

Así, una propuesta a la parte computacional en R sería la siguiente:

# 1. Simular datos
set.seed(123)

n   <- 100
eps <- rnorm(n, 1.7, 3)

summary(eps)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-5.2275  0.2184  1.8853  1.9712  3.7755  8.2620 

# 2. Escribir la función log-likelihood
logL<- function(theta) {
  mu    <- theta[1]
  sigma <- theta[2]
  return( sum(-log(2*pi*sigma^2)/2 - (eps-mu)^2/(2*sigma^2)) )
}

# 3. Optimizar la log-likelihood
#    Se puede usar optim directamente. Sin embargo, 
#    la libreria maxLik se encarga de computos adicionales 
#    por ejemplo, de la matriz hesiana para los s.e.
library(maxLik)
theta.guess <- c(mu = 1, sigma = 1)
mle         <- maxLik(logLik=logL, start=theta.guess)
summary(mle)

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/

--------------------------------------------
Maximum Likelihood estimation
Newton-Raphson maximisation, 8 iterations
Return code 8: successive function values within relative tolerance limit (reltol)
Log-Likelihood: -242.1305 
2  free parameters
Estimates:
      Estimate Std. error t value  Pr(> t)    
mu      1.9712     0.2724   7.235 4.65e-13 ***
sigma   2.7247     0.1925  14.151  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
--------------------------------------------

---

Ejemplo 2.

Bajo normalidad, MLE y OLS dan el mismo vector estimado de parámetros para el modelo lineal.

Sea,

\[Y=\boldsymbol{X}\boldsymbol{\beta}+u \hspace{0.5cm};\hspace{0.5cm}u_i\sim N(0,\sigma^2)\]

1. Planear la función log-verosimilitud

2. Encontrar \(\hat{\boldsymbol{\theta}}_{MLE}=(\hat{\boldsymbol{\beta}},\hat{\sigma}^2)\)

3. Encontrar la matriz de información

4. Realizar en el computador simulación de 100 datos para \(y_i=\beta_0+\beta_1x_{1i}+\beta_2x_{2i}+u_i\), donde \(x_1\sim U[0,10]\), \(x_2\sim\chi^2_1\), \(u_i\sim\mathbb{N}(0,\sigma^2=3^2)\) y \(\boldsymbol{\beta}'=(\beta_0,\beta_1,\beta_2)'=(1,0.75,3.2)\). Realizar estimación ML de los parámetros y comparar con OLS.

Respuesta sugerida.

Como vimos, para \(\boldsymbol{x}=(1,X_1,x_2)'\) y \(u_i=y_i-\boldsymbol{x}_i'\boldsymbol{\beta}\) con \(\boldsymbol{\beta}'=(\beta_0,\beta_1,\beta_2)'\) se tiene que

\[\ell(\boldsymbol{\theta})=-\frac{n\cdot log(2\pi\sigma^2)}{2}-\frac{(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta})'(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta})}{2\sigma^2}\]

y a partir de las condiciones de primer orden (F.O.C.) evaluadas en el estimador ML, se obtiene: $\(\hat{\boldsymbol{\beta}}_{MLE}=(X'X)^{-1}(X'Y)\hspace{0.5cm},\hspace{0.5cm}\hat{\sigma}^2_{MLE}=(Y-X\hat{\beta})'(Y-X\hat{\beta})\cdot n^{-1}\)$

Notar que la expresión para \(\beta\) corresponde a la misma que se obtiene por OLS.

Además, después de resolver, se obtiene la siguiente matriz de información

\[\begin{split}I(\boldsymbol{\beta}_0)= \left[\begin{array}{cc} \frac{1}{\sigma_0^2}(X'X) & 0\\ 0 & \frac{N}{2\sigma_0^4} \end{array}\right]\end{split}\]

Así, una propuesta para el código de programación en R es la siguiente:

# 1. Simulación de datos
n  <- 100 
s0 <- 3;
b0 <- c(1, 0.75, 3.2);

x1 <- runif(n,0,10);
x2 <- rchisq(n,1);
X  <- cbind(rep(1,n) , x1 , x2);
u  <- rnorm(n,0,s0);

y  <- X%*%b0 + u

# 2. Escribir la función log-likelihood
logLikFun <- function(param) {
    sigma <- param[1]
    b     <- param[-1]
    mu    <- X%*%b
    sum(dnorm(y, mean = mu, sd = sigma, log = TRUE))
}

# 3. Optimizar la log-likelihood
.guess <- c(sigma = 1, beta1 = 1, beta2 = 1, beta3=1)
library(maxLik)
mle     <- maxLik(logLik = logLikFun, start = .guess)
summary(mle)


# Comparación con OLS:
summary(lm(y~X-1))

--------------------------------------------
Maximum Likelihood estimation
Newton-Raphson maximisation, 10 iterations
Return code 1: gradient close to zero (gradtol)
Log-Likelihood: -246.4867 
4  free parameters
Estimates:
      Estimate Std. error t value Pr(> t)    
sigma  2.84604    0.20122  14.144  <2e-16 ***
beta1  0.53229    0.57205   0.930   0.352    
beta2  0.89770    0.09746   9.211  <2e-16 ***
beta3  2.86173    0.18127  15.787  <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
--------------------------------------------

Call:
lm(formula = y ~ X - 1)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.8445 -1.8828  0.0685  2.3769  7.8375 

Coefficients:
    Estimate Std. Error t value Pr(>|t|)    
X    0.53229    0.58040   0.917    0.361    
Xx1  0.89770    0.09886   9.081  1.3e-14 ***
Xx2  2.86173    0.18413  15.542  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.89 on 97 degrees of freedom
Multiple R-squared:  0.916,	Adjusted R-squared:  0.9134 
F-statistic: 352.6 on 3 and 97 DF,  p-value: < 2.2e-16