#source("c:/Work07/ShortCourse07/Programs/NormalHierarchy-1.txt",print.eval=TRUE)#
#Gibbs sampler for Girl's energy intake#
x<-c(91,504,557,609,693,727,764,803,857,929,970,1043,1089,1195,1384,1713)
x<-log(x)		#take logs
n<-length(x)
#---------------Initial values---------------------------
nsim<-5000
theta<-array(mean(x), c(nsim,1))
sigma2<-array(var(x),c(nsim,1))
theta0<-5
tau2<-10
a<-1;b<-1
#--------------------------------------------------------
for(i in 2:nsim)
{
mt<-(theta0/(1+n*tau2))+(n*tau2*mean(x)/(1+n*tau2))
vt<-sigma2[i-1]*tau2/(1+n*tau2)
theta[i]<-rnorm(1, mean=mt,sd=sqrt(vt))
sh<-(n+1)/2+a
sc<-1/(sum((x-theta[i])^2)/2+(theta[i]-theta0)^2/2+(1/b))
sigma2[i]<-1/rgamma(1,shape=sh,scale=sc)
}
#----------------Plot---------------------------------
par(mfrow=c(1,2))
hist(theta, col="grey")
hist(sigma2, col="grey")