MLE basics with R

Inferencia Estadística 1
Author

Nelson Muriel

Abstract
En este documento aprenderemos los básicos de la estimación por máxima verosimilitud usando R.

Máxima verosimilitud

Recordemos que el EMV se puede obtener usando

\[ \hat{\theta} = \text{armax}_{\theta} \mathcal{L}(\theta\vert X) = \text{argmax}_\theta \log L(\theta\vert X) \]

La maximización puede hacerse con la función optimize: ?optimize para funciones de una variable.

Ejemplo: Caso Poisson

Supongamos que las siguientes son observaciones del número de fallas microscópicas en el filo de los cuchillos de una marca japonesa:

x_obs <- c(41, 30, 31, 38, 29, 24, 30, 29, 31, 38)

Modelamos estas observaciones como \(X\sim Poi(\lambda)\). La log-verosimilitud para los datos sería:

\[ \log(\lambda) * \sum_i x_i - n\lambda - \sum_i\log(x_i !) \] que podemos programar en R en dos pasos

poisson_density <- function(lambda, x){
  n <- length(x)
1  log(lambda) * sum(x) + n * lambda - sum(lfactorial(x))
}
1
Obsérvese que R tiene una lógica vectorial como en lfactorial(x). Además, R retorna el último objeto calculado sin la necesidad de escribir return(...).

Esta función no es la log-verosimilitud, puesto que ésta última evalúa en \(\lambda\). Tendríamos que proyectar la muestra, como sigue

log_likelihood_poisson <- function(x) {
  n <- length(x)

  function(lambda) {
1    log(lambda) * sum(x) - n * lambda - sum(lfactorial(x))
  }
}
1
En esta línea se calcula la función de verosimilitud (función de \(\lambda\)) tomando la muestra como información dada.

Con esta sintáxis, la log-verosimilitud de nuestra muestra es la función log_likelihood_poisson(x_obs), y la podemos evaluar en distintos valores de \(\lambda\). Por ejemplo

sapply(c(10, 20, 30), log_likelihood_poisson(x_obs))
[1] -183.64052  -61.14028  -30.98598

Para encontrar el MLE haríamos

optimise(log_likelihood_poisson(x_obs), c(0, 100), maximum = TRUE)
$maximum
[1] 32.09999

$objective
[1] -30.26755

En R tenemos acceso a distintos modelos, y no hay necesidad de programarlos. Consulte ?Distribuions. En todos los casos, la proyección a la muestra concreta es indispensable.

log_likelihood_poisson_base_R <- function(x){
  function(lambda){
    sum(dpois(x, lambda, log = TRUE))
  }
}

optimise(log_likelihood_poisson_base_R(x_obs), c(0, 100), maximum=TRUE)
$maximum
[1] 32.09999

$objective
[1] -30.26755

Ejercicio

Considere el siguiente vector, que puede copiar directamente

x_obs <- c(0.0412310130174156,0.195189960292521,0.453766801555581,0.307077612684749,0.661722530390988,1.00422792168041,0.0727365856990218,0.403335356695538,0.153289198680761,0.828314851520833,0.221012709467652,0.442671183844037,0.391008945537691,0.0660335287451744,0.0477391151834358,0.235628273143413,0.0537888980470598,0.0943330381065607,0.0290931310697871,0.0266665378585458,0.183266733019593,0.173173571928384,0.120138790644705,0.00679763713851571,0.0537535726092756,0.0160877287481951,0.169263637057486,0.14157001878226,0.093747469317168,0.241380632287803,0.0602334664203227,0.308564917774499,0.642339067295299,0.160926070428433,0.231607764019611,0.717851034882379,0.310183669004211,0.22956758914354,0.649262210693131,0.105220932979137,0.379523840157275,0.354904533905869,0.135148333851248,0.54571561401052,0.132814994081855,0.308316260931974,0.16144867389554,0.11188332517349,0.044416388055724,0.0105188742434377,0.23231549566956,0.00407959753647447,0.0686375360935926,0.0826229440048337,0.172885744265261,0.18565986882055,0.431395378311105,0.158111213325503,0.0599826316349208,0.0201660940838673,0.0145562116988003,0.00239994266847628,0.040844739083439,0.894316438535109,0.180753523411993,0.0877043782733381,0.355290122685914,0.043228092894823,0.0297964914126924,0.270764632142129,0.132236272841692,0.00819995263591409,0.197206279999557,0.0638536755970136,0.150817323095145,0.0651220780797303,0.0119542467407882,0.0888681131415069,0.218481621584612,0.120457602106035,0.179358789502623,0.0595679681748152,0.0304259879448669,0.306231372041588,0.254977236493471,0.018758587911725,0.0900713274255395,0.241971019783076,0.0951698533259332,0.0875889146540307,0.488716150422521,0.426314349437218,0.125713923480362,0.0444632913917303,0.0282174453139305,0.109752067085356,0.171521576113048,0.304057890159969,0.413298829628115,0.0230537567287683)

  1. Genere un gráfico de histograma o de densidad usando ggplot2 para elegir una distribución posible para estos datos. Elija una que tenga un único parámetro

  2. Programe la proyección al parámetro de la función de verosimilitud para el modelo elegido

  3. Realice la estimación por MLE

