Ejercicio4

Ejercicios 04 - Ejemplo de Máxima verosimilitud

Estadística Computacional

José Carlos Marín García

Introducción

En el presente se estimará un parámetro de la Distribución Exponencial mediante el método de la máxima verosimilitud.

Distribución Exponencial

\[f(x)=\lambda*e^{-\lambda*x}\]

Función de Máxima verosimilitud

\[L(\lambda)=\lambda^{n}e(-\lambda\sum^n_{j=1}x_j)\]

Forma logarítmica

\[log(L(\lambda))=nln(\lambda)-\lambda \sum^n_{j=1}f(x_j)\]

Paso 1: Generar datos

N = 50           # Tamaño de la muestra
lambda = 3       # Valor de lambda, se puede cambiar para obtener distintas aproximaciones
distribucion = rexp(N,lambda) #Generar datos por la distribución exponencial

Al usar lambda igual a 3 significa que el resultado deberia darnos una aproximación de p cercana a 3.

Paso 2: Derivar la función de verosimilitud

#Función de verosimilitud en forma logarítmica
exp_log_likelihood = function(lambda, x, n){
  return(n*log(lambda)-lambda*sum(x))
}

#Valores de lambda para graficar
exp_values = seq(1,N)

#Graficar la función de verosimilitud
plot(exp_values, exp_log_likelihood(lambda = exp_values, x = distribucion, n = 200), main = 'Función de Verosimilitud con logaritmos', xlab = 'Valores de lambda', ylab = 'log-likelihood de lambda')

Paso 3: Estimación por máxima verosimilitud

exp_log_likelihood_neg <- function(lambda, x, n){
  return (-(n*log(lambda)-lambda*sum(x))) #Funcion de verosimilitud en su forma logarítmica multiplicada por -1 para optimizar el máximo con la función optime
}

MLE_estimates <- optim(fn=exp_log_likelihood_neg,  # Función de verosimilitud multiplicada por -1
                       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
                       n = N,                     
                       x = distribucion)

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

cat("- Valor Original de Lambda:",lambda)

## - Valor Original de Lambda: 3

cat("\n- Valor Estimado de Lambda:",MLE_par)

## 
## - Valor Estimado de Lambda: 2.89663

cat("\n- Desviación estandar:", MLE_SE)

## 
## - Desviación estandar: 0.4096453

Paso 4: Graficar Estimaciones

log_like_graph <- function(x = exp_values, n = N){ 
  lambda = MLE_par # valor de lamdba encontrado
  log_like_g <- (n*log(lambda)-lambda*sum(x)) # Función de máxima verosimilitud forma logarítmica
  return (log_like_g)
}
# Vectorizar
log_like_graph <- Vectorize(log_like_graph) 
# Graficar
lambda_graph <- ggplot(data = data.frame(lambda = 0), mapping = aes(lambda = lambda)) + stat_function(fun = log_like_graph) + xlim(0,N) + theme_bw() +xlab("lambda_vals") + ylab('log lik_vals')
lambda_graph