**A**s supplementary material to the ABC paper we just arXived, here is the R code I used to produce the Bayes factor comparisons between summary statistics in the normal versus Laplace example. *(***Warning:** running the R code takes a while!)

# ABC model comparison between Laplace and normal
nobs=10^4
nsims=100
Niter=10^5
sqrtwo=sqrt(2)
probA=probB=matrix(0,nsims,3)
dista=distb=rep(0,Niter)
pro=c(.001,.01,.1)
#A) Simulation from the normal model
for (sims in 1:nsims){
tru=rnorm(nobs)
#stat=c(mean(tru),median(tru),var(tru))
#stat=c(mean(tru^4),mean(tru^6))
stat=mad(tru)
mu=rnorm(Niter,sd=2)
for (t in 1:Niter){
#a) normal predictive
prop=rnorm(nobs,mean=mu[t])
#pstat=c(mean(prop),median(prop),var(prop))
#pstat=c(mean(prop^4),mean(prop^6))
pstat=mad(prop)
dista[t]=sum((pstat-stat)^2)
#b) Laplace predictive
prop=mu[t]+sample(c(-1,1),nobs,rep=TRUE)*rexp(nobs,rate=sqrtwo)
#pstat=c(mean(prop),median(prop),var(prop))
#pstat=c(mean(prop^4),mean(prop^6))
pstat=mad(prop)
distb[t]=sum((pstat-stat)^2)
}
epsi=quantile(c(dista,distb),prob=pro)
for (i in 1:3)
probA[sims,i]=sum(dista<epsi[i])/(2*Niter*pro[i])
}
#B) Simulation from the Laplace model
for (sims in 1:nsims){
tru=sample(c(-1,1),nobs,rep=TRUE)*rexp(nobs,rate=sqrtwo)
#stat=c(mean(tru),median(tru),var(tru))
stat=mad(tru)
mu=rnorm(Niter,sd=2)
for (t in 1:Niter){
#a) normal predictive
prop=rnorm(nobs,mean=mu[t])
#pstat=c(mean(prop),median(prop),var(prop))
#pstat=c(mean(prop^4),mean(prop^6))
pstat=mad(prop)
dista[t]=sum((pstat-stat)^2)
#b) Laplace predictive
prop=mu[t]+sample(c(-1,1),nobs,rep=TRUE)*rexp(nobs,rate=sqrtwo)
#pstat=c(mean(prop),median(prop),var(prop))
#pstat=c(mean(prop^4),mean(prop^6))
pstat=mad(prop)
distb[t]=sum((pstat-stat)^2)
}
epsi=quantile(c(dista,distb),prob=pro)
for (i in 1:3)
probB[sims,i]=sum(dista<epsi[i])/(2*Niter*pro[i])
}

Filed under: R, Statistics, University life Tagged: ABC, Bayesian model choice, Laplace distribution, R, summary statistics

*Related*

To

**leave a comment** for the author, please follow the link and comment on his blog:

** Xi'an's Og » R**.

R-bloggers.com offers

**daily e-mail updates** about

R news and

tutorials on topics such as: visualization (

ggplot2,

Boxplots,

maps,

animation), programming (

RStudio,

Sweave,

LaTeX,

SQL,

Eclipse,

git,

hadoop,

Web Scraping) statistics (

regression,

PCA,

time series,

trading) and more...

If you got this far, why not

__subscribe for updates__ from the site? Choose your flavor:

e-mail,

twitter,

RSS, or

facebook...

**Tags:** ABC, Bayesian model choice, Laplace distribution, R, statistics, summary statistics, University life