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!
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
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.
Curioso si pintas.
«`
noabc <- params [!validos]
qqplot(res$lambda, noabc]
« `
¿Qué esperabas que saliera?