# Obtener una muestra de p(theta | mu, tau, y)
conditional.theta=function(ybar,mu,tau,sigma){
theta=rep(0,nschools)
theta.hat=rep(0,nschools)
V.hat=rep(0,nschools)
for(j in 1:nschools){
V.hat[j]=1/(1/sigma[j]^2+1/(tau^2))
theta.hat[j]=(ybar[j]/sigma[j]^2+mu/tau^2)*V.hat[j]
theta[j]=rnorm(1,theta.hat[j],sqrt(V.hat[j]))
}
theta
}
# Obtener muestras de  p(mu | tau, y) de tama~no nsample
sample.mar.mu=function(ybar, tau, sigma,nsample)
{
V.mu.inv=sum(1/(sigma^2+tau^2))
mu.hat=sum((1/(sigma^2+tau^2))*ybar)/V.mu.inv
mu.sample=rnorm(nsample,mu.hat,sqrt(1/V.mu.inv))
mu.sample
}
# Evaluar p(tau | y)
marginal.tau=function(ybar,tau,sigma)
{
V.mu.inv=sum(1/(sigma^2+tau^2))
mu.hat=sum((1/(sigma^2+tau^2))*ybar)/V.mu.inv
eval=exp(-(ybar-mu.hat)^2/(2*(sigma^2+tau^2)))
eval=eval/sqrt(sigma^2+tau^2)
eval=sqrt(1/V.mu.inv)*prod(eval)
eval
}
########### Programa Principal ########################
# Lectura de los datos del archivo sa.scores
# School Treat.effect sd.effect
# A 28.39 14.9
# B  7.94 10.2 ...

sat.scores=read.table('sat.scores',header=TRUE)
ybar=sat.scores$Treat.effect
nschools=length(ybar)
sigma=sat.scores$sd.effect
# Grid to evaluate p(tau |y)
x.tau=seq(0.00001,40,length=1000)
# evaluar p(tau |y) en  1000 puntos en el
#intervalo [0.00001,40]
post.tau=apply(t(x.tau),2,marginal.tau,ybar=ybar,
sigma=sigma)
# simular 200 muestras de p(tau |y)
sample.tau=sample(x.tau,200,replace=TRUE,
prob=post.tau)
# simular 200 muestras de p(mu | tau, y)
sample.mu=apply(t(sample.tau),2,sample.mar.mu,
ybar=ybar, sigma=sigma,nsample=1)
# simular 200 muestras de p(theta | mu, tau, y)
sample.theta=matrix(0,ncol=nschools,nrow=200)
for (i in 1:200){
sample.theta[i,]=conditional.theta(ybar,
sample.mu[i], sample.tau[i],sigma)
}
# Medias esperadas a posteriori E(theta_j | tau, y)
# promediadas sobre mu
expected.theta=matrix(0,ncol=nschools,nrow=30)
x.tau.2=seq(0.00001,30,length=30)
for (i in 1:30){
sample.mu=sample.mar.mu(ybar,x.tau.2[i],sigma,
nsample=5000)
sample.theta.2=matrix(0,ncol=nschools,nrow=5000)
for (j in 1:5000){
sample.theta.2[j,]=conditional.theta(ybar,
sample.mu[j], x.tau.2[i],sigma)
}
expected.theta[i,]=apply(sample.theta.2,2,mean)
}
plot(x.tau,post.tau,type='l',ylab="",xlab="tau")
title('p(tau|y)')
dev.print(device=postscript,file="eje5.5.1.ps")
hist(sample.tau,ylab="",xlab="tau",main="p(tau|y)")

# 95\% P.I. for tau
sort(sample.tau)[5]
sort(sample.tau)[195]
for (i in 1:8){
a=sort(sample.theta[,i])[5]
b=sort(sample.theta[,i])[195]
c=mean(sample.theta[,i])
print(c(a,b,c))}
plot(x.tau.2,expected.theta[,1],ylim=c(-5,30),ylab="Estimated treatment effects",
xlab="tau",type='l')
text(x=11.64725,y=16.86634,"A",col=1)
lines(x.tau.2,expected.theta[,2],col=2)
text(x=11.89128,y=8.757426,"B",col=2)
lines(x.tau.2,expected.theta[,3],col=4)
text(x=12.98943,y=4.54703,"C",col=4)
lines(x.tau.2,expected.theta[,4],col=5)
text(x=10.30,y=7.0,"D",col=5)
lines(x.tau.2,expected.theta[,5],col=6)
text(x=9.93,y=3.37,"E",col=6)
lines(x.tau.2,expected.theta[,6],col="purple")
text(x=23,y=2.7,"F",col="purple")
lines(x.tau.2,expected.theta[,7],col="lightblue")
text(x=22,y=17.17,"G",col="lightblue")
lines(x.tau.2,expected.theta[,8],col="green")
text(x=22,y=12.17,"H",col="green")
dev.print(device=postscript,file="eje5.5.2.ps")