Estadística Bayesiana

En esta entrada voy a crear un conjunto de datos donde dos variables tienen una correlación muy alta, ajustar un modelo de regresión y obtener la siguiente representación de la distribución a posteriori de los coeficientes, donde se aprecia el efecto de la correlación entre x1 y x2. El código, library(mvtnorm) library(rstan) library(psych) n <- 100 corr_coef <- .9 x <- rmvnorm(n, c(0, 0), sigma = matrix(c(1, corr_coef, corr_coef, 1), 2, 2)) plot(x) x1 <- x[,1] x2 <- x[,2] x3 <- runif(n) - 0.5 y <- 1 + .4 * x1 - .2 * x2 + .1 * x3 + rnorm(n, 0, .1) summary(lm(y ~ x1 + x2 + x3)) stan_code <- " data { int N; vector[N] y; vector[N] x1; vector[N] x2; vector[N] x3; } parameters { real a; real a1; real a2; real a3; real sigma; } model { a ~ cauchy(0,10); a1 ~ cauchy(0,2.5); a2 ~ cauchy(0,2.5); a3 ~ cauchy(0,2.5); y ~ normal(a + a1 * x1 + a2 * x2 + a3 * x3, sigma); }" datos_stan <- list( N = n, y = y, x1 = x1, x2 = x2, x3 = x3 ) fit2 <- stan(model_code = stan_code, data = datos_stan, iter = 10000, warmup = 2000, chains = 2, thin = 4) res <- as.data.frame(fit2) pairs.panels(res[, c("a", "a1", "a2", "a3", "sigma")])

Más sobre lo de ayer. O más bien, una justificación por analogía. Con monedas. Tiras una moneda 100 veces y obtienes 60 caras. Tienes una priori $B(a,b)$ (beta). Tomas una muestra de valores $p_i$ con esa distribución y para cada una de ellas repites el experimento, es decir, obtienes lo que en R se expresaría de la forma rbinom(1, 100, p[i]) Si te quedas los valores $p_i$ tales que esa simulación es 60, enhorabuena, tienes una muestra de la distribución a posteriori. ...

Que quiere decir approximate Bayesian computation. Es un truco para pobres y desafortunados que no pueden quitarle la A a BC y usar directamente cosas como Stan o similares. El que no quiera prioris, además, puede usar el ABC para estimar la forma de la verosimilitud alrededor de una estimación puntual. Por supuesto, el objetivo es obtener una estimación de la posteriori para poder medir la incertidumbre de parámetros, etc. La idea es que se dispone de unos datos, $X$ y un mecanismo de generación de datos $X^\prime = f(\theta)$, donde $\theta$ es un vector de parámetros. ...

Existe un paquete muy curioso en CRAN, rSARP para diseñar, optimizar y comunicar la evolución de planes de búsqueda y/o rescate (p.e., de un niño desaparecido en un monte). Es particularmente interesante porque este tipo de problemas lo tienen todo: desde distribuciones a priori (sobre dónde es más probable encontrar lo que se busca) hasta la decisión final (explórese tanto aquí y tanto allá) teniendo en cuenta restricciones de tiempo y recursos. ...

A primeros de julio impartí un curso de estadística bayesiana aplicada con Stan. Tengo que examinar a los alumnos y he aquí el primero de los ejercicios: En un país, se extrae una muestra de 2000 hombres y mujeres con la siguiente distribución: men <- 170 + 3 * rt(1000, 6) women <- 160 + 2 * rt(1000, 5) heights <- c(men, women) Ajusta una distribución (una mezcla de dos distribuciones de Student) usando los datos anteriores, i.e., heights. Puedes suponer conocidos: Los pesos de la mezcla (0.5) cada uno. Que los grados de libertad de las t’s están entre 3 y 8 aproximadamente. Experimenta con otros tamaños muestrales y comenta los resultados obtenidos (y los tiempos de ejecución). Nota: este problema está motivado por una aplicación real: el ajuste de distribuciones de pérdida en banca y seguros. Típicamente, se mezclan dos distribuciones, una para la cola de la distribución y otra para el cuerpo. Hay técnicas frecuentistas (p.e., EM) para resolver estos problemas. Pero me parecen menos naturales y menos flexibles que la ruta 100% bayesiana.

Imagina que tienes que generar (reitero: generar) datos compatibles con el siguiente modelo: Tienes n sujetos a los que se proporciona un remedio para dormir en distintas dosis (conocidas) en distintos días. El número adicional de horas que duerme cada sujeto es lineal con una pendiente que depende de la dosis (una serie de dosis fijas). Esa recta tiene un término independiente (el número de horas que duerme el sujeto con una dosis igual a cero del remedio). Argumento que para generar los términos independientes usarías algo así como una normal de media igual a 8 horas. Seguro que usarías alguna otra distribución razonable para las pendientes (p.e., que prohibiese que con dosis pequeñas se durmiese, p.e., 80 horas). ...

A primeros de julio (de 2018) impartiré un curso de 15 horas de estadística bayesiana aplicada con Stan en la UPC (Barcelona). La información relevante está aquí y aquí. El proyecto y su definición es un tanto contradictorio en sus propios términos, lo reconozco. Es muy difícil hacer algo aplicado y, a la vez, bayesiano. Y más, con Stan. Además, podrían acusarme de hipócrita: ¿cuándo fue la última vez que facturé (recuérdese: facturable es el grado máximo de aplicado) por algo hecho con Stan? Porque la idea, en el fondo, es otra: esencialmente, cómo replantear modelos y estrategias de modelización, aunque se implenten con herramientas métodos de índole frecuentista, para enriquecerlos con la visión bayesiana. ...

ABC significa, entre otras cosas, approximate bayesian computation. Por lo que parece, consiste en calcular $P(\theta ,|, \text{datos})$ por el tradicional y directo método del rechazo. Es decir: Planteas un modelo generativo, con sus prioris y todo. Simulas casos, casos y casos. Te quedas con los que cumplen un criterio de aceptación. La distribución empírica de los parámetros en el subconjunto de los casos aceptados representa, en los libros está escrito, la distribución a posteriori. Sin MCMC ni historias. ...

Se ve que hay arqueólogos bayesianos. Un problema con el que se encuentran es que tropiezan con cacharros antiguos y quieren estimar su antigüedad. Así que prueban distintos métodos (¿químicos?), cada uno de los cuales con su precisión, y acaban recopilando una serie de estimaciones y errores. Obviamente, tienen que combinarlas de alguna manera. El modelo más simple es $$ M_i \sim N(\mu, \sigma_i)$$ donde $\mu$ es la antigüedad (desconocida) del artefacto y los $\sigma_i$ son las varianzas distintas de los distintos métodos de medida, que arrojan las estimaciones $M_i$. ...

Allá por el 2007 publicó The Independent una portada en que se retractaba. El diario había sido un histórico defensor de la legalización de la marihuana. Ese día hizo público su cambio de postura. Al parecer, motivada por las evidencias sobre los efectos sobre la salud mental. Este fin de semana he asistido a una serie de conferencias. En una de ellas participaba el representante de una organización que: Adoptaba de partida una posición militante, de parte, en cierto asunto de interés público. Se definía como data driven, evidence driven, etc. La pregunta obvia y que no tuve ocasión de plantear (por eso la traigo aquí) es la siguiente: si los datos y la evidencia se obstinaran en subrayar la bondad de la posición contraria a la que actualmente mantienen, ¿cuál de sus dos principios abandonarían primero? ...