Código
<- rexp(n = 100, rate = 1)
x <- length(x)
n <- sum(x)
y <- seq(0.5, 1.5, length.out = 1000) theta
En el capítulo anterior, nos enfocamos en construir una distribución posterior a partir de una distribución previa. Sin embargo, surge una pregunta:
¿Es posible estimar la distribución posterior sin una densidad previa?
Para responderla, primero debemos revisar algunas definiciones relacionadas con la muestra y la independencia condicional al valor de un parámetro.
Definición 5.1 (Función de verosimilitud) Dada una muestra
Si consideramos dos valores posibles del parámetro,
El principio de verosimilitud establece que si
Observación. Es más verosímil (realista) que el verdadero valor del parámetro sea
Definición 5.2 (Estimador de máxima verosimilitud (MLE)) Para cada observación
Ejemplo 5.1 Consideremos una muestra
La función de verosimilitud es:
Donde
La log-verosimilitud es una transformación monótona creciente de la función de verosimilitud y se define como:
Para encontrar el MLE de
Para verificar que
Por lo tanto
Ahora generemos 100 valores con una distribución exponencial con parámetro
<- rexp(n = 100, rate = 1)
x <- length(x)
n <- sum(x)
y <- seq(0.5, 1.5, length.out = 1000) theta
ggplot(
data.frame(
theta = theta,
L = theta^(-n) * exp(-y / theta)
),aes(x = theta, y = L)
+
) geom_line(linewidth = 2) +
geom_vline(aes(xintercept = mean(x)), linewidth = 2, color = "red") +
::theme_cowplot() +
cowplotlabs(
title = "Verosimilitud con la media de x",
x = "Theta",
y = "Verosimilitud"
)
ggplot(
data.frame(
theta = theta,
l = -n * log(theta) - y / theta
),aes(x = theta, y = l)
+
) geom_line(linewidth = 2) +
geom_vline(aes(xintercept = mean(x)), linewidth = 2, color = "red") +
::theme_cowplot() +
cowplotlabs(
title = "Log-verosimilitud con la media de x",
x = "Theta",
y = "Verosimilitud"
)
Ejemplo 5.2 Suponga que se está haciendo pruebas RT-PCR para detectar el virus SARS-CoV-2. La prueba es 90% confiable en el siguiente sentido:
En este caso nuestra variable aleatoria es
Nuestra muestra tendria la siguiente forma:
Resumiento el comportamiento de la prueba, tendríamos estas funciones de densidad
Si
Si
Veamos este gráfico que explica lo que está sucediendo:
Entonces los cuadros marcados en rojos, indican donde se maximiza la probabilidad para este ejemplo. Por lo tanto el MLE corresponde a
<- c(0.1, 0.9)
theta
## Define f(0|theta)
<- c(0.9, 0.1)
f_x0
## Define f(1|theta)
<- c(0.1, 0.9)
f_x1
<- data.frame(
df theta = rep(theta, 2),
likelihood = c(f_x0, f_x1),
x = factor(rep(0:1, each = 2))
)
# Gráfico
ggplot(df, aes(x = as.factor(theta), y = likelihood, fill = x)) +
geom_bar(stat = "identity", position = "dodge") +
labs(
title = "Función de verosimilitud para diferentes valores de x",
x = expression(theta),
y = "Verosimilitud",
fill = "Valor de x"
+
) theme_minimal() +
scale_fill_manual(values = c("blue", "red"), labels = c("x = 0", "x = 1"))
Ejemplo 5.3 (Caso 1: Varianza conocida) Consideremos una muestra
La función de verosimilitud para esta muestra es:
La log-verosimilitud es:
Para encontrar el estimador de máxima verosimilitud (MLE) de
Dado que la función es cuadrática en
Ejemplo 5.4 (Caso 2: Varianza desconocida) Si consideramos una muestra
Si fijamos
De la última igualdad, obtenemos que:
Esta es la fórmula para la varianza muestral.
Las condiciones de segundo orden para verificar que es un máximo se dejan como ejercicio para el lector.
Ejemplo 5.5 (Visualización de la función de verosimilitud) Ahora observe como quedaría la verosimilitud si variamos al mismo tiempo
Primero, generemos datos aleatorios de una distribución normal para simular una muestra:
set.seed(123)
<- 100
n <- rnorm(n, mean = 5, sd = 2) data
Luego, definimos la función de log-verosimilitud:
<- function(mu, sigma_sq, data) {
log_likelihood <- length(data)
n -n / 2 * log(2 * pi * sigma_sq) - (1 / (2 * sigma_sq)) * sum((data - mu)^2)
}
Ahora, generamos una grilla de valores para mu y sigma^2 y calculamos la log-verosimilitud para cada combinación:
<- seq(4, 6, length.out = 100)
mu_vals <- seq(2, 5, length.out = 100)
sigma_sq_vals
<- matrix(0, nrow = length(mu_vals), ncol = length(sigma_sq_vals))
ll_vals
for (i in seq_along(mu_vals)) {
for (j in seq_along(sigma_sq_vals)) {
<- log_likelihood(mu_vals[i], sigma_sq_vals[j], data)
ll_vals[i, j]
} }
Y finalmente gráficamos la función de verosimilitud:
<- expand.grid(mu = mu_vals, sigma_sq = sigma_sq_vals)
df $ll <- as.vector(ll_vals)
df
ggplot(df, aes(x = mu, y = sigma_sq, z = ll)) +
geom_raster(aes(fill = ll)) +
geom_contour(aes(z = ll), color = "white") +
geom_point(aes(x = mean(data), y = var(data)), color = "red", size = 4) +
scale_fill_viridis_c() +
labs(
title = "Log-verosimilitud en función de mu y sigma^2",
x = expression(mu),
y = expression(sigma^2),
fill = "Log-verosimilitud"
+
) ::theme_cowplot() cowplot
Ejemplo 5.6 Consideremos una muestra
La función de densidad para esta distribución es:
Donde
La función de verosimilitud para la muestra es:
Note que
El valor de la muestra en la
Por ejemplo, si definimos que
<- runif(100, 0, 2)
x <- length(x)
n
<- seq(1.5, 2.5, length.out = 1000)
theta
<- numeric(1000)
L for (k in 1:1000) {
<- 1 / theta[k]^n * prod(x < theta[k])
L[k]
}
ggplot(data.frame(theta, L), aes(x = theta, y = L)) +
geom_line(linewidth = 2) +
geom_vline(
aes(xintercept = max(x)),
color = "red",
linetype = "dashed",
linewidth = 2
+
) labs(
title = "Función de verosimilitud para diferentes valores de theta",
x = expression(theta),
y = "Verosimilitud"
+
) theme_minimal()
Teorema 5.1 Si
Prueba. Sea
El MLE de
Ejemplo 5.7 Suponga que para el caso exponencial se tiene que el MLE de
¿Qué pasa si
no es biyectiva?
Definición 5.3 (Generalización del MLE) Si
Teorema 5.2 Si
Prueba. Basta probar que
Ejemplo 5.8 Considere una muestra
Observación. Si
Definición 5.4 Un estimador
Denotamos esta convergencia como
Informalmente, la definición anterior dice que a medida que el tamaño de la muestra se vuelve infinito (y la información de la muestra se vuelve cada vez mejor), el estimador estará arbitrariamente cerca del parámetro con alta probabilidad, una propiedad eminentemente deseable. O, cambiando las cosas, podemos decir que la probabilidad de que una secuencia consistente de estimadores no alcance el parámetro verdadero es pequeña.
Recordemos que los estimadores bayesianos son de la forma
El método de los momentos es una técnica de estimación que utiliza los momentos muestrales para estimar los parámetros de una distribución. A continuación, se presentan algunas propiedades y ejemplos de esta técnica.
Ejemplo 5.9 Consideremos una muestra
La función de verosimilitud asociada es:
Para encontrar el estimador de máxima verosimilitud (MLE) de
El problema con este cálculo es que no existe una expresión cerrada para la derivada de la función Gamma. Usualmente se hacen aproximaciones numéricas.
Definición 5.5 Supongamos que tenemos una muestra
y definimos los momentos empíricos como:
El estimador de máxima verosimilitud (MLE) según el método de los momentos es:
Ejemplo 5.10 Del ejemplo anterior,
Ejemplo 5.11 Consideremos una muestra
Para encontrar los estimadores de
De
Teorema 5.3 Supongamos que
Suponga que tenemos una tabla con los siguientes datos, los cuales representan la cantidad de giros hacia la derecha en cierta intersección.
<- c(
(X rep(0, 14),
rep(1, 30),
rep(2, 36),
rep(3, 68),
rep(4, 43),
rep(5, 43),
rep(6, 30),
rep(7, 14),
rep(8, 10),
rep(9, 6),
rep(10, 4),
rep(11, 1),
rep(12, 1)
))
[1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1
[26] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2
[51] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[76] 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
[101] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
[126] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 4 4
[151] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
[176] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5
[201] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
[226] 5 5 5 5 5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
[251] 6 6 6 6 6 6 6 6 6 6 6 6 6 6 7 7 7 7 7 7 7 7 7 7 7
[276] 7 7 7 8 8 8 8 8 8 8 8 8 8 9 9 9 9 9 9 10 10 10 10 11 12
Queremos ajustar esta tabla a una distribución Poisson con función de densidad
Se puede probar que teórico de máxima verosimilitud para
Primero veamos la forma de los datos,
hist(X, main = "Número de giros a la derecha", right = FALSE, prob = TRUE)
Definamos la función correspondiente como
<- length(X)
n <- function(lambda) {
negloglike * lambda - sum(X) * log(lambda) + sum(log(factorial(X)))
n }
Para encontrar el parámetro deseado, basta minimizar la función negloglike
usando el la instrucción de optimización no lineal optim
.
<- optim(par = 0.5, fn = negloglike, hessian = TRUE) lambda_hat
Aquí el valor par = c(0.5)
representa un valor inicial de búsqueda y hessian = TRUE
determina el cálculo explícito de la segunda derivada.
Compare el resultado de lambda.hat$estimate
con mean(X)
el cual sería el estimador por el método de los momentos.
$par lambda_hat
[1] 3.892969
mean(X)
[1] 3.893333