#source("c:/Work07/ShortCourse07/Programs/NormalHierarchy-2.txt",print.eval=TRUE)#
#Gibbs sampler for Girl's energy intake#
#Flat prior on theta
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))
a<-1;b<-1
#--------------------------------------------------------
for(i in 2:nsim)
{
mt<-mean(x)
vt<-sigma2[i-1]/n
theta[i]<-rnorm(1, mean=mt,sd=sqrt(vt))
sh<-(n/2)+a
sc<-1/((sum((x-theta[i])^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")