Ejercicios 05 - Ejemplo de Máxima verosimilitud

logo

Instrucciones

El trabajo debe ser desarrollado en el lenguaje de programación R (https://cran.r-project.org/).
Trabaje individualmente.
Cree un archivo en Rpubs para documentar cada respuesta y procedimiento asociado a la resolución de este ejercicio.
La dirección (url) respectiva deberá ser subida en la sección de Moodle habilitada para este fin.
Es suficiente con que sólo un o una integrante suba la dirección a Moodle.

Ejemplo de máxima verosimilitud

Se solicita crear un ejemplo de estimación de máxima verosimilitud en el que explique cada paso. Puede utilizar como base el siguiente código utilizando la siguiente función de densidad, suponga que \(\delta\) es conocido: \[ \begin{align} f(x;\beta,\delta)=\beta \delta x^{\delta-1}e^{-\beta x^\delta} && \text{para x>0, 0 en otro caso} \end{align} \]

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)

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)

kable(MLE)

param	estimates	sd
beta	4.963228	0.0192996
sigma_2	4.617199	0.2920173

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