Gibbs sampler examples in R How can I implement Gibbs sampler for the posterior distribution, and estimating the marginal posterior distribution by making histogram? 
 A: Problem
Suppose $Y \sim \text{N}(\text{mean}  = \mu, \text{Var} = \frac{1}{\tau})$.
Based on a sample, obtain the posterior distributions of $\mu$ and $\tau$ using the Gibbs sampler.
Notation
$ \mu$ = population mean 
$ \tau$ = population precision (1 / variance)
$n$ = sample size
$\bar{y}$ = sample mean
$s^2$ = sample variance
Gibbs sampler
[Casella, G. & George, E. I. (1992). Explaining the Gibbs Sampler. The American Statistician, 46, 167–174.]

At iteration $i$ ($i = 1, \dots, N$):
  
  
*
  
*sample $\mu^{(i)}$ from $f(\mu \,|\, \tau^{(i - 1)}, \text{data})$   (see below)
  
*sample $\tau^{(i)}$ from $f(\tau \,|\, \mu^{(i)}, \text{data})$  (see below)
  

The theory ensures that after a sufficiently large number of iterations, $T$, the set $\{(\mu^{()},  \tau^{()} ) : i = T+1, \dots,  \}$ can be seen as a random sample from the joint posterior distribution.
Priors
$f(\mu, \tau) = f(\mu) \times f(\tau)$, with 
$f(\mu) \propto 1$
$f(\tau) \propto \tau^{-1}$
Conditional posterior for the mean, given the precision
$$(\mu \,|\, \tau, \text{data}) \sim \text{N}\Big(\bar{y}, \frac{1}{n\tau}\Big)$$
Conditional posterior for the precision, given the mean
$$(\tau \,|\, \mu, \text{data}) \sim \text{Gam}\Big(\frac{n}{2}, \frac{2}{(n-1)s^2 + n(\mu - \bar{y})^2} \Big)$$
(quick) R implementation
# summary statistics of sample
n    <- 30
ybar <- 15
s2   <- 3

# sample from the joint posterior (mu, tau | data)
mu     <- rep(NA, 11000)
tau    <- rep(NA, 11000)
T      <- 1000    # burnin
tau[1] <- 1  # initialisation
for(i in 2:11000) {   
    mu[i]  <- rnorm(n = 1, mean = ybar, sd = sqrt(1 / (n * tau[i - 1])))    
    tau[i] <- rgamma(n = 1, shape = n / 2, scale = 2 / ((n - 1) * s2 + n * (mu[i] - ybar)^2))
}
mu  <- mu[-(1:T)]   # remove burnin
tau <- tau[-(1:T)] # remove burnin

$$
$$
hist(mu)
hist(tau)