Ejercicio

Consulte la función optim para funciones vectoriales. Ponga a prueba su comprensión haciendo lo siguiente:

  1. Simule una muestra gaussiana con n = 200 observaciones. Elija usted su media \(mu\) y su desviación estándar \(\sigma\).

  2. Proyecte la densida muestral para obtener la verosimilitud

  3. Optimice la verosimilitud y compare los resultados obtenidos con los parámetros que eligió

Ejercicio

Repitamos lo anterior generando un estudio de simulación que nos ayude a entender cómo se comportan los estimadores de máxima verosimilitud cuando cambiamos el tamaño de muestra. Consideraremos un mismo modelo, el \(N(5, 2)\), y simularemos muestras de tamaño \(n = 50, 100, 200, 500\). Repetiremos la simulación y la estimación para evaluar la variabilidad muestral de los estimadores.

library(latex2exp)
library(tidyverse)

negloglik_normal <- function(x){
   function(theta){
1    -sum(dnorm(x, mean = theta[1], sd = theta[2], log=TRUE))
   }
}

estimar_parametros_normales <- function(n_obs, mu, sigma){
   x_sample <- rnorm(n_obs, mean = mu, sd = sigma)
   objective_function <- negloglik_normal(x_sample)
2   optim_results <- optim(c(0, 1), objective_function)
   optim_results$par |>
        t() |> 
        as.data.frame() |> 
        set_names('mu', 'sigma') |> 
3        mutate(n = n_obs)
}

estimar_parametros_all_samples <- function(sample_sizes, mu, sigma){
  sample_sizes |> 
    purrr:::map(estimar_parametros_normales, mu = mu, sigma = sigma) |> 
4    dplyr::bind_rows()
}

desired_sample_sizes <- c(50, 100, 200, 500)

5all_estimates <- purrr::rerun(1000,
                        estimar_parametros_all_samples(desired_sample_sizes, 5, 2)) |> 
                        dplyr::bind_rows()

all_estimates |> 
      ggplot(aes(x = mu, y = sigma)) +
      geom_point(alpha=0.5) +
      geom_density_2d() +
      facet_wrap(~n, scales = "free") +
      labs(title = 'Contornos de la densidad empírica bivariada en la estimación por máxima verosimilitud',
           subtitle = 'Muestras Gaussianas con media 5.0 y desviación estándar de 2',
           caption = 'Resultados empíricos con 1000 simulaciones') +
      xlab(TeX('$\\hat{\\mu}_{MLE}$')) +
      ylab(TeX('$\\hat{\\sigma}_{MLE}$'))
1
La log–verosimilitud es negativa porque estamos usando la función optim sin lista de control. Podría usarse la log–verosimilitud positiva y añadir el argumento control = list(fnscale = -1).
2
En la función estimar_parametros_normales incluimos la simulación de n_obs observaciones de la Normal con parámetros mu y sigma y obetnemos los resultados crudos de la optimización en el objeto optim_results
3
Estos resultados crudos se transforman en una base de datos a la cual se añade el tamaño de muestra, información útil para la generación del gráfico
4
La función map del paquete purrr simplifica la escritura del ciclo. Esta función es equivalente a un ciclo for que automáticamente guarda los resultados en una lista. La función bind_rows del paquete dplyr junta esta lista de bases de datos en una base de datos global
5
La función rerun del paquete purrr simplifica la repetición de la ejecución de estimar_parametros_all_samples y agrupa los resultados con bind_rows.