La funcion de densidad de la distribucion de Poisson esta dada por la siguiente formula:

\[ f(x) = \frac{e^{-\lambda}\lambda^{x}}{x!} \space\space,\space\space x>0 \] Utilizando el logaritmo de la productoria de la funcion podemos obtener la siguiente funcion de logaritmo de la funcion 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!}) \\ \space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space =\sum_{i=1}^{n} (-\lambda + x\log(\lambda) - \log(x!)) \\ \space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space = -n\lambda + \log(\lambda) \sum_{i=1}^{n} x \space -\sum_{i=1}^{n} \log(x!) \]

Paso 1: Generacion de Datos

Se generaran datos con un valor arbitrario de λ, el cual representa esperanza de repeticiones del evento en el intervalo dado:

N = 500      # Tamaño de la muestra
lambda = 4   #lambda arbitrario
data <- rpois(N,lambda)
DT <- data.table(X = data)

Paso 2: Derivar la Función de Verosimilitud

Se utiliza el logaritmo de la funcion de verosimilitud para simplificar el calculo de los datos.

#Definicion de la Funcion de Verosimilitud para Poisson
verosimilitud <- function(n,lambda,x){
  a <- -n*lambda
  b <- log(lambda)*sum(x)
  c <- sum(log(factorial(x)))
  return(-(a+b-c))
}
#Graficando la Funcion
lambdas <- seq(1,50)
plot(lambdas,verosimilitud(N,lambdas,data) , xlab = "Lambda", ylab = "log(L(lambda))" , main = "Funcion de Verosimilitud Poisson", sub = "[1,50]")

Paso 3: Estimación por Maxima Verosimilitud

Ahora utilizamos la funcion optim() para comprobar el valor de λ para nuestra muestra

MLE_estimates <- optim(fn=verosimilitud,                   # Función de verosimilitud
                       par=c(1),                    # Estimación inicial
                       lower = c(0,0),         # 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
                       x = data,
                       n = N)

# Examinar estimaciones
MLE_par <- MLE_estimates$par
MLE_SE <- sqrt(diag(solve(MLE_estimates$hessian)))
## El valor estimado para lambda es de 4.004 ± 0.08948742 . El valor indicado para la generacion de la muestra fue un lambda de 4

Este valor no esta lejos del indicado. La diferencia con el valor inicial de la generacion de la muestra se debe a la aleatoriedad al generar datos, ya que estos tienen una tendencia a mantenerse cerca del valor que le dimos a la funcion rpois().

Paso 4: Graficar Estimación

graphlambda <- function(x){
  lambda <- MLE_par
  a      <- -N*lambda
  b      <- log(lambda)*sum(x)
  c      <- sum(log(factorial(x)))
  return(a+b-c)
}
graphlambda <- Vectorize(graphlambda)

glambda <- ggplot(data = data.frame(lambda = 0), mapping = aes(lambda = lambda)) + 
    stat_function(fun = graphlambda) +
    xlim(0,50) + theme_bw() +xlab("lambda") + ylab("log(L(lambda))")

glambda