Máxima Verosimilitud (Likelihood)
Objetivos:
- Relacionar el likelihood de datos con parámetros de una distribución
- Familiarizarse con curvas de likelihood y log-likelihood
- Comprobar que la forma de la curva de likelihood depende del tamaño muestreal
Vimos que las funciones de densidad nos dan la probabilidad de observar un valor determinado de una variable aleatoria (para variables continuas nos dan la probabilidad alrededor de un valor). Por ejemplo, para una distribución de Poisson, la probabilidad para un valor \(y\) es \(\frac{\lambda^{y} \exp(- \lambda)}{y!}\) (donde \(!\) es el factorial de \(y\)). Entonces, para una variable aleatoria que suponemos tiene una distribución de Poisson con \(\lambda=2\), la probabilidad de observar un cero (\(y=0\)) es de \(\frac{2^{0} \exp(-2)}{0!} = \exp(-2)\). En R usamos la función dpois para hacer este tipo de cuentas. Para este ejemplo escribimos dpois(0, lambda = 2).
Normalmente tenemos un conjunto de datos, por ejemplo números de renovales por metro cuadrado de muestreo en un bosque: \(y = 1, 0, 2, 1, 0, 0, 1, 1, 2, 2\). Si suponemos que estos datos se pueden describir como el resultado de un proceso aleatorio de tipo Poisson, y que además las observaciones son independientes, podemos calcular la probabilidad del conjunto de datos como el producto de las probabilidades de cada una de las observaciones:
\[ p(y)= \prod_i^n{ \frac{\lambda^{y_i} \exp(- \lambda)} {y_i!} } \]
Podemos ver que el producto (\(\prod\)) se calcula usando \(i\) como índice para cada observación del set de datos \(y\), desde \(i=1\) hasta \(i=n\) que es el tamaño de \(y\). Esta probabilidad de observar los datos es lo que llamamos likelihood o verosimilitud. Para hacer la cuenta de arriba tenemos que definir además el parámetro \(\lambda\). La probabilidad de los datos depende entonces no sólo de los valores presentes en el set de datos sino también de la distribución utilizada y el valor del o los parámetros de la distribución. De manera más general, la expresión es un modelo de datos (data model, data generating process), es decir una descripción probabilística de cómo se generan los datos. Más adelante veremos como combinando procesos determinísticos y estocásticos podemos armar modelos de datos realistas y útiles. Claramente, esto es siempre una simplificación de la realidad!
Veamos ahora como calcular la probabilidad de estos datos en R suponiendo que \(\lambda = 1\):
y <- c(1, 0, 2, 1, 0, 0, 1, 1, 2, 2)
lambda <- 1
p.y <- dpois(y, lambda = lambda)
p.y## [1] 0.3678794 0.3678794 0.1839397 0.3678794 0.3678794 0.3678794 0.3678794
## [8] 0.3678794 0.1839397 0.1839397prod(p.y)## [1] 5.674991e-06Como las probabilidades son números que van entre \(0\) y \(1\), el producto de varias probabilidades resulta en números muy pequeños. Vemos en este caso que la probabilidad del conjunto de datos es muy baja (\(5.6749912\times 10^{-6}\)). Esto no debería sorprendernos porque esa es la probabilidad de encontrar exactamente esos datos en particular (sin importar el orden). Entonces, para evitar problemas numéricos que pueden aparecer cuando la computadora trabaja con números muy chicos, preferimos trabajar con la suma de los logaritmos de las probabilidades, es decir con el logaritmo del likelihood. En R, podemos usar la opción log = TRUE dentro de las funciones de densidad.
y <- c(1, 0, 2, 1, 0, 0, 1, 1, 2, 2)
lambda <- 1
ll <- dpois(y, lambda = lambda, log = TRUE)
sum(ll)## [1] -12.07944Una vez que logramos hacer estas cuentas, podemos preguntarnos cómo cambia la probabilidad de observar el set de datos si cambiamos el valor de \(\lambda\). Probemos por ejemplo con \(\lambda = 2\)
y <- c(1, 0, 2, 1, 0, 0, 1, 1, 2, 2)
lambda <- 2
sum(dpois(y, lambda = lambda, log = TRUE))## [1] -15.14797Con un \(\lambda\) de \(2\) el log likelihood de los datos empeora (es más chico), de manera que nos hace pensar que los datos son más consistentes con \(\lambda=1\) que con \(\lambda=2\). Siguiendo este razonamiento, podemos continuar probando distintos valores de \(\lambda\) hasta encontrar el que hace más probable al set de datos. Esta es la idea detrás de los análisis de máxima verosimilitud (a falta de mejor nombre) o de maximum likelihood (si estamos dispuestos a usar palabras en inglés). Veamos gráficamente cómo cambia el log-likelihood a medida que cambia \(\lambda\):
y <- c(1, 0, 2, 1, 0, 0, 1, 1, 2, 2)
lambda <- seq(0, 3, by = 0.1)
ll <- numeric(length(lambda))
for (i in 1:length(lambda)) {
ll[i] <- sum(dpois(y, lambda = lambda[i], log = TRUE))
}
op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)
plot(lambda, ll, pch = 21, bg = "grey")
par(op)Podemos ver que la curva de logL llega a un máximo cuando \(\lambda = 1\). Podemos decir entonces que \(\lambda = 1\) maximiza la probabilidad de observar nuestro set de datos.
Normalmente tenemos más de un parámetro en nuestro modelo de datos, y usamos algún algoritmo de búsqueda para encontrar la combinación de parámetros que hace más probables a nuestros datos (por tradición, los algoritmos buscan minimizar el menos logaritmo de la verosimilitud de nuestros datos). Esta idea de buscar combinaciones de parámetros que maximizan la verosimilitud de los datos se conoce como maximum likelihood estimation y muchas veces aparece abreviado como MLE.
Ejemplo: Remoción de Frutos
Mariano Rodriguez-Cabal estudió la remoción de frutos de quintral por parte del monito del monte en \(3\) sitios pareados de bosque continuo y fragmentado con la idea de ver si la fragmentación afectaba la interacción entre este marsupial y el muérdago. El trabajo está disponible aquí
Los datos son de número de frutos removidos en la temporada reproductiva de un total disponible. Podemos usar una distribución Binomial para capturar la variabilidad (estocasticidad) en el proceso de remoción de frutos. Es decir, suponemos que la distribución Binomial tiene sentido para representar el proceso que genera los datos. En este ejemplo vemos otra manera de cargar los datos desde un repositorio online.
url <- "https://sites.google.com/site/modelosydatos/quintral.txt"
quintral <- read.table(url, header = TRUE)Para ver como están organizados los datos podemos hacer:
str(quintral)## 'data.frame': 70 obs. of 4 variables:
## $ bosque : Factor w/ 2 levels "c","f": 1 1 1 1 1 1 1 1 1 1 ...
## $ sitio : Factor w/ 3 levels "Campanario","Llao-Llao",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Removidos: int 30 40 11 33 35 25 48 13 25 31 ...
## $ Frutos : int 46 51 32 49 46 37 53 31 33 44 ...Vemos que se trata de un data frame que es la manera preferida de R para organizar datos. Los data frames tienen las observaciones en filas y las variables en columnas. Las variables individuales se pueden acceder escribiendo el nombre del data frame seguido del signo $ y el nombre de la variable. Por ejemplo, si queremos ver los valores de frutos disponibles en las tres primeras observaciones hacemos:
quintral$Frutos[1:3]## [1] 46 51 32Ahora calculemos el likelihood para la primera observación del set de datos considerando que la tasa de remoción es de \(0.5\) (si no nos acordamos cómo es la distribución Binomial o cómo usarla en R podemos pedir ayuda con ?dbinom). En lo que sigue, vamos a llamar \(\theta\) a la tasa (probabilidad) de remoción de frutos.
theta <- 0.5
dbinom(quintral$Removidos[1], size = quintral$Frutos[1], prob = theta)## [1] 0.01408998Normalmente trabajamos con log-likelihoods, por eso especificamos log = TRUE:
dbinom(quintral$Removidos[1], size = quintral$Frutos[1], prob = theta, log = TRUE)## [1] -4.262292¿Cómo cambia la probabilidad de observar ese dato cuando cambiamos el parámetro de probabilidad de éxito (remoción) en la Binomial? Como en este caso el likelihood depende de un solo parámetro, podemos responder a esta pregunta gráficamente.
op <- par(mfrow = c(1, 2), lwd = 2, bty = "n", cex.lab = 1.5, font.lab = 2,
cex.axis = 1.3, las = 1)
theta <- seq(0, 1, length = 1000)
plot(theta, dbinom(quintral$Removidos[1], size = quintral$Frutos[1], prob = theta,
log = FALSE), type = "l", ylab = "Likelihood", xlab = expression(theta))
plot(theta, -dbinom(quintral$Removidos[1], size = quintral$Frutos[1], prob = theta,
log = T), type = "l", ylab = "-logL", xlab = expression(theta))
par(op)También podemos ver cómo cambia el Likelihood cuando cambiamos el número de observaciones que consideramos.
op <- par(mfrow = c(2, 3), lwd = 2, bty = "n", cex.lab = 1.5, las = 1, mar = c(5,
4, 1, 1), font.lab = 2, cex.axis = 1.3)
theta <- seq(0, 1, length = 1000)
plot(theta, dbinom(quintral$Removidos[1], size = quintral$Frutos[1], prob = theta,
log = FALSE), type = "l", xlab = "", ylab = "Likelihood", main = "n = 1")
nll2 <- array(0, length(theta))
for (i in 1:length(theta)) {
nll2[i] <- -sum(dbinom(quintral$Removidos[1:2], size = quintral$Frutos[1:2],
prob = theta[i], log = TRUE))
}
plot(theta, exp(-nll2), type = "l", xlab = "", ylab = "", main = "n = 2")
nll10 <- array(0, length(theta))
for (i in 1:length(theta)) {
nll10[i] <- -sum(dbinom(quintral$Removidos[1:10], size = quintral$Frutos[1:10],
prob = theta[i], log = TRUE))
}
plot(theta, exp(-nll10), type = "l", xlab = "", ylab = "", main = "n = 10")
plot(theta, -dbinom(quintral$Removidos[1], size = quintral$Frutos[1], prob = theta,
log = TRUE), type = "l", xlab = "", ylab = "-logL", main = "n = 1")
plot(theta, nll2, type = "l", xlab = expression(theta), ylab = "", main = "n = 2",
ylim = c(0, 100))
plot(theta, nll10, type = "l", xlab = "", ylab = "", main = "n = 10", ylim = c(0,
100))
par(op)Pregunta (1): ¿Cómo interpretan estos cambios?
Buscando el máximo (o más bien el mínimo…)
Queremos encontrar el valor de tasa de remoción por fruto que maximiza el likelihood (o minimiza el negative log-likelihood) del conjunto de datos. Una opción es hacerlo a lo bruto:
theta <- seq(0, 1, length = 50)
nll <- array(NA, length(theta))
for (i in 1:length(theta)) {
nll[i] <- -sum(dbinom(quintral$Removidos, size = quintral$Frutos, prob = theta[i],
log = TRUE))
}
theta[which(nll == min(nll))]## [1] 0.5918367Para este ejemplo en que usamos la distribución Binomial y queremos estimar la probabilidad de éxito, tenemos una solución analítica y podemos compararla con la que acabamos de calcular:
sum(quintral$Removidos)/sum(quintral$Frutos)## [1] 0.5968072Pregunta (2): ¿Qué pasa si cambiamos la resolución del vector de tasas de remoción que usamos para calcular el likelihood?
Veamos ahora un ejemplo con dos parámetros. Queremos modelar la distancia diaria de desplazamiento (en km) de unos elk a los que les pusimos collares con GPS. Si consideramos que la distribución Weibull es una buena opción para estos datos, tenemos que estimar un parámetro de forma y uno de escala. Entonces, la búsqueda la tenemos que hacer en dos dimensiones.
elk <- read.table("https://sites.google.com/site/modelosydatos/aelk.txt", header = TRUE)
step <- elk$km_day
# removemos NAs
step = step[!is.na(step)]
scale <- seq(0.4, 1.4, length.out = 100)
shape <- seq(0.4, 1.2, length.out = 100)
nll <- matrix(NA, 100, 100)
for (i in 1:100) {
for (j in 1:100) {
nll[i, j] <- -sum(dweibull(step, shape = shape[i], scale = scale[j],
log = TRUE))
}
}
contour(shape, scale, nll, xlab = "Shape", ylab = "Scale", levels = seq(0, 3000,
by = 10))
tmp <- which(nll == min(nll), arr.ind = TRUE)
points(shape[tmp[1]], scale[tmp[2]], pch = 19)
Búsquedas con Algoritmos
En general, queremos estimar modelos que tienen varios parámetros y usamos algoritmos de búsqeda en lugar de hacer estas cuentas directas. Si tienen tiempo y ganas, pueden ver cómo usar algoritmos de búsqueda aquí
Juan Manuel Morales. 6 de Septiembre del 2015. Última actualización: 2019-02-21