Distribución elegida

La distribución que se elige que ya fue usada en la tarea anterior es la distribución Delaporte del paquete Delaporte, esta tiene definida su función por:

\[P(n) = ∑_{i=0}^{n} \frac{Γ(α+i)β^{i}λ^{n-i}e^{- λ}}{Γ(α)i!(1 + β)^{α + i}(n-i)!} \] Al generar la función de verosimilitud de esta (log(L(\(\theta\)))), se obtiene: \[ln(L(\alpha, \beta, \lambda; n_{j})) = ∑_{j=0}^{n}ln(∑_{i=0}^{n_{j}} \frac{Γ(α+i)β^{i}λ^{n_{j}-i}e^{- λ}}{Γ(α)i!(1 + β)^{α + i}(n_{j}-i)!}) \]

Paso 1: Se generan los datos

Importación de los paquetes

library(Delaporte)
library(ggplot2)
library(plot3D)
library(data.table)
library(knitr)

Se define la semilla y se generan los datos

set.seed(123)
N = 5000 # Cantidad de los datos
alpha = 5 # Se define el alpha
beta = 5 # Se define el beta
lambda = 5 # Se define el lambda
datos <- data.table(X = rdelap(N, alpha = alpha, beta = beta, lambda = lambda))

Paso 2: Se deriva la función de Verosimilitud

#Funcion de Verosimilitud 

log_likelihood <- function(theta, data) { # Theta es un vector con alpha, beta y lambda
  pmf <- ddelap(data, theta[1], theta[2], theta[3])
  log_pmf <- log(pmf)
  resultado = (-(sum(log_pmf)))
  return (resultado)
}

Graficos en donde se fijan los parametros de alpha, beta y lambda en 5 para cada uno. Esto se hace ya que no se pueden graficar en 4 dimensiones.

Grafico para valores de alpha de 1 a 10 y beta de 1 a 10, lambda se define en 5

log_likelihood_graf <- function(alpha, beta)
{
  data <- datos$X
  pmf <- ddelap(data, alpha, beta, 5)
  log_pmf <- log(pmf)
  resultado = (sum(log_pmf))
  return (resultado)
}
log_likelihood_graf <- Vectorize(log_likelihood_graf)
  
valores_alpha = seq(1, 10)
valores_beta = seq(1, 10)



valor <- outer(valores_alpha, valores_beta, log_likelihood_graf)


persp(valores_alpha, valores_beta, valor, theta = 7, phi =8 , r= 500)

Grafico para valores de alpha de 1 a 10 y lambda de 1 a 10, beta se define en 5

log_likelihood_graf <- function(alpha, lambda)
{
  data <- datos$X
  pmf <- ddelap(data, alpha, 5, lambda)
  log_pmf <- log(pmf)
  resultado = (sum(log_pmf))
  return (resultado)
}
log_likelihood_graf <- Vectorize(log_likelihood_graf)
  
valores_alpha = seq(1, 10)
valores_lambda = seq(1, 10)


valor <- outer(valores_alpha, valores_lambda, log_likelihood_graf)

persp(valores_alpha, valores_lambda, valor, theta = 7, phi =8 , r= 500)

Grafico para valores de beta de 1 a 10 y lamda de 1 a 10, alpha se define en 5

log_likelihood_graf <- function(beta, lambda)
{
  data <- datos$X
  pmf <- ddelap(data, 5, beta, lambda)
  log_pmf <- log(pmf)
  resultado = (sum(log_pmf))
  return (resultado)
}
log_likelihood_graf <- Vectorize(log_likelihood_graf)
  
valores_beta = seq(1, 10)
valores_lambda = seq(1, 10)


valor <- outer(valores_beta, valores_lambda, log_likelihood_graf)

persp(valores_beta, valores_lambda, valor, theta = 7, phi =8 , r= 500)

Paso 3: Estimación por máxima Verosimilitud

MLE_estimates <- optim(fn=log_likelihood,                   # Función de verosimilitud
                       par=c(1,1,1),                    # Estimación inicial
                       lower = c(-Inf, -Inf, -Inf),         # Límite inferior de los parámetros
                       upper = c(Inf, Inf, Inf),           # Límite superior de los parámetros
                       hessian=TRUE,                  # Devuelve el Hessiano 
                       method = "L-BFGS-B",
                       # Entradas personalizadas
                       data = datos$X)

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

kable(MLE)
param estimates sd
alpha 5.209515 0.5347753
beta 4.861958 0.2911534
lambda 4.528150 1.1555432

Conclusiones:

Se puede notar que gracias a este método se puede llegar a los valores reales de lo que sería la población, sin embargo, para esta distribución no es del todo útil ya que la cantidad de datos que se necesitan son bastantes altos. Se hicieron pruebas para un N igual a 500 y los resultados fueron bastante alejados de los parámetros fijos que se utilizan en el comienzo y para valores de N 1000 y 5000 se fue acercando mucho más a los valores reales teniendo unos pocos decimales de diferencia, para el primer caso tiene varios enteros de diferencia.

Graficos de estimaciones

log_like_graph_alpha <- function(alpha) {
  data <- datos$X
  pmf <- ddelap(data, alpha, MLE_par[2], MLE_par[3])
  log_pmf <- log(pmf)
  resultado = (sum(log_pmf))
  return (resultado)
}
log_like_graph_alpha <- Vectorize(log_like_graph_alpha)


log_like_graph_beta <- function(beta) {
  data <- datos$X
  pmf <- ddelap(data, MLE_par[1], beta, MLE_par[3])
  log_pmf <- log(pmf)
  resultado = (sum(log_pmf))
  return (resultado)
}
log_like_graph_beta <- Vectorize(log_like_graph_beta)

log_like_graph_lambda <- function(lambda)
{
  data <- datos$X
  pmf <- ddelap(data, MLE_par[1], MLE_par[2], lambda)
  log_pmf <- log(pmf)
  resultado = (sum(log_pmf))
  return (resultado)
}
log_like_graph_lambda <- Vectorize(log_like_graph_lambda)

Gráfica de alpha

alpha <- ggplot(data = data.frame(alpha = 0), mapping = aes(alpha = alpha)) + 
    stat_function(fun = log_like_graph_alpha) +
    xlim(1,10) + theme_bw() +xlab("alpha") + ylab("log lik")

plot(alpha)

Gráfica de beta

beta <- ggplot(data = data.frame(beta = 0), mapping = aes(beta = beta)) + 
    stat_function(fun = log_like_graph_beta) +
    xlim(1,10) + theme_bw() +xlab("beta") + ylab("log lik")

plot(beta)

Gráfica de lambda

lambda <- ggplot(data = data.frame(lambda = 0), mapping = aes(lambda = lambda)) + 
    stat_function(fun = log_like_graph_lambda) +
    xlim(1,10) + theme_bw() +xlab("lambda") + ylab("log lik")

plot(lambda)