n=30; % sample size
 %Random Number Seed
s = RandStream('mt19937ar','Seed',1);
RandStream.setGlobalStream(s);

 y = 4 * randn(1,n) + 1;
%------------------------------------------
NN = 10000;
mus = []; taus = [];
sumdata = sum(y);
mu = 0;  % set the parameters as prior means
tau = 2; %
for i = 1 : NN
  newmu  = sqrt(1/(1+n*tau)) * randn + (tau * sumdata)/(1+n*tau);
  par   = 1+1/2 * sum ( (y - mu).^2);
  newtau = gamrnd(2 + n/2, 1/par); %par is rate
  mus = [mus newmu];
  taus = [taus newtau];
  mu=newmu;
  tau=newtau;
end
mean(mus)
1/mean(taus)