Ejercicios 05 - Ejemplo de Máxima verosimilitud

Se realiza un ejemplo de estimación de máxima verosimilitud para la función de densidad de la distribución binomial. \[ p(x, θ)=θ^x (1 - θ)^{1-x} \] Entonces, su función de verosimilitud corresponde a: \[ L(x_i, θ) = \prod_{i=1}^{N}p_i(x_i, θ) = \prod_{i=1}^{N}θ^x_i (1 - θ)^{1-x_i} \] \[ log(L(x_i, θ)) = log(θ) \sum_{i=1}^{N} x_i + log(1-θ)(N-\sum_{i=1}^{N} x_i ) \]

Paso 1: Generar los datos

# Tamaño de la muestra
N = 200 
# Probabilidad 
Prob = 0.1   
rango = seq(1,200)
# Distribución binomial
distribucion = dbinom(rango, size = N,prob = Prob)

Paso 2: Derivar la función de verosimilitud

# Función de verosimilitud
log_like <- function(theta, n, x){
  return(-log(theta)*sum(x) - log(1-theta)*(n - sum(x)))
}
# Valores para la probabilidad
set.seed(1)
valores = runif(n = N, min = 0, max = 1)

plot(valores, log_like(valores, N, distribucion), main = "Función de log verosimilitud")

Paso 3: Estimación por máxima verosimilitud

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

# Examinar estimaciones
MLE_par <- MLE_estimates$par
MLE_SE <- sqrt(diag(solve(MLE_estimates$hessian)))


MLE <- data.table(param = c("Probabilidad"),
                  estimates = MLE_par,
                  sd = MLE_SE)

kable(MLE)

param	estimates	sd
Probabilidad	0.1	0.0537837

Paso 4: Graficar estimaciones

graph_prob <- function(x=valores, n=N){
  prob = MLE_par
  return(log(prob)*sum(x) + log(1-prob)*(n - sum(x))) 
}

# Vectorizar
graph_prob <- Vectorize(graph_prob)


# Graficar
p <- ggplot(data = data.frame(prob = 0), mapping = aes(prob = prob)) + 
    stat_function(fun = graph_prob) +
    xlim(0,20) + theme_bw() +xlab("Probabilidad") + ylab("log like")

p