5  Estimación por máxima verosimilitud

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 X1,,Xn que son variables aleatorias independientes e idénticamente distribuidas (i.i.d.) con distribución f(Xθ) y un θ fijo, la función de verosimilitud se define como:

fn(Xθ)=i=1nf(Xiθ)=L(θX).

Si consideramos dos valores posibles del parámetro, θ1 y θ2, y evaluamos la función de verosimilitud para cada uno de ellos, obtenemos: fn(Xθ1)=i=1nf(Xiθ1)=L(θ1X)fn(Xθ2)=i=1nf(Xiθ2)=L(θ2X)

El principio de verosimilitud establece que si fn(Xθ1)>fn(Xθ2), entonces L(θ1X)>L(θ2X).

Observación. Es más verosímil (realista) que el verdadero valor del parámetro sea θ1 que θ2 dada la muestra observada.

Definición 5.2 (Estimador de máxima verosimilitud (MLE)) Para cada observación x en el espacio muestral X, si existe un estimador δ(x) de θ que maximiza la función de verosimilitud fn(xθ), entonces a δ(x) se le denomina estimador de máxima verosimilitud (MLE).

Ejemplo 5.1 Consideremos una muestra X1,,Xn con distribución exponencial Exp(θ), donde la función de densidad es 1θexθ si x>0. Estamos interesados en estimar θ.

La función de verosimilitud es:

fn(Xθ)=1θnexp(i=1nxiθ)=θnexp(yθ).

Donde y=i=1nxi.

La log-verosimilitud es una transformación monótona creciente de la función de verosimilitud y se define como:

(θX)=lnfn(Xθ)=nlnθyθ.

Para encontrar el MLE de θ, derivamos la log-verosimilitud con respecto a θ y resolvemos para θ:

θ(θX)=θ[nlnθyθ]=nθ+yθ2=01θ(n+yθ)=0θ^=yn=X¯n.

Para verificar que θ^ es un máximo, derivamos la log-verosimilitud con respecto a θ dos veces y evaluamos en θ=θ^:

2θ2(θX)=nθ22yθ3|θ=yn=n(yn)22y(yn)3=n3y22n3y2=n3y2<0.

Por lo tanto θ^=X¯n es el MLE de θ.

Ahora generemos 100 valores con una distribución exponencial con parámetro θ=1 y grafiquemos la función de verosimilitud y la log-verosimilitud:

Código
x <- rexp(n = 100, rate = 1)
n <- length(x)
y <- sum(x)
theta <- seq(0.5, 1.5, length.out = 1000)
Código
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") +
  cowplot::theme_cowplot() +
  labs(
    title = "Verosimilitud con la media de x",
    x = "Theta",
    y = "Verosimilitud"
  )

