Utilizando la función de la distribución de poisson: \[\frac{\lambda^xe^{-\lambda}}{x!}\] Generamos la productoria (Likelihood Function): \[\Pi_{i=1}^{n}\frac{\lambda^{x_i}e^{-\lambda}}{x_i!}\]

\[\log (\Pi_{i=1}^{n}\frac{\lambda^{x_i}e^{-\lambda}}{x_i!}) = \sum_{i=1}^{n}\log(\frac{\lambda^{x_i}e^{-\lambda}}{x_i!}) = \sum_{i=1}^n x_i\log(\lambda) - \sum_{i=1}^n \lambda - \sum_{i=1}^n \log (x_i!)\] \[=\ - n\lambda + \sum_{i=1}^nx_ilog(\lambda) - \sum_{i=1}^n \log (x_i!)\] Lo anterior corresponde a la Función de Verosimilitud. Ahora, podemos realizar los pasos para obtener la estimación de Lambda, basandonos en el enunciado.

PASO 1: Generar los datos

Definimos las librerias y datos a utilizar.

library(pacman)
library(knitr)
p_load(data.table, fixest, lattice, magrittr, ggplot2, kableExtra,dplyr)
N = 200 # Tamano de la muestra
lambda = 6 # ESTE VALOR SE PUEDE CAMBIAR PARA OBTENER DISTINTAS APROXIMACIONES
poison_data = rpois(n = N, lambda =lambda) # Generamos una distribucion de poisson
data_frame_poison = data.frame(poison_data) # Convertimos lo anterior en data frame

Vamos a utilizar Lambda = 6, lo cual significa que el resultado deberia darnos una aproximación de lambda cercana a 6. Si se desea, el valor de lambda se puede cambiar para obtener distintos resultados y aproximaciones.

PASO 2: Derivar la función de verosimilitud

Definimos la función anteriormente obtenida

poisson_log_likelihood = function(lambda, x, n){ # Definiendo la funcion de verosimilitud con logaritmos
  return(-(-n * lambda + sum(x) * log(lambda) - sum(log(factorial(x)))))
}

lambda_values = seq(1,20) # Definimos distintos valores de lambda, para hacer la grafica

plot(lambda_values, poisson_log_likelihood(lambda = lambda_values, x = poison_data, n = 200), type='o', main = 'Función de Verosimilitud con logaritmos, considerando lambda = 6', ylab = 'log-likelihood de lambda', col = "blue")

# Lo anterior se puede transformar a un dataframe para visualizar la información.
df = data.frame(lambda_values, poisson_log_likelihood(lambda_values, poison_data, n= 200))
df
##    lambda_values poisson_log_likelihood.lambda_values..poison_data..n...200.
## 1              1                                                   1592.4851
## 2              2                                                    959.3222
## 3              3                                                    671.9531
## 4              4                                                    526.1593
## 5              5                                                    457.9407
## 6              6                                                    438.7902
## 7              7                                                    453.5011
## 8              8                                                    492.9963
## 9              9                                                    551.4211
## 10            10                                                    624.7778
## 11            11                                                    710.2150
## 12            12                                                    805.6273
## 13            13                                                    909.4160
## 14            14                                                   1020.3382
## 15            15                                                   1137.4087
## 16            16                                                   1259.8334
## 17            17                                                   1386.9626
## 18            18                                                   1518.2582
## 19            19                                                   1653.2694
## 20            20                                                   1791.6149

Notar que la función se retorna con negativo. Esto es para poder optimizarla correctamente en el siguiente paso.

PASO 3: Estimación por máxima verosimilitud

Utilizamos optim para optimizar la funcion. Esto nos retorna la estimación del valor de lambda y el valor de sd.

MLE_estimates <- optim(fn=poisson_log_likelihood,   # Función de verosimilitud
                       par=c(1),                    
                       lower = c(-Inf, -Inf),       # Límite inferior de los parámetros
                       upper = c(Inf, Inf),         # Límite superior de los parámetros
                       hessian=TRUE,                # Hessiano
                       method = "L-BFGS-B",
                       n = 200,                     # Parametro n
                       x = poison_data)             # Parametro x

MLE_par <- MLE_estimates$par # Valor de lambda estimado
MLE_SE <- sqrt(diag(solve(MLE_estimates$hessian))) # Valor de sd
cat("El valor original de Lambda era: Lambda = ", lambda)
## El valor original de Lambda era: Lambda =  6
cat("La estimación de Lambda utilizando Maxima Verosimilitud es: Lambda = ", MLE_par)
## La estimación de Lambda utilizando Maxima Verosimilitud es: Lambda =  6.01
cat("El sd asociado a la estimación es: sd = ", MLE_SE)
## El sd asociado a la estimación es: sd =  0.1733494

Notar que debido a que rpois() genera información diferente cada vez que se ejecuta, al ejecutar todo el código el valor de la estimación variará. Sin embargo, siempre dará una aproximación cercana al valor real de lambda

PASO 4: Graficar estimaciones

Utilizamos toda la información anterior para graficar la función de verosimilitud con el valor de lambda encontrado

log_like_graph <- function(x = poison_data, n = N){ # Misma funcion de antes
  lambda = MLE_par # Utilizamos el valor de lamdba encontrado
  log_like_g <- -n * lambda + sum(x) * log(lambda) - sum(log(factorial(x)))
  return (log_like_g)
}
log_like_graph <- Vectorize(log_like_graph) # Vectorizamos la funcion para poder graficarla

# Graficamos la función
lambda_graph <- ggplot(data = data.frame(lambda = 0), mapping = aes(lambda = lambda)) + 
    stat_function(fun = log_like_graph) +
    xlim(0,40) + theme_bw() +xlab("lambda") + ylab("log lik")

lambda_graph

Con esto, se concluye el proceso de encontrar una estimación de maxima verosimilitud para la distribución de poisson. Finalmente, fue posible encontrar una aproximación muy cercana al valor real de lambda escogido.

Referencias utilizadas:

1. Taboga, Marco (2021). “Poisson distribution - Maximum Likelihood Estimation”, Lectures on probability theory and mathematical statistics. Kindle Direct Publishing. Online appendix. https://www.statlect.com/fundamentals-of-statistics/Poisson-distribution-maximum-likelihood.

2. Singla, AanishS (2018). “An Introductory Guide to Maximum Likelihood Estimation (with a case study in R)”. Lectures on probability theory and mathematical statistics. https://www.analyticsvidhya.com/blog/2018/07/introductory-guide-maximum-likelihood-estimation-case-study-r/

3. Zach (2020). “MLE for a Poisson Distribution (Step-by-Step)”. Lectures on probability theory and mathematical statistics. https://www.statology.org/mle-poisson-distribution/