Stan

Esta entrada muestra cómo afrontar (con Stan) un problema que encontré el otro día en un lugar que no puedo mencionar pero en el que sé que me leen (y los destinatarios sabrán que va por ellos). El contexto es el siguiente: se hace un test A/B donde la variable de interés son unos conteos. Hay varios grupos (aquí los reduciré a dos) y los datos siguen aproximadamente (aquí omitiré la parte de la inflación de ceros) una distribución de Poisson. Pero solo aproximadamente: existe sobredispersión, es decir, la varianza de los datos excede su media. ...

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")])

Hay una discusión imperdible sobre datos anchos y largos (y modelos con y sin estructura) a partir del minuto 4:20 de https://www.youtube.com/watch?v=VnNdhsm0rJQ El resto también vale muchísimo la pena.

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. ...

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.

Breiman habló de las dos. Dice, y tiene razón, que: Según él, la estadística tradicional rellena la caja negra con: ¡Aburrido, aburrido, aburrido! Aburrido y limitado (aunque, hay que admitirlo, útil en ocasiones muy concretas). Breiman sugiere sustituir las cajas negras que encontramos en la naturaleza por otras cajas negras conceptuales: Que es aún más aburrido y patrimonio, además, de toda suerte de script kiddies. La tercera cultura reemplaza la caja negra por un modelo generativo que simula el comportamiento de la naturaleza (i.e., del sistema generador de números aleatorios pero con estructura). Y usa Stan (o sus alternativas) para estimar, predecir y, en última instancia, facilitar decisiones informadas.

Este mes de julio, cuórum mediante, impartiré en la UPC un curso que he maltitulado, mor de brevedad, Estadística Bayesiana Aplicada. Los cursos de estadística bayesiana son teoría, mucha teoría, y unos ejemplos tontos que quieren justificarla. Del tipo: hagamos lo que ya sabemos hacer de otra manera más; busquemos una alternativa molona al p-valor (y usémosla como usar íamos un p-valor, por supuesto), etc. Mi curso debería haberse titulado algo así como: Problemas reales (aunque simplificados por motivos estrictamente pedagógicos) resueltos con tecnología bayesiana porque, si no, dígame Vd. cómo lo haría: ¿con optim? Jajajajaja… ...

Pues no se sabe bien. Además, habrá quién pudiéndola haber averiguado, prefirió dejarse llevar por la intuición y errar. Pero volvamos a los hechos. Dado En un país hipotético, las familias tienen críos hasta que nace el primer varón. En un año, en promedio, nacen: — Carlos Gil Bellosta (@gilbellosta) December 10, 2017 la pregunta urgente es: ¿cuántos podrían haber conocido la respuesta? Suponiendo que el conocimiento de la respuesta es algo binarizable (¿lo es?), la distribución del número de respuestas correctas sería $pN + X$, donde $N$ es el número total de respuestas, $p$ es la proporción de quienes sabe la respuesta y $X \sim B(N - pN, 1/3)$, suponiendo siempre que $pN$ es entero. ...

Envalentonado por el comentario de Iñaki Úcar a mi entrada del otro día, que me remitía a este artículo, decidí rizar el rizo y crear intervalos de confianza no ya discontinuos sino con otra propiedad topológica imposible: homeomorfos con un toro. Y aquí está: El modelo, el código y demás, library(rstan) library(ggplot2) n <- 100 a1 <- 1 a2 <- 1 sigma <- 0.4 datos <- data.frame(x1 = rnorm(n, 2, 0.1), x2 = rnorm(n, 2, 0.1)) datos$y <- a1^datos$x1 + a2^datos$x2 + rnorm(n, 0, sigma) codigo <- " data { int<lower=1> N; real y[N]; real x1[N]; real x2[N]; } parameters { real<lower=-3, upper="3"> a1; real<lower=-3, upper="3"> a2; real<lower=0, upper="3"> sigma; } model { for (n in 1:N) y[n] ~ normal(fabs(a1)^x1[n] + fabs(a2)^x2[n], sigma); }" fit <- stan(model_code = codigo, data = list(N = length(datos$y), y = datos$y, x1 = datos$x1, x2 = datos$x2), iter=40000, warmup=2000, chains=1, thin=10) res <- as.data.frame(fit) ggplot(res, aes(x = a1, y = a2)) + geom_point(alpha = 0.1) De nuevo, no son intervalos propiamente dichos, lo convengo. Pero son configuraciones más fieles al espíritu de lo que un intervalo de confianza es y representa que su(s) letra(s) I N T E R V A L O. ...

Un año, el 2016, mueren 1160 personas en accidentes de tráfico. El anterior, 1131, i.e., 29 menos. Ruido estadístico aparte, ¿aumentan? Comenzamos a optar. Primera elección subjetiva: son muestras de una Poisson de parámetro desconocido. La pregunta: ¿el mismo? Una manera de estudiar lo anterior es plantear 1160 ~ poisson(lambda * (1 + incr)) 1131 ~ poisson(lambda) y estudiar la distribución de incr. Que a saber qué distribución tendrá (teóricamente). Pero, ¿importa? Mejor que rebuscar a ver qué distribución podría tener la cosa, basta con envolverlo en un poco de seudo-C++, ...