Código
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") +
  cowplot::theme_cowplot() +
  labs(
    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:

  • Si la persona es positiva, entonces la probabilidad que una persona sea efectivamente positiva es 0.9. Si la persona es negativa, entonces la prueba dirá con probabilidad 0.1 que positiva.
  • El mismo razonamiento aplica si una persona es negativa. La prueba dirá con 0.9 probabilidad que es negativa; mientras que si es realmente es positiva dirá con 0.1 probabilidad que es negativa.

En este caso nuestra variable aleatoria es Bernoulli(θ),donde el espacio total parametrico consiste de θ{0.9,0.1}. En este ejemplo 0.1 significa negativo y 0.9 significa positivo.

Nuestra muestra tendria la siguiente forma:

x={1si la prueba es positiva0si no

Resumiento el comportamiento de la prueba, tendríamos estas funciones de densidad

Si x=1, entonces f(1θ)=P(θx=1)={0.1si θ=0.10.9si θ=0.9.

Si x=0, entonces f(0θ)=P(θx=0)={0.9si θ=0.10.1si θ=0.9.

Veamos este gráfico que explica lo que está sucediendo:

Inicio

Resultado de la Prueba

Positiva (x=1)

Negativa (x=0)

P(θ=0.1|x=0) = 0.9

(P(θ=0.9|x=0) = 0.1)

P(θ=0.1|x=1) = 0.1

P(θ=0.9|x=1) = 0.9

Entonces los cuadros marcados en rojos, indican donde se maximiza la probabilidad para este ejemplo. Por lo tanto el MLE corresponde a θ^={0.1si x=00.9si x=1

Código
theta <- c(0.1, 0.9)

## Define f(0|theta)
f_x0 <- c(0.9, 0.1)

## Define f(1|theta)
f_x1 <- c(0.1, 0.9)

df <- data.frame(
  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 X1,,Xn que sigue una distribución normal N(μ,σ2) con σ2 conocida. Estamos interesados en estimar el parámetro μ.

La función de verosimilitud para esta muestra es:

fn(Xμ)=i=1nf(Xiμ)=i=1n12πσ2exp((xiμ)22σ2)=(2πσ2)n/2exp(12σ2i=1n(xiμ)2).

La log-verosimilitud es:

(μX)=lnfn(Xμ)=n2ln(2πσ2)12σ2i=1n(xiμ)2.

Para encontrar el estimador de máxima verosimilitud (MLE) de μ, derivamos la log-verosimilitud con respecto a μ y resolvemos para μ:

μ(μX)=μ[n2ln(2πσ2)12σ2i=1n(Xiμ)2]=1σ2i=1n(Xiμ)=0i=1n(Xiμ)=0i=1nXinμ=0μ^=1ni=1nXi=X¯n.

Dado que la función es cuadrática en μ, tiene un único máximo. Por lo tanto, el MLE de μ es:

μ^=X¯n.

Ejemplo 5.4 (Caso 2: Varianza desconocida) Si consideramos una muestra X1,,Xn que sigue una distribución normal N(μ,σ2) con σ2 desconocida, la función de log-verosimilitud es:

(μ,σ2X)=lnfn(Xμ,σ2)=n2ln(2πσ2)12σ2i=1n(Xiμ)2.

Si fijamos μ=X¯n y derivamos con respecto a σ2, obtenemos:

σ2(μ,σ2X)=σ2[n2ln(2πσ2)12σ2i=1n(XiX¯)2]=n2σ2+12(σ2)2i=1n(XiX¯n)2=0

De la última igualdad, obtenemos que:

n2σ2=12(σ2)2i=1n(XiX¯n)2σ^2=1ni=1n(XiX¯n)2

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 μ y σ.

Primero, generemos datos aleatorios de una distribución normal para simular una muestra:

Código
set.seed(123)
n <- 100
data <- rnorm(n, mean = 5, sd = 2)

Luego, definimos la función de log-verosimilitud:

Código
log_likelihood <- function(mu, sigma_sq, data) {
  n <- length(data)
  -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:

Código
mu_vals <- seq(4, 6, length.out = 100)
sigma_sq_vals <- seq(2, 5, length.out = 100)

ll_vals <- matrix(0, nrow = length(mu_vals), ncol = length(sigma_sq_vals))

for (i in seq_along(mu_vals)) {
  for (j in seq_along(sigma_sq_vals)) {
    ll_vals[i, j] <- log_likelihood(mu_vals[i], sigma_sq_vals[j], data)
  }
}

Y finalmente gráficamos la función de verosimilitud:

Código
df <- expand.grid(mu = mu_vals, sigma_sq = sigma_sq_vals)
df$ll <- as.vector(ll_vals)


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"
  ) +
  cowplot::theme_cowplot()

Ejemplo 5.6 Consideremos una muestra X1,,Xn que son variables aleatorias independientes e idénticamente distribuidas (i.i.d.) con distribución uniforme en el intervalo [0,θ], donde θ>0. Asumimos que cada observación Xi>0.

La función de densidad para esta distribución es:

f(Xθ)=1θ1[0,θ](x)

Donde 1[0,θ](x) es una función indicadora que toma el valor de 1 si x está en el intervalo [0,θ] y 0 en caso contrario.

La función de verosimilitud para la muestra es:

fn(θθ)=i=1nf(Xiθ)=1θni=1n1[0,θ](Xi)=1θni=1n1{0Xiθ}0Xiθi.

Note que fn(θθ) es positivo si y solo si 0X(n)θ. donde X(n)=max{X1,,Xn} es el estadístico de orden máximo.

El valor de la muestra en la i-ésima posición, cuando los datos se ordenan de menor a mayor, se denota X(i) (estadístico de orden). En este caso, X(n)=max{X1,,Xn}. Entonces θ^MLE=X(n).

Por ejemplo, si definimos que θ=2 entonces en ese punto existe una singularidad de la verosimilitud. En términos prácticos, esto significa que la función de verosimilitud no está definida o no es finita en ese punto. Es esencial ser consciente de estas singularidades al trabajar con estimaciones de máxima verosimilitud, ya que pueden afectar los resultados y la interpretación.

Código
x <- runif(100, 0, 2)
n <- length(x)

theta <- seq(1.5, 2.5, length.out = 1000)

L <- numeric(1000)
for (k in 1:1000) {
  L[k] <- 1 / theta[k]^n * prod(x < theta[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()

5.1 Propiedades de los estimadores de máxima verosimilitud

5.1.1 Propiedad de invarianza

Teorema 5.1 Si θ^ es el MLE de θ^ y si g es biyectiva, entonces g(θ) es el MLE de g(θ).

Prueba. Sea Ω el espacio paramétrico g(Ω). Dado que g es biyectiva entonces defina la inversa de g como θ=g1(ψ),ψΩ. Al reparametrizar la función de verosimilitud, tenemos:

fn(xθ)=fn(xg1(ψ))

El MLE de ψ es ψ^ y satisface que fn(xg1(ψ^)) es máximo. Dado que fn(xθ) se maximiza cuando θ=θ^, entonces fn(xg1(ψ)) se maximiza cuando θ^=g1(ψ) para algún ψ. Por lo tanto, si θ^=g1(ψ^) entonces ψ^=g(θ^).

Ejemplo 5.7 Suponga que para el caso exponencial se tiene que el MLE de θ es θ^=X¯n. Si se quiesiera el MLE de 1θ entonces se tiene que g(θ)=1θ y g1(ψ)=1ψ. Por lo tanto, el MLE de 1θ es 1X¯n.

¿Qué pasa si h no es biyectiva?

Definición 5.3 (Generalización del MLE) Si g es una función de θ y G la imagen de Ω bajo g. Para cada tG defina Gt={θ:g(θ)=t} Defina L(t)=maxθGtlnfn(x|θ). El MLE de g(θ) es igual a un valor t^ y satisface L(t^)=maxtGL(t).

Teorema 5.2 Si θ^ es el MLE de θ entonces g(θ^) es el MLE de g(θ) (g es arbitraria).

Prueba. Basta probar que L(t^)=lnfn(x|θ^). Se cumple que θ^Gt^. Como θ^ maximiza fn(x|θ) θ, también lo hace si θGt^. Entonces t^=g(θ^) (no pueden existir 2 máximos en un conjunto con la misma imagen).

Ejemplo 5.8 Considere una muestra X1,,XnN(μ,σ2).

  • El MLE de μ es X¯n.
  • El MLE de σ2 es σ^2=1ni=1n(XiX¯n)2.
  • Si g(μ,σ2)=σ2 entonces el MLE de σ2 es σ^2=1ni=1n(XiX¯n)2.
  • Si g(μ,σ2)=σμ entonces el MLE de σμ es σ^X¯n.
  • Si g(μ,σ2)=μ2+σ2 entonces el MLE de μ2+σ2 es X¯n2+σ^2.

Observación. Si θMLE de θ, entonces h(θMLE) es el MLE de h(θ).

5.1.2 Consistencia

Definición 5.4 Un estimador θ^n proveniente de la muestra X1,,Xn se dice consistente si para cada ϵ>0 y θΩ limnP(|θ^nθ|<ϵ)=1.

Denotamos esta convergencia como

θ^nnPθ.

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 W1E[Previa]+W2θ^MLE. El estimador bayesiano “combina” la esperanza de la previa y el θ^MLE. El θ^MLE hereda la consistencia del estimador bayesiano si n tiende a infinito. Eso si, se debe garantizar que limnVar(θ^n)=0 y limnBias(θ^n)=0. Bajo estas condiciones entonces

θ^MLEnPθ.

5.2 Cálculo numérico

5.2.1 Método de los momentos

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 X1,,Xn que sigue una distribución Gamma Γ(α,1). Estamos interesados en estimar el parámetro α. La función de densidad de la distribución Gamma es dada por:

f(xα)=1Γ(α)xα1ex

La función de verosimilitud asociada es:

fn(xα)=1Γ(α)ni=1nxiα1exi

Para encontrar el estimador de máxima verosimilitud (MLE) de α, derivamos la log-verosimilitud con respecto a α y la igualamos a cero:

α(αX)=α[nlnΓ(α)+(α1)i=1nlnXii=1nXi]=n1Γ(α)αΓ(α)+i=1nlnXi=0

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 X1,,Xn que sigue una distribución F, indexada por un parámetro θ que reside en Rk, y que tiene al menos k momentos finitos. Para j=1,,k, definimos μj(θ)=E[X1jθ]. Asumimos que la función μ(θ)=(μ1(θ),,μ2(θ)) es biyectiva. Definimos M como la inversa de μ, es decir,

M(μ(θ))=θ=M(μ1(θ),,μ2(θ))

y definimos los momentos empíricos como:

mj=1ni=1nXij,j=1,,k.

El estimador de máxima verosimilitud (MLE) según el método de los momentos es:

θ^=M(m1,,mk).

Ejemplo 5.10 Del ejemplo anterior, μ1(α)=E[x1|α]=α.Dado que m1=x^n, el sistema por resolver es μ1(α)=m1α=x¯n El estimador por método de momentos es α^=X¯n.

Ejemplo 5.11 Consideremos una muestra X1,,Xn que sigue una distribución Gamma(α,β). Sabemos que la media es μ=α/β y la varianza es σ2=α/β2. Por lo tanto el segundo momento es

E[X2]=Var(X)+E[X]2=αβ2+(αβ)2=α(α+1)β2

Para encontrar los estimadores de α y β, necesitamos resolver el siguiente sistema de ecuaciones:

{μ1(θ)=αβ=X¯n=m1(1)μ2(θ)=α(α+1)β2=m2(2)

De (1), α=m1β. Sustituyendo en (2),

m2=m1β(m1β+1)β2=m12+m1β=m2m2m12=m1β. De esta manera, β^=m1m2m12,α^=m12m2m12

Teorema 5.3 Supongamos que X1,,Xn son variables aleatorias independientes e idénticamente distribuidas (i.i.d) con una distribución indexada por θRk. Si los k momentos teóricos son finitos para todo θ y M es una función continua, entonces el estimador obtenido mediante el método de los momentos es consistente.

5.3 Laboratorio

Suponga que tenemos una tabla con los siguientes datos, los cuales representan la cantidad de giros hacia la derecha en cierta intersección.

Código
(X <- c(
  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

P(X=x)=λxeλx!

Se puede probar que teórico de máxima verosimilitud para λ es X (Tarea). Queremos estimar este parámetro alternativamente maximizando la función de verosimilitud.

Primero veamos la forma de los datos,

Código
hist(X, main = "Número de giros a la derecha", right = FALSE, prob = TRUE)

Definamos la función correspondiente como log(P(X=x))

Código
n <- length(X)
negloglike <- function(lambda) {
  n * lambda - sum(X) * log(lambda) + sum(log(factorial(X)))
}

Para encontrar el parámetro deseado, basta minimizar la función negloglike usando el la instrucción de optimización no lineal optim.

Código
lambda_hat <- optim(par = 0.5, fn = negloglike, hessian = TRUE)

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.

Código
lambda_hat$par
[1] 3.892969
Código
mean(X)
[1] 3.893333