ABC

23 octubre, 2018 Carlos J. Gil Bellosta 5 comentarios

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.

Supongamos que tenemos una estimación puntual de $\hat{\theta}$ de $\theta$ y queremos construir una posteriori. La idea es obtener valores $\theta_i$ próximos a $\hat{\theta}$ y generar muestras $X_i$ de acuerdo con $f(\theta_i)$ . Entonces se consideran buenas las simulaciones $X_i$ tales que $d(X, X_i) \le \epsilon$ para una cierta distancia $d$ y un valor dado de $\epsilon$ . Los valores $\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)

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, $\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 = "")

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

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!

5 comentarios sobre “ABC”

Jose Luis Cañadas Reche 23 octubre, 2018 11:23

Pues parece un “truco” útil para no tener que esperar si Stan tarda eones. Aunque me sigue pareciendo algo similar al LRT. https://en.m.wikipedia.org/wiki/Likelihood-ratio_test
En LRT entre dos modelos comparo las verosimilitudes y rechazo uno de los modelos si su verosimilitud es menor que la otra en cierto grado. En ABC coges thetas próximas a la estimada y comparas si los datos generados por ese nuevo modelo generativo son próximos a los actuales en cierto grado. Me parecen conceptos similares vistos desde ópticas distintas
Jose Luis Cañadas Reche 23 octubre, 2018 11:29

Salvo por la diferencia que en la óptica frecuentista el parámetro se entiende como fijo y en la bayesiana el parámetro a estimar es una v.a.
Pingback: ABC (2) – datanalytics
Jose Luis Cañadas Reche 24 octubre, 2018 19:30

Curioso si pintas.
“`
noabc <- params [!validos]
qqplot(res$lambda, noabc]

“ `
Carlos J. Gil Bellosta 29 octubre, 2018 0:23

¿Qué esperabas que saliera?

Comenta Cancela la respuesta

Este sitio usa Akismet para reducir el spam. Aprende cómo se procesan los datos de tus comentarios.