Gibbs sampler with R

## Gibbs sampler

What’s a Gibbs Sampler

• It’s an algorithm for simulating Markov chains.

• Let $$\theta=(\theta_1,\ldots,\theta_k)$$ be a random vector with joint distribution $$f$$,

• Aim: simulate a sample with distribution $$f$$.

How it works?

1. Choose arbitrary starting values: $$\theta^{(0)}=(\theta_1^{(0)},\ldots,\theta_k^{(0)})$$

2. Sample new values for each element of $$\theta$$ by cycling through the following steps:

• Sample a new value for $$\theta_1^{(1)}$$ from the full conditional distribution of $$\theta_1$$ given the most recent values of all other elements of $$\theta$$: $\theta_1^{(1)}\sim p\left(\theta_1 \mid\theta_2^{(0)},\ldots,\theta_k^{(0)} \right)$

• Sample a new value $$\theta_2^{(1)}$$ for the 2nd component of $$\theta$$, from its full conditional distribution $\theta_2^{(1)}\sim p\left(\theta_2 \mid\theta_1^{(1)},\theta_3^{(0)}\ldots,\theta_k^{(0)} \right)$

• $$\theta_k^{(1)}\sim p\left(\theta_1 \mid\theta_1^{(1)},\theta_2^{(1)}\ldots,\theta_{k-1}^{(1)} \right)$$

This completes one iteration of the Gibbs sampler and generates a new realization of the vector of unknowns $$\theta^{(1)}$$.

1. Repeat Step 2 many times and obtain a sample $$\theta^{(1)},\ldots, \theta^{(T)}$$

To better understand the Gibbs Sampler, I suggest to read the paper by Casela and George 1992

## Markov Chain

• $$\theta^{(1)},\ldots, \theta^{(T)}$$ is a Markov Chain

• If $$\theta^{(1)},\ldots, \theta^{(T)}$$ converges, then from $$N$$, $$\theta^{(N)},\ldots, \theta^{(N)}$$ is a sample with $$f$$ distribution

• coda and codadiags two R packages for convergence diagnostics

## Example 1: Posterior of the couple mean and precision of a Normal distribution

Problem

Suppose that we had observed a sample $$y_1,\ldots,y_n$$ generated from random variable $$Y\sim\mathcal{N}\left(\mbox{mean}=\mu,\,\mbox{Var}=1/\tau\right)$$ where $$\mu$$ and $$\tau$$ are unknown. Suppose also that we had defined the following prior distribution $f(\mu,\,\tau) \propto 1/\tau$

Based on a sample, obtain the posterior distributions of $$\mu$$ and $$\tau$$ using the Gibbs sampler and write an R code simulating this posterior distribution.

Solution

We can first prouve that

• Conditional posterior for the mean, given the precision is given by

$\mu\mid y_1,\ldots,y_n,\tau \sim \mathcal{N}\left(\overline{y},\displaystyle\frac{1}{n\tau}\right)$

• Conditional posterior for the precision, given the mean

$\tau\mid y_1,\ldots,y_n,\mu \sim \mbox{Gamma}\left(\displaystyle\frac{n}{2},\displaystyle\frac{2}{(n-1)s^2+n(\mu-\overline{y})^2}\right)$ where $\overline{y}=\displaystyle\frac{1}{n}\sum_{i=1}^n y_i$ is the sample mean and $s^2=\displaystyle\frac{1}{n-1}\sum_{i=1}^n (y_i-\overline{y})^2$ is the sample variance.

Let’s then implement then the Gibbs Sampler into R. Here’s the code

> # The data
> n    <- 40
> ybar <- 35
> s2   <- 5
>
> # sample from the joint posterior (mu, tau | data)
> mu     <- rep(NA, 11000)
> tau    <- rep(NA, 11000)
> T      <- 1000    # burnin
> tau[1] <- 1  # initialisation
> mu[1] <- 1
> 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

Let’s quickly visualize the marginal distribution obtained from the sample runned above.

> summary(mu)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
33.66   34.75   35.00   35.00   35.24   36.49 
> hist(mu)

> summary(1/tau)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
2.364   4.353   5.061   5.241   5.934  13.035 
> hist(1/tau)

## Example 2: Bivariate Gaussian Distribution

Let $\theta=(\theta_1,\,\theta_2)\sim \mathcal{N}_2(0,I_2)$ with density function $f(\theta_1,\,\theta_2)=\displaystyle\frac{1}{\sqrt{2\pi}}\exp\left(-\displaystyle\frac{1}{2}(\theta_1^2+\theta_2^2)\right)$

1. Compute the distributions of $$\theta_1 \mid \theta_2$$ and $$\theta_2 \mid \theta_1$$

2. Write an R code based on the Gibbs algorithm to simulate a sample with a Standard bivariate Gaussian distribution.