Resumen

A continuación se encuentran las respuestas y su código correspondiente a los enunciados presentados en la quinta sesión de Ejercicio del Ramo de Estadística Computacional.

Paso 1: Generar los datos

library(pacman)
p_load(data.table, ggplot2)

Utilizando una función de probabilidad de distribución de Poisson

x = 50 #Tamaño muestra
lambda = 3 #Valor esperado
datos <- data.table(rpois(n = x, lambda = lambda)) #Dist. de Poisson

Paso 2: Derivar la función de verosimilitud

La función de verosimilitud es: \[ -n\lambda - \sum_{i=1}^{n} ln(x_i !) + ln(\lambda)\sum_{i=1}^{n} x_i \] Fuente: 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.

#Función de verosimilitud

funcion_ver = function(x, lambda, n){
  func_ver = -n*lambda - sum(log(factorial(x))) + log(lambda)*sum(x)
  return(-func_ver) #Valor negativo de la función para la función optim
}

#Graficar la función de verosimilitud

Lambdas <- seq(0,10, by =1) #Valores para graficar función
Funcion = funcion_ver(lambda=Lambdas, x=datos, n=x)
df <- data.frame(Lambdas, Funcion)
ggplot(data = df, aes(x=Lambdas, y=Funcion)) + geom_point() + xlab("Valores de Lambda") + ylab("Función de verosimilitud")

Paso 3: Estimación por máxima verosimilitud

MLE_estimates <- optim(fn=funcion_ver,                # Función de verosimilitud
                       par=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
                       x = datos,
                       n = x)

# Examinar estimaciones
MLE_par <- MLE_estimates$par
MLE_SE <- sqrt(diag(solve(MLE_estimates$hessian)))
MLE <- data.table(param = "lambda",
                  estimates = MLE_par,
                  sd = MLE_SE)

knitr::kable(MLE)
param estimates sd
lambda 3.2 0.2529822