(This article was first published on

**Xi'an's Og » R**, and kindly contributed to R-bloggers)**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

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...