Paso 0: Parte teórica
Si tenemos la función de densidad de una distribución exponencial
\[ f(x;\theta)=\theta e^{-\theta x}
\]
Además de una muestra independiente e identicamente distrubuida como la siguiente:
\[x_1,x_2,...,x_n \]
Podemos construir la siguiente función de verosimilitud:
\[ L(\theta)=\prod_{i=1}^{n}\theta
e^{-\theta x_i} \]
Lo cual se puede simplificar a una función con una sola sumatoria de la siguiente forma:
\[ L(\theta)=(\theta )^{n}e^{-\sum_{i=1}^n
\theta * x_i} \]
Para simplificarlo más podemos aplicar logaritmo quedandonos:
\[ ln(L(\theta))=n*ln(\theta ) - \sum_{i=1}^n \theta * x_1 \]
Derivando nos queda:
\[
\frac{dln(\theta)}{d\theta}=\frac{n}{\theta}- \sum_{i=1}^n x_1 \]
Que igualado a 0 para maximizar nos queda:
\[ \theta = \frac{n}{\sum_{i=1}^n x_1}
\]
Reordenamos la fracción de manera conveniente quedando:
\[ \theta = \frac{1}{\frac{\sum_{i=1}^n
x_1}{n}} \]
Pero sabemos que la sumatoria de los elementos de una muestra, dividida en el número de elementos de la misma, es el promedio; debido a todo lo anterior el estimador de máxima verosimilitud será:
\[ \hat{\theta}_{MV} = \frac{1}{\bar{x}}
\]
Ya habiendo obtenido lo anterior, volvamos a lo que nos compete, la estimación del mismo en r.
Paso 1: Generación de datos
N = 1500 # Tamaño de la muestra
theta = 5
DT <- data.table (X = rexp(N,16000))Paso 2: Derivar función de verosimilitud
\[ ln(L(\theta))=n*ln(\theta ) - \sum_{i=1}^n \theta * x_1 \]
ll <- function(theta,X){
N <- length(X)
-1 * N * log(theta) + (theta * sum(X)) # Notar que la función de ln(L(theta)) se multiplico por -1 para el posterior uso de optim para buscar un máximo.
}
log_like_graph <- function(theta){
X <- as.matrix(DT$X)
N <- nrow(X)
N * log(theta) + theta * sum(X)
}
log_like_graph <- Vectorize(log_like_graph)
theta_vals <- seq(0,1,by=0.01)
log_vals <- log_like_graph(theta_vals)
plot(theta_vals, log_vals)
Paso 3: Estimación por máxima verosimilitud
MLE_estimates <- optim(fn=ll, # Función a minimizar (con el -1 puesto previamente estariamos maximizando la real)
par = 0, # Valor inicial para theta
lower = 0.000000000001, # Valor mínimo para theta
upper = Inf, # Valor máximo para theta
method = "L-BFGS-B",
hessian = TRUE,
X = DT$X)
MLE_par <- MLE_estimates$par
MLE_SE <- sqrt(diag(solve(MLE_estimates$hessian)))
MLE <- data.table(theta = c("theta"),
estimates = MLE_par,
sd = MLE_SE)
kable(MLE)| theta | estimates | sd |
|---|---|---|
| theta | 15744.55 | 468.9374 |
Paso 4: Graficar estimaciones
log_like_graph_theta <- function(theta){
theta = MLE_par
X <- as.matrix(DT$X)
N <- nrow(X)
N * log(theta) + theta * sum(X)
}
# Vectorizar
log_like_graph_theta <- Vectorize(log_like_graph_theta)
# Gráficar
theta <- ggplot(data = data.frame(theta=0), mapping = aes(theta = theta)) +
stat_function(fun = log_like_graph_theta) + xlim(-1,1) + theme_bw() + xlab("theta") + ylab("log_like")
print(theta)