4
$\begingroup$

To understand how ABC works, I like to use it in a simple linear regression model. I am using EasyABC package in R. My problem is that the results I get are not good. In particular, the posterior distributions are extremely wide. I think this is partly because I am using prior distributions with a very large variance. In fact, if I use different priors (e.g., Unif(0,20) --- something much tighter), results are better. But I am using those priors to mean something commonly used as non-informative priors in Bayesian method.

My question is that is there anything wrong in the way ABC is used in this example (the code is given below)? Why is the posterior so bad?

library(EasyABC)

## data generation
n <- 100
intercept <- 10
slope     <- 1
s         <- 2.5
x <- runif(n,10,20)
y <- rnorm(n,intercept+slope*x,s)

model <- lm(y~x)
intercept.data <- coef(model)[1]
slope.data <- coef(model)[2]
sd.data <- summary(model)$sigma

sum_stat_obs <- c(intercept.data,slope.data,sd.data)

sim.abc <- function(par){ 
   samples <- rnorm(n, mean =par[1]+par[2]*x, sd = par[3])
   model0 <- lm(samples~x)
   intercept.sim <- coef(model0)[1]
   slope.sim <- coef(model0)[2]
   sd.sim <- summary(model0)$sigma
  return(c(intercept.sim,slope.sim,sd.sim))
}

model.abc <-ABC_mcmc(method="Marjoram", model=sim.abc,
  prior=list(c("normal",0,1000),c("normal",0,1000),c("exponential",1)), 
  summary_stat_target=sum_stat_obs, n_rec = 10000)

par(mfrow=c(2,2))
plot(model.abc$param[5000:10000,1],type="l")
hist(model.abc$param[5000:10000,1], main = "Posterior for intercept")
plot(model.abc$param[5000:10000,2],type="l")
hist(model.abc$param[5000:10000,2], main = "Posterior for slope")

enter image description here

$\endgroup$

1 Answer 1

3
$\begingroup$

The main problem is the tolerance together with the terrible choice of priors (recall that ABC does not work as other MCMC methods since it is based on simulations from the priors). You need to reduce the tolerance_quantile in order to reduce the distance between simulated and observed summary statistics.

See for instance (better tolerance):

model.abc <-ABC_mcmc(method="Marjoram", model=sim.abc,
                     prior=list(c("normal",0,100),c("normal",0,100),c("exponential",1)), 
                     summary_stat_target=sum_stat_obs, n_rec = 10000, tolerance_quantile=0.001)

apply(model.abc$param[5000:10000,],2,median)
apply(model.abc$param[5000:10000,],2,mean)

and (better tolerance and priors):

model.abc <-ABC_mcmc(method="Marjoram", model=sim.abc,
                     prior=list(c("normal",10,5),c("normal",1,3),c("exponential",2)), 
                     summary_stat_target=sum_stat_obs, n_rec = 10000, tolerance_quantile=0.0005)

apply(model.abc$param[5000:10000,],2,median)
apply(model.abc$param[5000:10000,],2,mean)

par(mfrow=c(3,2))
plot(model.abc$param[5000:10000,1],type="l")
hist(model.abc$param[5000:10000,1], main = "Posterior for intercept")
plot(model.abc$param[5000:10000,2],type="l")
hist(model.abc$param[5000:10000,2], main = "Posterior for slope")
plot(model.abc$param[5000:10000,3],type="l")
hist(model.abc$param[5000:10000,3], main = "Posterior for sigma")
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.