ABC

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, y un mecanismo de generación de datos , donde es un vector de parámetros.

Supongamos que tenemos una estimación puntual de de y queremos construir una posteriori. La idea es obtener valores próximos a y generar muestras de acuerdo con . Entonces se consideran buenas las simulaciones tales que para una cierta distancia y un valor dado de . Los valores 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, , es:

stan_code <- "
  data {
    int suma;
    int n;
  }

  parameters {
    real 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!