ABC (I)

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.

Supongamos que tenemos una estimación puntual de $latex \hat{\theta}$ de $latex \theta$ y queremos construir una posteriori. La idea es obtener valores $latex \theta_i$ próximos a $latex \hat{\theta}$ y generar muestras $latex X_i$ de acuerdo con $latex f(\theta_i)$. Entonces se consideran buenas las simulaciones $latex X_i$ tales que $latex d(X, X_i) \le \epsilon$ para una cierta distancia $latex d$ y un valor dado de $latex \epsilon$. Los valores $latex \theta_i$ que pasan el filtro deberían ser, aproximadamente, una muestra de la distribución a posteriori.

Véamoslo en acción. Primero, creamos un mecanismo generador de datos:

set.seed(155)
my_lambda <- 7
n_trials <- 10
generador <- function(lambda) mean(rpois(n_trials, lambda))

mis_datos <- generador(my_lambda)

La distribución a posteriori de nuestro valor de interés, $latex \lambda$, es:

stan_code <- "
data {
    int<lower = 0> suma;
    int<lower = 0> n;
}

parameters {
    real<lower = 0, upper = 20> lambda;
}

model {
    suma ~ poisson(n * lambda);
}
"

library(rstan)
fit <- stan(model_code = stan_code,
            data = list(suma = round(mis_datos * n_trials),
                        n = n_trials),
            chains = 1, iter = 12000,
            warmup = 2000, thin = 10)

res <- as.data.frame(fit)
hist(res$lambda, density = TRUE,
    breaks = 50, main = "distribución de lambda",
    xlab = "", ylab = "")

Es decir,

Veamos lo que podemos hacer con ABC:

params <- runif(1e6, 0, 20)
generados <- sapply(params, generador)
validos <- abs(generados - mis_datos) < 0.1
abc <- params[validos]

qqplot(res$lambda, abc, main = "stan vs abc")
abline(a = 0, b = 1, col = "red")

Tachán:

Pas mal!