x<-rnorm(100,3,sqrt(12))

I<-500

sumx<-sum(x)
m0<-0
s0<-100
a<-1
b<-0.001
n<-length(x)

mu<-tau<-rep(NA,I+1)
tau[1]<-1000
mu[1]<-1

for (i in (1:I)+1)
{
  tau[i]<-rgamma(1,a+n/2,b+0.5*sum((x-mu[i-1])^2))
  m<-tau[i]*sumx+m0/s0
  s<-n*tau[i]+1/s0
  mu[i]<-rnorm(1,m/s,1/s)
}

tau<-tau[-1]
mu<-mu[-1]

plot(mu,type="l")
plot(1/tau,type="l")
plot(tau,type="l")
plot(mu,tau)
plot(mu[1:10],(1/tau)[1:10],type="s")
plot(density(mu))
plot(density(1/tau))