# set.seed(11);

alpha1 <- 3;
beta1<- 2;
mu1 <- 5;
lambda1 <- 10;

numIter <- 5;

tau <- rep(1,numIter);
mu <- rep(0,numIter);

for (i in 1:numIter) {
   mu[2*i] <- rnorm(1,mu1,sqrt(1/(tau[(2*i)-1]*lambda1)));
   tau[2*i] <- tau[(2*i)-1];
   mu[(2*i)+1] <- mu[2*i]
   tau[(2*i)+1] <- rgamma(1,alpha1+0.5, beta1+((lambda1/2)*((mu[2*i]-mu1)^2)));
}

plot(mu,tau,type="o");
cat(mean(tau^3),"\n");