Hace un montón que no me paso por aquí. Últimamente la vida me ha regalado muchas cosas y, por mucho que predique sobre aprendizaje Bayesiano, me temo que no hay forma de que yo aprenda a decir que no. Y de eso justo venia a hablaros, no de como decir que no, que no soy yo la más idónea, sino de Estadística Bayesiana. Pero ahora sí, ya no me iré por las ramas, hoy vengo a contaros cómo hacer un análisis Bayesiano completo.
¿Me dejas que te cuente?
Empezaré diciendo que esta no será la entrada más fácil de leer de las que he escrito o escribiré. Aun así espero que os resulte útil. Y ahora sí, sin más dilación, vayamos a ello.
Empecemos por lo sencillo
Para poder contarte como va esto de hacer un análisis Bayesiano vamos a recurrir a un modelo sencillo. Vamos a empezar por un modelo que solo tenga un parámetro desconocido.
Imagina que una tarde de primavera, mientras los niños juegan en el parque, te dedicas a observar el tiempo que pasa entre que dos pájaros que descansan sobre un cable de la luz emprenden el vuelo.
Este tipo de variables suelen tener una distribución conocida como exponencial y será lo que denominaremos nuestro modelo. La forma de esta distribución va a depender de un parámetro al que llamamos tasa y que será la inversa del tiempo esperado entre dos sucesos. Por entender un poco lo que significa, si el tiempo esperado (también llamado media) es de 6 minutos, la tasa será 1/6 y estamos diciendo que los tiempos observados se moverán alrededor de los 6 minutos, siguiendo la siguiente distribución:
Ahora bien, esa tasa es un parámetro desconocido que nos interesa determinar. Vamos a llamarlo 
Sin embargo, la visión bayesiana de la estadística sí nos permite tratar a
¿Qué sabemos?
El primer paso será establecer sobre
Durante mucho tiempo, la única opción valida eran lo que se conocen como previas conjugadas. Este tipo de distribuciones aprovechan la forma matemática del modelo de forma que al aplicar el teorema de bayes, la distribución a posteriori sea de la misma familia.
En el caso concreto que nos ocupa, la distribución conjugada es una gamma. Esta distribución depende dos parámetros, 
Rara, no me lo vais a negar… pero es la conjugada, la que es más fácil utilizar.
¿Qué nos dicen los datos?
Ahora que ya tenemos nuestra previa lo que hacemos es observar. Tras hora y media de parque has conseguido ver unos 15 pájaros alzar el vuelo, lo que te da un total de 14 tiempos entre sus salidas.
Una vez observados los datos, toca resumir la información que estos contienen a través de lo que llamamos verosimilitud y que ya os conté en este otro post. Cuando las observaciones se obtienen de forma independiente, como es este caso (suponiendo que el vuelo de cada pájaro, no está relacionado con el resto), la verosimilitud se obtiene al multiplicar la función de densidad para cada dato
exp_lik <- function(data, tasa) {
aux <- lapply(data, dexp, rate=tasa)
#Se calcula la densidad para cada dato
exp(sum(log(unlist(aux))))
#Se aplica el truco de sumar los logaritmos en lugar de multiplicar directamente los valores
}
#exp_lik(data=datos$tiempo, mean=10)
C <- 10^18
# Utilizamos una constante para poder dibujarla. No está escogida siguiendo criterios matemáticos, solo estéticos
m <- seq(0,1,by=0.01)
y <- C*unlist(lapply(m, exp_lik,data =datos$tiempo))En concreto, lo que hacemos al calcular la verosimilitud es estudiar lo bien que se adecuan diferentes exponenciales, entre las que cambia la media, a los datos observados (figura de la izquierda). Esa información se resume en la figura de la derecha y se observa que la máxima verosimilitud se obtiene para la tasa observada que es 0.162 y que equivale a un tiempo medio de 6.157.
Ahora ya tenemos todo lo necesario y solo nos falta unirlo mediante el teorema de Bayes para llegar a la distribución posterior.
Acumulando información.
En particular, el teorema de Bayes en este caso podemos escribirlo como:
Aquí, 
La tercera cantidad en discordia, que en la fórmula aparece como 
Esta función es el verdadero quebradero de cabeza de toda esta historia la que, más allá del tabú impuesto al teorema de Bayes, propició que se tardase mucho mucho tiempo en llegar a poder aplicar la estadística bayesiana tal y como hoy se hace.
En el caso que nos ocupa, y siempre que usemos previas conjugadas, la ventaja es que la distribución posterior, es decir 
¿Qué sabemos ahora?
En particular, se puede ver matemáticamente que los nuevos parámetros serán
y
Con ello, el resultado final puede verse reflejado en la siguiente gráfica, donde la curva roja sera la distribución a posteriori del parámetro
Así sabemos, por ejemplo, que el valor esperado para la tasa a posteriori es 0.163, que equivale a un tiempo medio entre vuelos de 6.15 y está entre la tasa observada 0.162 y la tasa que consideramos a priori, 1/5 = 0.2. Además, utilizando la función qgamma de R podemos calcular los cuantiles de la distribución. Con ello, estaremos dando un intervalo de credibilidad que nos diga por donde anda el valor de 
Vale, hasta aquí bien, ya podemos contar más cosas de
Pero lo más bonito e interesante de esta metodología es que da igual lo complicado que sea el modelo que queramos poner sobre los datos. Solo tenemos que identificar los parámetros desconocidos, poner una previa sobre ellos que capte la información disponible y dejar que sean los datos los que actualicen dicha información.
Pero claro… las cosas no son siempre tan fáciles.
¿Y si no puedo usar conjugadas?
Existen muchos ejemplos donde las previas conjugadas no son viables, no existen o, incluso, no son las más adecuadas para representar el conocimiento que tenemos sobre los parámetros de interés (que cuando son más de uno la cuestión puede complicarse). De hecho, si pensamos en nuestro ejemplo, hemos estado a vueltas con la tasa que es mucho menos intuitivo que usar la media directamente… pero claro, la previa conjugada era justo para la tasa. Y si queremos la media, ¿qué pasaría?
Pues lo cierto es que en la mayoría de los casos la parte del numerador del teorema de Bayes sigue siendo fácil de obtener. Sin embargo, lo que sucede es que la distribución que se obtiene no es una de las conocidas y la constante del denominador no se puede calcular de forma analítica (esto quiere decir que no tenemos una fórmula bonita para ella). Tenemos entonces un problema pues ya no podemos obtener de manera automática la media, la varianza o los intervalos de credibilidad… ¿qué hacemos entonces? Pues recurrir a una gran aliada: La simulación.
En términos generales, simular consiste en obtener valores de una variable de forma que, a partir de ellos, podamos entender cómo se comporta su distribución.
En el caso de la distribución posterior de nuestro parámetro o parámetros existen diferentes formas de simular, pero quizás la más conocida son los Métodos Montecarlo por Cadenas de Markov (MCMC). No entraré en detalle hoy en estos métodos, me los reservo para otro día, pero vamos a ver como usar algún software de los que permiten obtener estas simulaciones de forma bastante intuitiva.
En concreto, vamos a ver como podríamos abordar nuestro ejemplo pero utilizando un programa llamado Just another Gibbs Sampler (JAGS) –haciendo referencia a uno de los métodos usados en MCMC, el método Gibbs (Gelfand, 2000)– a través del paquete de R rjags.
Para usar JAGS, debemos descargarlo de la siguiente página web y después instalamos la libreria rjags desde R. Existen otras librerías para manejar el paquete y podéis encontrar un completo manual de uso aquí.
Una vez hecho esto, lo primero que tendremos que hacer será especificar el modelo con el que vamos a trabajar así como las distribuciones a priori sobre los parámetros desconocidos. Esto lo haremos sobre un archivo de texto que podemos escribir directamente desde R.
A continuación mostramos como sería si en lugar de trabajar con la tasa trabajásemos directamente con el tiempo medio. Además, a priori, vamos a pensar que puede tomar cualquier valor entre 1 y 10 minutos para lo que utilizamos una previa uniforme en ese intervalo.
library(rjags)
modelString = "
model{
# Verosimilitud
for(i in 1:n){
tiempo[i] ~ dexp(lambda)
#Modelo exponencial
}
# Distribuciones a priori
lambda <- 1/media
media ~ dunif(1,10)
}"
writeLines( modelString , con="model.bug" )Por supuesto el modelo puede ser todo lo complicado que necesitemos. Por poner un ejemplo, si tuviésemos una recta de regresión con dos variables explicativas, el modelo podría escribirse como:
modelString = "
model{
# Verosimilitud
for(i in 1:n){
y[i] ~ dnorm(media[i], precision)
#la normal se expresa en terminos de la precisión
media[i] <- b0 + b1*X1[i] + b2*X2[i]
}
# Distribuciones a priori
# Para los coeficientes ponemos normales con precisión baja, varianza alta, indicando que sabemos poco sobre ese parámetro:
b0 ~ dnorm(0.0, 0.0001)
b1 ~ dnorm(0.0, 0.0001)
b2 ~ dnorm(0.0, 0.0001)
desviacion <- 1.0/sqrt(precision)
#transformamos la precisión en desviación estándar
#Le ponemos una distribución previa uniforme
desviación ~ dunif(0,3)
}"
writeLines( modelString , con="model.bug" )Pero vamos al lío. Una vez especificado el modelo y las previas, ahora tenemos que poner el modelo a simular. Para ello debemos establecer unos valores de partida para las simulaciones.
Llegados a este punto, cabe mencionar que los métodos MCMC necesitan un periodo para que las simulaciones que realizan sean realmente representativas de la distribución a posteriori. A este periodo se le conoce como burning. Además, para estar seguras de que esto ha sucedido, hace falta simular al menos un par de conjuntos de valores, a los que llamamos cadenas, empezando en valores iniciales diferentes. Si todas las cadenas convergen a valores similares, ya lo tendremos.
Una de las formas de establecer esos valores iniciales distintos para cada cadena es simulándolos mediante una función como la siguiente.
#Valores iniciales
init_values <- function(){
list(media = runif(1,4,6))
#Simulo un valor entre 4 y 6 (esto no es la previa)
}
Lo siguiente será configurar los datos para introducirlos en el modelo:
# Datos
jags_data <- list(tiempo = datos$tiempo,
n = length(datos$tiempo))Y ahora toca comunicarle a JAGS cual será el modelo, los datos, los valores iniciales y el número de cadenas, en este caso 3. En este momento, que llamamos habitualmente compilar, la librería se encarga de comprobar que todo está correcto:
model <- jags.model(file = "model.bug",data = jags_data,
inits = init_values,
n.chains = 3)Y llegados a este punto podemos empezar a simular. Primero unas cuantas simulaciones que supondrán el proceso de burning:
update(model, n.iter=5000)A continuación, algunas simulaciones más que serán con las que nos quedemos.
samples <- coda.samples(model, variable.names=c("media"), n.iter=10000, thin = 10)Cabe mencionar que lo que aparece como thin=10 es lo que se conoce como proceso de thining o adelgazamiento. La idea es quedarse solo con algunas de las simulaciones. En este caso, 1 de cada 10. Esto es porque cada nueva simulación depende de la anterior y podría ser que existiese lo que se conoce como autocorrelación. Cuando esto sucede, puede que los valores se queden estancados durante un tiempo en valores muy similares, haciendo que el resultado no sea todo lo representativo que debiese de la distribución que estamos intentando simular.
El proceso de simulación puede ser más o menos largo dependiendo, entre otras cuestiones, de la complejidad del modelo y del número de observaciones.
Al terminar tendremos 3000 observaciones, 1000 por cada una de las 3 cadenas, y es importante que observemos si estas han convergido. Para ello podemos utilizar el comando plot sobre el objeto samples
plot(samples)
El resultado son dos gráficos. El de la izquierda muestra como las 3 cadenas simuladas (roja, negra y verde) mezclan perfectamente ocupando una especie de rectángulo horizontal. El de la derecha, por su parte, intenta aproximar como sería la distribución posterior para la media.
Además, podemos resumir todas las características de esta distribución mediante el comando
summary(samples)##
## Iterations = 6010:16000
## Thinning interval = 10
## Number of chains = 3
## Sample size per chain = 1000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## 6.67235 1.47780 0.02698 0.02699
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## 4.129 5.557 6.577 7.702 9.641o directamente usando los valores simulados
sim <- unlist(samples)
# la media de la posterior será aproximadamente la media de las simulaciones
mean(sim)## [1] 6.67235# lo mismo para la desviación estándar
sd(sim)## [1] 1.477799# o los cuantiles
quantile(sim)## 0% 25% 50% 75% 100%
## 2.684764 5.557457 6.577261 7.701866 9.994424Por último, podemos volver a visualizar la previa la posterior y la verosimilitud, para entender que está pasando:
Y, con todo esto, hemos llegado al final. Y gracias a la distribución a posteriori, no solo sabemos que el tiempo medio que separa el vuelo de dos pájaros está entorno a los 6.672 minutos, si no que tenemos una idea aproximada de como se comporta este tiempo medio y la probabilidad de que se encuentre en un rango de valores dado.
Si has llegado leyendo hasta aquí. Gracias. Espero que te sea útil y te sirva para empezar a aplicar la Estadística Bayesiana en tu investigación. Estaré deseando leer vuestros progresos.
Referencias
- Gelfand, A. E. (2000). Gibbs Sampling. Journal of the American Statistical Association, 95(452), 1300–1304. https://doi.org/10.2307/2669775
- Manual de uso de JAGS