STAT 238 - Bayesian Statistics Lecture Twenty Six

Spring 2026, UC Berkeley

The Gibbs Sampler¶

Let us introduce some notation. The goal is to obtain samples from a target distribution with density $\pi$ (in our applications, $\pi$ will be the posterior density). $\pi$ represents the density on $\R^d$ of a random vector $\Theta$ of size $d \times 1$:

We assume that $\Theta$ can be decomposed as:

for $k$ subvectors $\Theta_{(1)}, \dots, \Theta_{(k)}$. The $j^{th}$ subvector is $\Theta_{(j)}$. We shall also denote the vector obtained by dropping the $j^{th}$ subvector $\Theta_{(j)}$ from $\Theta$ by $\Theta_{-(j)}$. For example,

Gibbs sampling works according to the following algorithm.

1. Initialize with an arbitrary value $\theta^{(0)}$.

2. Repeat the following for $t = 1, \dots, N$ (for a large number $N$):

   1. Pick $J$ uniformly at random from $1, \dots, k$. Suppose $J$ turned out to be $j$.

   2. Set $\theta^{(t+1)}_{-(j)} = \theta^{(t)}_{-(j)}$ i.e., $\theta^{(t+1)}_{(i)} = \theta^{(t)}_{(i)}$ for $i \neq j$.

   3. $\theta^{(t+1)}_{(j)}$ is randomly drawn according to the conditional distribution of $\Theta_{(j)}$ given $\Theta_{-(j)} = \theta^{(t)}_{-(j)}$.

This algorithm assumes that we have some way of generating observations from the conditional distribution of $\Theta_{(j)}$ given $\Theta_{-(j)}$ for each $j$.

Examples of the Gibbs Sampler¶

In the first two examples below, the posterior can be written in closed form so MCMC is not really necessary. These are still useful for evaluating the performance of the Gibbs sampler. The third example is more non-trivial and illustrates the usefulness of the Gibbs sampler.

Bivariate Normal¶

We want to sample $(\theta_1, \theta_2)$ from the following bivariate distribution:

Here $y = (y_1, y_2)^T$ is fixed. One can, of course, directly sample from this normal distribution without any MCMC. We can use this example to illustrate the performance of the Gibbs sampler. To use the Gibbs sampler, observe that

The Gibbs sampler algorithm then would be the following:

1. Initialize at arbitrary $\theta_1^{(0)}$ and $\theta_2^{(0)}$.

2. Repeat the following for $t = 1, \dots, N$:

   1. Pick $J$ to be either 1 or 2 with probability 0.5.

   2. If $J = 1$, take $\theta_2^{(t+1)} = \theta_2^{(t)}$ and generate $\theta_1^{(t+1)}$ from the normal distribution with mean $y_1 + \rho (\theta_2^{(t)} - y_2)$ and variance $1 - \rho^2$.

   3. If $J = 2$, take $\theta_1^{(t+1)} = \theta_1^{(t)}$ and generate $\theta_2^{(t+1)}$ from the normal distribution with mean $y_2 + \rho (\theta_1^{(t)} - y_1)$ and variance $1 - \rho^2$.

This example is also useful for illustrating one potential problem with the Gibbs sampler. Suppose $\rho = 1$. In this case $\theta_1 - y_1 \sim N(0, 1)$ and $\theta_2 - y_2 = \theta_1 - y_1$. The Gibbs iterates in this case will be

and

This results in the chain not moving after the first iterate. For example, if the first step is:

then the chain will stay at $\begin{pmatrix} y_1 - y_2 + \theta_2^{(0)}\\ \theta_2^{(0)}\end{pmatrix}$ forever. On the other hand, if the first step is

then the chain will stay forever at $\begin{pmatrix} \theta_1^{(0)}\\ y_2 - y_1 + \theta_1^{(0)}\end{pmatrix}$. This chain will obviously not converge unless the initial iterate already has the right distribution. A similar problem arises when $\rho = -1$.

The lesson to be drawn from this is that the Gibbs sampler has issues of convergence when there is high correlation between the variables.

Linear Regression¶

Consider the linear regression model that we discussed previously. We observe data $(y_1, x_1), \dots, (y_n, x_n)$ with $y_1, \dots, y_n$ real-valued and $x_1, \dots, x_n \in \R^p$. We treat $x_1, \dots, x_n$ as non-random and $y_1, \dots, y_n$ as random. The parameter is $\theta \in \R^{p+1}$ with $\theta = (\beta, \sigma)$ (here $\beta \in \R^p$ represents the vector of regression coefficients and $\sigma$ is the noise standard deviation). Consider the model

with $y_1, \dots, y_n$ being independent conditional on $\theta$ along with the (improper) prior:

The exact joint posterior is given by

We can directly sample from this posterior (e.g., by first sampling from the marginal posterior of $\sigma$ and then from the conditional posterior of $\beta$ given $\sigma$. Alternatively, one can first sample from the marginal posterior of $\beta$ and then from the conditional posterior of $\sigma$ given $\beta$).

Now let us use the Gibbs sample to simulate from (1). Note that:

which gives

Here $\hat{\beta}$ is the least squares estimator. For the conditional posterior of $\sigma$ given $\beta$, note that

Using the formula $f_{1/\sigma^2}(x) = f_{\sigma}(x^{-1/2}) x^{-3/2}$, it is now easy to check that the density of $1/\sigma^2$ conditional on $\beta$ and the data is proportional to

which means that

The Gibbs sampler algorithm is thus:

1. Initialize at arbitrary $\beta^{(0)}$ and $\sigma^{(0)}$.

2. Repeat the following for $t = 1, \dots, N$:

   1. Pick $J$ to be either 1 or 2 with equal probability.

   2. If $J =1$, take $\beta^{(t+1)} \sim N_p(\hat{\beta}, (\sigma^{(t)})^2 (X' X)^{-1})$ and $\sigma^{(t+1)} = \sigma^{(t)}$.

   3. If $J = 2$, take $\beta^{(t+1)} = \beta^{(t)}$. Generate an observation $x$ from the Gamma distribution with parameters $n/2$ and $\|Y - X\beta^{(t)}\|^2/2$. Then take $\sigma^{(t+1)} = x^{-1/2}$.

Probit Regression¶

Our next application of the Gibbs sampler is to Probit Regression which was introduced in Homework Four. The model is given by the following. We observe data $(y_1, x_1), \dots, (y_n, x_n)$ with $y_1, \dots, y_n$ binary and $x_1, \dots, x_n \in \R^p$. $x_1, \dots, x_n$ are treated non-random and $y_1, \dots, y_n$ random. The parameter is $\beta \in \R^p$. The probit regression model assumes that $y_1 \dots, y_n$ are independent given $\beta$ with

where $\Phi(\cdot)$ is the standard normal cdf. Let $\pi(\cdot)$ be our prior for $\beta$. The posterior is then given by:

Suppose we take the (improper) uniform prior for $\beta$ so that the posterior becomes

It is not actually clear how to use the Gibbs sampler on the above posterior. We can try decomposing $\beta$ as $(\beta_1, \dots, \beta_p)$ and then attempt to obtain samples from the conditional:

I am not sure how to obtain samples from the above conditional.

There is a trick that people use to run the Gibbs sampler for probit regression. Let us introduce random variables $w_1, \dots, w_n$ which, conditional on $\beta$, are independent with

We then let

It is then clear that $y_1, \dots, y_n$ defined thus have the same likelihood as in the probit regression model. Let $w = (w_1, \dots, w_n)$. The trick is to use the Gibbs sampler with $(w, \beta)$. We shall see how this is done in the next class.