ABC in Roma [R lab #1]
[This article was first published on Xi'an's Og » R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Here are the R codes of the R labs organised by Serena Arima in supplement of my lectures. This is quite impressive and helpful to the students, as illustrated by the first example below (using the abc software).
### Example 1: Conjugate model (Normal-Inverse Gamma) ### y1,y2,...,yn ~N(mu,sigma2) ### mu|sigma2 ~ N(0,sigma2), sigma2 ~IG(1/2,1/2) library(abc) ### Iris data: sepal width of Iris Setosa data(iris3) y=iris3[,2,1] ### We want to obtain the following quantities ### "par.sim" "post.mu" "post.sigma2" "stat.obs" "stat.sim" ## STAT.OBS: stat.obs are mean and variance (log scale) of the data mean(y) log(var(y)) stat.obs stat.oss=c(mean(y),log(var(y))) ### PAR.SIM: par.sim simulated values from the prior distribution n.sim=10000 gdl=1 invsigma.sim=rchisq(n.sim,df=gdl) sigma.sim=gdl/invsigma.sim m=3 mu.sim=c() for(i in 1:length(sigma.sim)){ mu.sim=c(mu.sim,rnorm(1,mean=m,sd=sqrt(((sigma.sim[i])))))} prior.sim=data.frame(mu.sim,sigma.sim) ### STAT.SIM: for mu and sigma simulated from the prior, ### generate data from the model y ~ N(mu,sigma^2) stat.simulated=matrix(NA,nrow=length(mu.sim),ncol=2) for(i in 1:length(mu.sim)){ y.new=rnorm(length(y),mu.sim[i],sqrt(sigma.sim[i])) stat.simulated[i,1]=mean(y.new) stat.simulated[i,2]=log(var(y.new))} ### Obtain posterior distribution using ABC post.value=abc(target=stat.oss, param=prior.sim, sumstat=data.frame(stat.simulated),tol=0.001,method="rejection") summary(post.value) posterior.values=post.value$unadj.values mu.post=posterior.values[,1] sigma.post=posterior.values[,2] ### True values, thanks to conjugancy post.mean.mu=(length(y)/(length(y)+1))*mean(y) post.a.sigma=length(y)/2 post.b.sigma=0.5+0.5*sum((y-mean(y))^2) hist(mu.post,main="Posterior distribution of mu") abline(v=post.mean.mu,col=2,lwd=2) hist(sigma.post,main="Posterior distribution of sigma2") abline(v=post.b.sigma/(post.a.sigma-1),col=2,lwd=2)
I am having a great time teaching this “ABC in Roma” course, in particular because of the level of interaction and exchange with the participants (after, if not during, the classes).
Filed under: R, Statistics, University life Tagged: ABC, La Sapienza, PhD course, R, Roma
To leave a comment for the author, please follow the link and comment on their blog: Xi'an's Og » R.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.