La función de densidad de la distribución de Poisson está dada por la siguiente formula: \(f(x)=\frac{e^{-\lambda}\lambda^{x}}{x!},\)\(x>0\)

Usando el logaritmo de la productoria de la función se genera la siguiente función de Verosimilitud:

\(L(\lambda) = \prod_{i=1}^{N} \frac{e^{-\lambda}\lambda^{x}}{x!}\\\)
\(\log(L(\lambda)) = \sum_{i=1}^{n}log(\frac{e^{-\lambda}\lambda^{x}}{x!})\\\)
\(=\sum_{i=1}^{n} (-\lambda + x\log(\lambda) - \log(x!))\\\)
\(= -n\lambda + \log(\lambda) \sum_{i=1}^{n} x \space -\sum_{i=1}^{n}\log(x!)\\\)

Paso 1: Generar los datos

Se define el tamaño de la muestra y un valor arbitrario de \(\lambda\), el cual representará la esperanza de repeticiones del evento de un intervalo dado:

N = 500                        # Número de observaciones   
lambda = 4                     # Lambda arbitrario, se puede cambiar
data = rpois(N, lambda)        # Se obtienen los datos con rpois para el N y lambda ingresado
D_table = data.table(X = data) # Se almacenan los datos obtenidos en un data table

Paso 2: Derivar la función de Verosimilitud

Hacemos uso del logaritmo de la función de Verosimilitud para realziar el calculo de los datos

Vero_funct = function(n, lambda, x){
  a = sum(log(factorial(x)))
  b = log(lambda)*sum(x)
  c = -n*lambda
  return(-(c+b-a))
}
# Se genera el gráfico de la función
elements = seq(1,60)
plot(elements,Vero_funct(N,elements,data) , xlab = "Lambda", ylab = "Log(L(λ))", type='o', 
     main = "Funcion de Verosimilitud Poisson con λ = 4", sub = "[1,60]", col=" darkgreen")

Paso 3: Estimación de Máxima Verosimilitud

Se hace uso de la función optim() para optimizar la función. Esto nos retornará la estimación del valor de λ

MLE_stimates = optim(fn=Vero_funct,      # Función de verosimilitud
                     par=c(1),           # Estimación inicial
                     lower = -Inf,       # Límite inferior de los parámetros
                     upper = Inf,        # Límite superior de los parámetros
                     hessian=TRUE,       # Hessiano
                     method= "L-BFGS-B",
                     # Entradas personalizadas
                     n = N,              # Parámetro n
                     x = data)           # Parámetro x
# Examinar los resultados
MLE_par = MLE_stimates$par
MLE_SE = sqrt(diag(solve(MLE_stimates$hessian)))
cat("El valor original de Lambda era: ", lambda)
## El valor original de Lambda era:  4
cat("La estimación de Lambda calculada con la Máxima Verosimilitud es: ", MLE_par)
## La estimación de Lambda calculada con la Máxima Verosimilitud es:  3.892

El valor obtenido no está lejos del indicado, esta diferencia ocurre debido a que rpois() posee una aleatoriedad al generar datos. Sin embargo, siempre se dará una aproximación cercana al valor real de λ.

Paso 4: Graficar Estimación

lambda_graph = function(x = data, n = N){
  lambda = MLE_par
  a = sum(log(factorial(x)))
  b = log(lambda)*sum(x)
  c = -n*lambda
  return(c+b-a)
}
lambda_graph = Vectorize(lambda_graph)

graph <- ggplot(data = data.frame(lambda = 0), mapping = aes(lambda = lambda)) +
  stat_function(fun = lambda_graph) + 
  xlim(0,60) + theme_bw() + xlab("λ") + ylab("Log(L(λ))")

graph

Se da entonces por finalizado el proceso de lograr calcular una estimación de Máxima Verosimilitud para la distribución de Poisson. Se logró encontrar una aproximación muy cercana al valor real del lambda escogido.