Estadística Bayesiana

Leyendo sobre si dizque PyTorch le siega la hierba debajo de los pies a TensorFlow, averigué la existencia de Pyro.

Pyro se autopresenta como Deep Universal Probabilistic Programming, pero aplicando métodos porfirianos (ya sabéis: género próximo y diferencia específica), es, o pretende ser, Stan en Python y a escala.

Aquí van mis dos primeras impresiones, basadas en una inspección superficial de los tutoriales.

En primer lugar, aunque Pyro permite usar (distintas versiones de) MCMC, parece que su especialidad es la inferencia variacional estocástica. Que parece funcionar de la siguiente manera. En el MCMC tradicional uno obtiene una muestra de la distribución (a posteriori, para los amigos) de los parámetros de interés. Eso es todo: vectores de puntos. En la inferencia variacional estocástica, uno preespecifica la forma paramétrica de la posteriori y el algoritmo calcula sus parámetros a partir de los valores simulados. Por ejemplo, uno va y dice: me da que la distribución del término independiente de mi regresión lineal va a ser normal. Entonces, Pyro responde: si es normal, la mejor media y desviación estándar que encuentro son tal y cual.

Rootclaim es un portal donde la gente plantea preguntas como

plantea hipótesis como

se recogen evidencias y usando este método (leedlo, es sumamente aprovechable: usa la palabra bayesian 23 veces), llega a conclusiones tales como

Muestreo. Se trata de seleccionar unas unidades experimentales (proceso caro) y tratar de estimar una proporción (p.e.) en la población total.

Existen técnicas para estimar el valor N mínimo para garantizar cierto margen de error. Pero dichas técnicas requieren conocer (algo d-) el resultado del experimento para estimar N (p.e. una estimación de la proporción que cabe esperar).

Circulus in demonstrando.

Bayes. Ve examinando unidades y actualiza tus intervalos de credibilidad hasta que tengan la anchura solicitada.

Aquí se describe una suerte de recíproco para el teorema de Bernstein–von Mises. Aquí se resume de esta manera:

The celebrated Aumann’s Agreement Theorem shows that two rational agents with the same priors on an event who make different observations will always converge on the same posteriors after some civilized conversation over tea.

En resumen:

• B-vM: frente a la misma evidencia, observadores con prioris distintas tienen posteriores similares.
• Aumann: frente a evidencias disímiles, observadores con las mismas prioris pueden acordar posterioris similares.

Las dimensiones altas son un campo minado para la intuición. Hace poco (y he perdido la referencia) leí a un matemático que trabajaba en problemas en dimensiones altas decir que le gustaba representar y pensar en las bolas (regiones del espacio a distancia <1 de 0) en esos espacios usando figuras cóncavas, como las que aparecen a la izquierda de

precisamente porque una de las propiedades más fructíferas de las bolas en altas dimensiones es que apenas tienen interior. De hecho, es trivial probar que la proporción del volumen de una bola a distancia mayor que $latex \epsilon$ de su borde tiende a cero con la dimensión.

[Esta entrada recoge la pregunta y la duda que motivó una conversación con Javier Nogales en Twitter hace unos días.]

Citaba (él) un resultado de Theobald de 1974 (¿tanto lleva ridge entre nosotros? ¡habría jurado que menos!) que viene a decir que siempre existe un peso $latex \lambda$ para el que ridge es mejor que OLS.

Ves el álgebra y piensas: verdad será.

Pero te fías de tu propia intuición y piensas: ¡vaya un resultado contraintuitivo si no contradictorio! Porque:

Para muchos, el futuro de la llamada ciencia de datos seguirá la estela dejada por

y sus continuadores usando cosas deep. Pero a la vez, sin tanto estruendo y con una mucho menor cobertura mediática, otros están trazando una ruta alternativa que ilustran artículos como Bayes and Big Data: The Consensus Monte Carlo Algorithm (atención todos a lo que hace uno de sus coautores, Steven L. Scott, que convierte en oro todo lo que toca). Como abrebocas, su resumen (con mi subrayado):

Esa frase la he pronunciado en alguna ocasión y no sé si la habré escrito en este blog. La reescribo porque hace apenas unas horas he leído un artículo en el que un tipo ha redescubierto el partial pooling (quien lo ignore lea esto urgentemente). Claro, proponía unas cosas tan raras como ocurrentes que se reducían en la estrategia que he contado: tengo cierta intuición de una idea genial que no llego a aprehender enteramente y procedo a moverme dando tumbos y a golpe de ocurrencias en la difusa dirección en la que parece apuntar.

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 $latex B(a,b)$ (beta). Tomas una muestra de valores $latex 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, $latex X$ y un mecanismo de generación de datos $latex X^\prime = f(\theta)$, donde $latex \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.