ndraws <- 10000

mu0 <- 0
eta0 <- 1
alpha0 <- 2
beta0 <- 3

data <- scan()
   2.7785780  3.1950704  1.2179930  3.6112721  2.3768019  0.5801824
   1.6884122  0.9132051  2.1997683  0.7394964  2.3064339 -0.6005665
   3.6549066  0.8674690  1.6731264  0.8864085  3.7975642  3.1209520
   1.8015749  2.3796327

draws <- data.frame(matrix(NA,nrow=ndraws,ncol=2))
names(draws) <- c("mu","lambda")

draws[1,"mu"] <- mu0
draws[1,"lambda"] <- alpha0/beta0
for ( i in 2:ndraws ) {
  pre <- eta0 + length(data)*draws[i-1,"lambda"]
  mea <- ( eta0*mu0 + draws[i-1,"lambda"]*sum(data) ) / pre
  draws[i,"mu"] <- rnorm(1,mean=mea,sd=1/sqrt(pre))
  alp <- alpha0 + 0.5*length(data)
  bet <- beta0 + 0.5*sum((data-draws[i,"mu"])^2)
  draws[i,"lambda"] <- rgamma(1,shape=alp,rate=bet)
}

mean(draws$mu)                           # Point estimate of mu
quantile(draws$mu,c(0.025,0.975))        # 95% credible interval for mu

mean(draws$lambda)                       # Point estimate of lambda
quantile(draws$lambda,c(0.025,0.975))    # 95% credible interval for lambda

hist(draws$mu)
hist(draws$lambda)
plot(draws$mu,draws$lambda)
acf(draws$mu)
acf(draws$lambda)