La función de distribución exponencial se define como:
\(\lambda^{n}e^{-\lambda\sum_{i = 1}^{n}x_{j}}\)

Mientras que su log-likelihood function se define como:
\(n*ln(\lambda)-\lambda\sum_{i = 1}^{n}x_{j}\)
= \(ln(\lambda^{n}e^{-\lambda\sum_{i = 1}^{n}x_{j}})\)
= \(ln(\lambda^{n})+ln(e^{-\lambda\sum_{i = 1}^{n}x_{j}})\)
= \(n*ln(\lambda)-\lambda\sum_{i = 1}^{n}x_{j}\)

Paso 1: Generar los datos

Para poder calcular la función de máxima verosimilitud, se especifica el valor que se desea estimar y el valor n de la función
N = 25      #Tamaño de la muestra
n = 200     #Valor de n de la funcion log-likelihood
lambda = 5  #Valor que se quiere aproximar
datos = rexp(n=N, rate=lambda)  #Obtiene n observaciones de una distribucion exponencial

Paso 2: Derivar la función de verosimilitud

Ahora, para poder derivar la función de verosimilitud, se crea una función que la calcule junto a una otra que entrega el mismo valor pero negativo. Además se grafica la primera función dado un intervalo de valores lambda.
exp_log_likelihood = function(lambda, n, x){
  return (n * log(lambda) - lambda*sum(x))
}
#Copia de la funcion anterior pero retorna negativo
exp_copia = function(lambda, n, x){
  return ((n * log(lambda) - lambda*sum(x))*(-1))
}

#Gráfico
intervalo_lambda = seq(1,N)
plot(intervalo_lambda, exp_log_likelihood(intervalo_lambda, n, datos), type='l', main = 'Función de Verosimilitud', ylab = 'Exponencial log-likelihood', xlab = 'lambda')

Paso 3: Estimación por máxima verosimilitud

La función que retorna el valor negativo, se utiliza para calcular el óptimo que se presenta a continuación
datos = rexp(n=N, rate=lambda)
MLE_estimates <- optim(fn=exp_copia,                # Función de verosimilitud
                       par=c(1),                    # Estimación  
                       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",
                       n = N,                    
                       x = datos)                   

MLE_par <- MLE_estimates$par 
MLE_SE <- sqrt(diag(solve(MLE_estimates$hessian)))
# Tabula el valor calculado estimado y la desviacion estandar
MLE <- data.table(param = c("Lambda"),
                  estimates = MLE_par,
                  sd = MLE_SE)
kable(MLE)
param estimates sd
Lambda 5.483088 1.096618
Al estimar el valor de lambda, se muestra una tabla que contiene el valor estimado calculado y su desviación estandar

Paso 4: Graficar estimacion

Finalmente, se grafica la estimación obtenida
log_like_graph_beta <- function(x = datos, N = n){
  lambda = MLE_par
  loglik = N * log(lambda) - lambda*sum(x)
  return(loglik) 
}

# Vectorizar
log_like_graph_beta <- Vectorize(log_like_graph_beta)

# Graficar
lambda = ggplot(data = data.frame(lambda = 0), mapping = aes(lambda = lambda)) + 
    stat_function(fun = log_like_graph_beta) +
    xlim(-40,40) + theme_bw() +xlab("lambda") + ylab("log lik")

lambda