##hyperparameter 
mu0<-1.9;//prior mean for theta, normal distribution
t20<-0.95^2;//prior variance for theta, normal distribution
s20<-0.01;//prior a for gamma distribution.
nu0<-1;//prior b for gamma distribution

##data
y<-c(1.64,1.70,1.72,1.74,1.82,1.82,1.82,1.90,2.08);
mean.y<-mean(y);//sample mean
var.y<-var(y);//sample variance
n<-length(y);//sample size

##initial values
T<-1000;
PHI<-matrix(nrow=T,ncol=2);
PHI[1,]<-phi<-c(mean.y,var.y);

###Gibbs Sampling
set.seed(3)
for(s in 2:T)
{

#sample a new theta value from its full conditional distribution(normal)
mu.new<-(mu0/t20 + n*mean.y*phi[2])/(1/t20 + n*phi[2])
t2.new<-1/(1/t20+n*phi[2])
phi[1]<-rnorm(1,mu.new,sqrt(t2.new))

#genereate a new precision value from full conditional distribution(gamma)
nu.new<-nu0+n
s2.new<-(nu0*s20+(n-1)*var.y + n*(mean.y-phi[1])^2)/nu.new;
phi[2]<-rgamma(1,nu.new/2,nu.new*s2.new/2)

#populate the parameter sets
PHI[s,]<-phi

}