STAT 238 - Bayesian Statistics Lecture Twenty Seven

Spring 2026, UC Berkeley

Gibbs Sampler for 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 posterior of $(w, \beta)$ given the data ($\text{data} = y = (y_1, \dots, y_n)$) is

where $X$ is, as usual, the design matrix with rows $x_1', \dots, x_n'$.

We shall now run the Gibbs sampler with the decomposition $\theta = (w, \beta)$. For this, we need to sample from the conditionals $\beta|w, \text{data}$ and $w | \beta, \text{data}$.

For a fixed $w \in \R^n$, let

so that

Thus

In other words,

The next task is to figure out the conditional $w | \beta, \data$:

Thus $w_1, \dots, w_n$ are independent conditional on $\beta$ and $\data$. Further

It makes sense now to separate the two cases $y_i = 1$ and $y_i = 0$. For $y_i = 1$, we get

Thus, conditional on $\beta$ and $y_i = 1$, the random variable $w_i$ has the $N(x_i' \beta, 1)$ distribution trucated to the positive half-line.

Here is a fact (from the wikipedia article for truncated normal distributions: Truncated normal distribution#Generating values from the truncated normal distribution) that is useful for simulating from truncated normal distributions.

Using this fact, it is straightforward to verify that $w_i|\{\beta, y_i = 1\}$ can be simulated by

This follows directly from the fact above by taking $\mu = x_i' \beta$, $\sigma = 1$, $a = 0$ and $b = \infty$ (this implies that $a_0 = -x_i' \beta$ and $b_0 = \infty$).

One can argue analogously to deduce that $w_i | \{\beta, y_i = 0\}$ can be simulated by

The Gibbs sampler for the probit model can be therefore implemented by the following model:

1. Initialize at $\beta^{(0)}$ and $w^{(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((X'X)^{-1} X' w^{(t)}, (X' X)^{-1})$ and $w^{(t+1)} = w^{(t)}$.

   3. If $J = 2$, take $\beta^{(t+1)} = \beta^{(t)}$ and

      $$w_i = \begin{cases} x_i' \beta^{(t)} + \Phi^{-1} \left((1 - U_i) \Phi(-x_i' \beta^{(t)}) + U_i \right) & \text{if } y_i = 1 \\ x_i' \beta^{(t)} + \Phi^{-1} \left(U_i \Phi(-x_i' \beta^{(t)}) \right) & \text{if } y_i = 0 \end{cases}$$

      (4)

      where $U_i \sim Unif(0, 1)$.