R
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?
Choose arbitrary starting values: \(\theta^{(0)}=(\theta_1^{(0)},\ldots,\theta_k^{(0)})\)
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)}\).
To better understand the Gibbs Sampler, I suggest to read the paper by Casela and George 1992
\(\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
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
\[\mu\mid y_1,\ldots,y_n,\tau \sim \mathcal{N}\left(\overline{y},\displaystyle\frac{1}{n\tau}\right)\]
\[\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)
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)\]
Compute the distributions of \(\theta_1 \mid \theta_2\) and \(\theta_2 \mid \theta_1\)
Write an R
code based on the Gibbs algorithm to simulate a sample with a Standard bivariate Gaussian distribution.