# Gibbs sampling for the model:
#
#    x1, x2, ..., xn ~ N(mu,1/tau)
#  
#    mu  ~ N(0,1)
#    tau ~ Exp(1)
#
# The function below takes the data as its first argument, followed by initial
# values for mu and tau, and the number of iterations to do.  It returns a list
# containing vectors of values for mu and tau at each iteration.
#
# The function could be speeded up by pre-computing some sufficient statistics,
# but this isn't done here for simplicity.

gibbs <- function (x,mu0,tau0,k)
{
  n <- length(x)

  mu.vec  <- numeric(k)
  tau.vec <- numeric(k)

  mu  <- mu0
  tau <- tau0

  for (i in 1:k)
  { 
    mu  <- rnorm (1, (n*tau*mean(x))/(1+n*tau), 1/sqrt(1+n*tau))
    tau <- rgamma (1, 1+n/2) / (1 + sum((x-mu)^2)/2)

    mu.vec[i]  <- mu
    tau.vec[i] <- tau
  }

  list (mu=mu.vec, tau=tau.vec)
}