Ejercicio 05 - Ejemplo de Máxima verosimilitud

# Máxima verosimilitud

Se tiene una muestra independiente e idénticamente distribuída. \[ y_1,y_2,...,y_N \]

Se asume que los datos son generados por alguna función de densidad paramétrica. En particular, los datos se generan con una función de densidad (o función de masa de probabilidad) conocida: \[ y_i \sim f(y|\theta) \] Considerando el hecho de que los datos son independientes e idénticamente distribuidos, se puede construir la siguiente función de verosimilitud:

\[ L(\theta) = \prod^{N}_{i=1}f(y_i|\theta)\\ \log(L(\theta)) = \sum^{N}_{i=1} \log f(y_i |\theta) \]

Es posible que exista un número infinito de parámetros, se desea elegir los que mejor se “ajusten” a los datos. Para hacer esto es necesario maximizar la función anterior.

Ejemplo utilizando mínimos cuadrados

\[ \begin{align} Y_i = \beta X_i+ U_i && \text{(modelo lineal)}\\ U_i = Y_i - \beta X_i \\ U_i \sim N(0, \sigma^2) && \text{(normalmente distribuido)}\\ Y_i \sim N(\beta X_i, \sigma^2) && \text{(normalmente distribuido)} \end{align} \] La función de densidad esta dada por: \[ f(Y_i|\sigma^2,\beta) \sim \frac{1}{\sqrt{2\pi\sigma^2}}\exp\Big(-\frac{Y_i -\beta X_i}{2\sigma^2}\Big) \] Función de verosimilitud:

\[ L(\sigma^2, \beta) = \prod^{N}_{i=1}f(Y_i|\sigma^2, \beta) = \prod^{N}_{i=1} \frac{1}{\sqrt{2\pi\sigma^2}}\exp\Big(-\frac{Y_i -\beta X_i}{2\sigma^2}\Big) \] \[ \log L(\sigma^2, \beta)=\sum^{n}_{i=1} -\frac{n}{2}\log(2\pi)-\frac{n}{2}\log(\sigma^2) -\frac{1}{2\sigma^2}(Y_i - \beta X_i)^2 \]

Paso 1: Generar los datos

library(pacman)
p_load(data.table, fixest, lattice, magrittr, ggplot2, kableExtra,dplyr)

## Installing package into '/home/juan/R/x86_64-redhat-linux-gnu-library/4.3'
## (as 'lib' is unspecified)

## 
## kableExtra installed

## Warning in p_load(data.table, fixest, lattice, magrittr, ggplot2, kableExtra, : Failed to install/load:
## kableExtra

N = 500      # Tamaño de la muestra
beta = 5     # Beta  
sigma_2 = 5  # sigma^2 (Distribución de U)
DT <- data.table(X = rnorm(N, 0, 5),
                 U = rnorm(N, 0, sqrt(sigma_2))) %>%
    .[, Y := beta*X + U]

Paso 2: Derivar la función de verosimilitud

# Función de verosimilitud
log_like <- function(theta, Y, X){   # theta en este caso es un vector con beta y sigma_2
  X <- as.matrix(X); Y <- as.matrix(Y);
  N       <- nrow(X)
  beta    <- theta[1]
  sigma_2 <- theta[2]
  e       <- Y - beta*X  # Cálculo de residuales
  loglik  <- -.5*N*log(2*pi)-.5*N*log(sigma_2) - ( (t(e) %*% e)/ (2*sigma_2) ) # Calcula log L, ver ecuación
  return(-loglik)  # Se utiliza el resultado negativo para la función optim 
}


# Graficar la función de verosimilitud

log_like_graph <- function(beta, sigma_2){
  X <- as.matrix(DT$X); Y <- as.matrix(DT$Y);
  N       <- nrow(X)
  e       <- Y - beta*X
  loglik  <- -.5*N*log(2*pi)-.5*N*log(sigma_2) - ( (t(e) %*% e)/ (2*sigma_2) )
  return(loglik)
}
log_like_graph <- Vectorize(log_like_graph)

# Crear matriz de valores para graficar beta y sigma
beta_vals <- seq(-10,10, by =1)
sigma2_vals <- seq(1,10, by =1)
log_vals <- outer(beta_vals, sigma2_vals, log_like_graph)

persp(beta_vals, sigma2_vals, log_vals, theta = 7, phi =8 , r= 500)

Paso 3: Estimación por máxima verosimilitud

MLE_estimates <- optim(fn=log_like,                   # Función de verosimilitud
                       par=c(1,1),                    # Estimación inicial
                       lower = c(-Inf, -Inf),         # Límite inferior de los parámetros
                       upper = c(Inf, Inf),           # Límite superior de los parámetros
                       hessian=TRUE,                  # Devuelve el Hessiano 
                       method = "L-BFGS-B",
                       # Entradas personalizadas
                       Y = DT$Y,                     
                       X = DT$X)

# Examinar estimaciones
MLE_par <- MLE_estimates$par
MLE_SE <- sqrt(diag(solve(MLE_estimates$hessian)))
MLE <- data.table(param = c("beta", "sigma_2"),
                  estimates = MLE_par,
                  sd = MLE_SE)

knitr::kable(MLE)

param	estimates	sd
beta	4.971642	0.0202098
sigma_2	4.951841	0.3131818

Paso 4: Graficar estimaciones

log_like_graph_beta <- function(beta){
  sigma_2 = MLE_par[2]
  X = as.matrix(DT$X)
  Y = as.matrix(DT$Y)
  N = nrow(X)
  e = Y - beta*X
  loglik  <- -.5*N*log(2*pi)-.5*N*log(sigma_2) - ( (t(e) %*% e)/ (2*sigma_2) )
  return(loglik) 
}

log_like_graph_sigma2 <- function(sigma_2){
  beta = MLE_par[1]
  X = as.matrix(DT$X)
  Y = as.matrix(DT$Y)
  N = nrow(X)
  e = Y - beta*X
  loglik  <- -.5*N*log(2*pi)-.5*N*log(sigma_2) - ( (t(e) %*% e)/ (2*sigma_2) )
  return(loglik) 
}

# Vectorizar
log_like_graph_beta <- Vectorize(log_like_graph_beta)
log_like_graph_sigma2 <- Vectorize(log_like_graph_sigma2)

# Graficar
beta <- ggplot(data = data.frame(beta = 0), mapping = aes(beta = beta)) + 
    stat_function(fun = log_like_graph_beta) +
    xlim(-40,40) + theme_bw() +xlab("beta") + ylab("log lik")

sigma2 <- ggplot(data = data.frame(sigma2 = 0), mapping = aes(sigma2 = sigma2)) + 
    stat_function(fun = log_like_graph_sigma2) +
    xlim(1,20) + theme_bw() +xlab("sigma2") + ylab("log lik")

beta

sigma2

Referencia: https://www.youtube.com/watch?v=w3drLH-DFpE&t=77s&ab_channel=CalebLikesR

Links de utilidad:

https://www.analyticsvidhya.com/blog/2018/07/introductory-guide-maximum-likelihood-estimation-case-study-r/ https://stats.stackexchange.com/questions/88697/sample-from-a-custom-continuous-distribution-in-r