#source("c:/Work07/ShortCourse07/Programs/NormalHierarchy-3.txt",print.eval=TRUE)#
#Gibbs sampler for Girl's  and Boy's energy intake#
y1<-c(91,504,557,609,693,727,764,803,857,929,970,1043,1089,1195,1384,1713)
y2<-c(457,645,736,790,899,991,1104,1154,1203,1320,1417,1569,1685,1843,2296,2710)
y1<-log(y1);y2<-log(y2)		#take logs
n1<-length(y1);n2<-length(y2)
SSE<-sum((y1-mean(y1))^2)+ sum((y2-mean(y2))^2)
k=2				#number of treatments
#---------------Initial values---------------------------
nsim<-5000
mu<-array(mean(y), c(nsim,1))
alpha1<-array(0, c(nsim,1));alpha2<-array(0, c(nsim,1))
sigma2<-array((var(y1)+var(y2))/2,c(nsim,1))
tau2<-array((var(y1)+var(y2))/2,c(nsim,1))
a<-2;b<-2
#--------------------------------------------------------
for(i in 2:nsim)
{
#----------------------mu----------------------------------
m<-mean(c(y1,y2))-mean(c(alpha1[i-1],alpha2[i-1]))
v<-sigma2[i-1]/(n1+n2)
mu[i]<-rnorm(1, mean=m,sd=sqrt(v))
#---------------------alpha--------------------------------
v1<-sigma2[i-1]*tau2[i-1]/(sigma2[i-1]+n1*tau2[i-1])
alpha1[i]<-rnorm(1,mean=n1*v1*(mean(y1)-mu[i]), sd=sqrt(v1))
v2<-sigma2[i-1]*tau2[i-1]/(sigma2[i-1]+n2*tau2[i-1])
alpha2[i]<-rnorm(1,mean=n2*v2*(mean(y2)-mu[i]), sd=sqrt(v2))
#---------------------sigma---------------------------------
sh1<-((n1+n2)/2)+a
sc1<-1/(.5*(sum((y1-alpha1[i]-mu[i])^2)+sum((y2-alpha2[i]-mu[i])^2))+1/b)
sigma2[i]<-1/rgamma(1,shape=sh1,scale=sc1)
#---------------------tau------------------------------------
sh2<-(k/2)+a
sc2<-1/(.5*alpha1[i]^2+.5*alpha2[i]^2+1/b)
tau2[i]<-1/rgamma(1,shape=sh2,scale=sc2)
}
#----------------Plot---------------------------------
par(mfrow=c(1,3))
hist(mu, col="grey")
hist(tau2, col="grey",xlim=c(0,1),breaks=50)
hist(sigma2, col="grey")