Selecting statistics for ABC model choice [R code]

November 1, 2011
By

(This article was first published on Xi'an's Og » R, and kindly contributed to R-bloggers)

As 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

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: , , , , , ,

Comments are closed.