Here you will find daily news and tutorials about R, contributed by over 573 bloggers.
There are many ways to follow us - By e-mail:On Facebook: If you are an R blogger yourself you are invited to add your own R content feed to this site (Non-English R bloggers should add themselves- here)

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])
}