Archive for February, 2012

ABC in Roma [R lab #1]

Posted in R, Statistics, University life with tags , , , , on February 29, 2012 by xi'an

Here are the R codes of the first R lab written and ran 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).

A Roma [3]

Posted in pictures, Travel on February 29, 2012 by xi'an

Elsevier’s letter to mathematicians

Posted in Books, University life on February 29, 2012 by xi'an

Just like Timothy Gower and thousands of others, I have received in the past days an email from Elsevier about their pricing, agressive marketing policy towards University libraries, and support of the Research Work Act. This email was an answer to the growing action of the mathematical community to boycott Elsevier through [not] publishing/refereeing/aediting in Elsevier’s journals, an action started by an earlier blog of Timothy Gower. He has now posted his extensive reaction to this email and, as it perfectly fits mine, I see no point in writing another post (esp. with very limited posting abilities this week in Roma…)

winter trees (3)

Posted in pictures, Travel, University life with tags , , , , on February 28, 2012 by xi'an

A Roma [2]

Posted in pictures, Travel on February 28, 2012 by xi'an

Follow

Get every new post delivered to your Inbox.

Join 598 other followers