STAT 238 - Bayesian Statistics Lecture Thirty

Spring 2026, UC Berkeley

The Gibbs Sampler for Gaussian Mixture Models¶

We observe real-valued data $y_1, \dots, y_n$ and we consider the model:

with unknown parameters $(w_1, \dots, w_k)$ and $(\mu_j, \sigma_j^2)$ for $j = 1, \dots, k$. $(w_1, \dots, w_k)$ is a probablity vector (i.e., $w_j \geq 0$ and $\sum_j w_j = 1$).

We use the following priors:

The default choices are $a_j = 1$ for all $j$, $m = 0$ and $s^2$ to be large, and $\alpha, \beta$ to be near zero. Note that the density of the Inverse Gamma distribution $IG(\alpha, \beta)$ is given by $x^{-\alpha - 1} \exp(-\beta/x) I\{x > 0\}$. It is easy to check that this is simply the density of $X^{-1}$ where $X \sim \text{Gamma}(\alpha, \beta)$.

To use the Gibbs sampler, we introduce latent variables $z_1, \dots, z_n$ which represent group memberships. Specifically, $z_i$ represents which of the $k$ Gaussians the observation $y_i$ is coming from. More precisely, $z_i$ takes the values $1, \dots, k$ with probabilities $w_1, \dots, w_k$ (i.e., $z_i \sim \text{Categorical}(w)$), and

The Gibbs sampler will then jointly sample from all the unknown variables:

The full conditionals corresponding to all the variables can be written in closed form as described below. We will use the notation $w = (w_1, \dots, w_k)$, $\mu = (\mu_1, \dots, \mu_k)$ and $\sigma = (\sigma_1, \dots, \sigma_k)$, also $y = (y_1, \dots, y_n)$ and $z = (z_1, \dots, z_n)$.

For the full conditional corresponding to $z_1, \dots, z_n$ (i.e., the conditional distribution of $z_1, \dots, z_n$ given $y_1, \dots, y_n$ as well as $w_1, \dots, w_k, \mu_1, \sigma_1^2, \dots, \mu_k, \sigma_k^2$) is given by:

In other words,

The full conditional of $w = (w_1, \dots, w_k)$ (i.e., the conditional distribution of $w$ given $z, y, \mu, \sigma$) is:

The full conditional of $\mu_j$ (i.e., the conditional distribution of $\mu_j$ given $z, y, \sigma$) is:

Note that we are also conditioning on $\sigma_j$ above. If we do not condition on $\sigma_j$, then the conditional distribution of $\mu_j$ will be $t$-distributed and not normal distributed.

The full conditional of $\sigma_j^2$ (i.e., the conditional distribution of $\sigma_j^2$ given $z, y, \mu$) is:

where again $n_j = \sum_{i=1}^n I\{z_i = j\}$.

So the Gibbs sampler algorithm is given by:

1. Initialize the parameters $w^{(0)}, \mu^{(0)}, \sigma^{(0)}$.

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

   1. Generate $z_1^{(t)}, \dots, z_n^{(t)}$ via:

      $$\begin{align*} z_i^{(t)} \sim \text{Categorical} \left(\frac{w_j^{(t)} \phi \left(y_i, \mu^{(t)}_j, (\sigma_j^{(t)})^2 \right)}{\sum_{j=1}^k w_j^{(t)} \phi \left(y_i, \mu^{(t)}_j, (\sigma_j^{(t)})^2 \right)}, j = 1, \dots, k \right) \end{align*}$$

      (10)

   2. Calculate

      $$\begin{align*} n_j^{(t)} := \sum_{i=1}^n I\{z_i^{(t)} = j\} ~~ \text{ for $j = 1, \dots, k$}. \end{align*}$$

      (11)

   3. Generate $w^{(t+1)} = (w_1^{(t+1)}, \dots, w_k^{(t+1)})$ via

      $$\begin{align*} w^{(t+1)} \sim \text{Dirichlet}(a_1 + n_1^{(t)}, \dots, a_k + n_k^{(t)}). \end{align*}$$

      (12)

   4. Generate $\mu_1^{(t+1)}, \dots, \mu_k^{(t+1)}$ via

      $$\begin{align*} \mu_j^{(t+1)} \sim N \left(\frac{m/s^2 + \sum_{i : z^{(t)}_i = j} y_i/(\sigma_j^{(t)})^2}{1/s^2 + n_j^{(t)}/(\sigma_j^{(t)})^2}, \frac{1}{1/s^2 + n^{(t)}_j/(\sigma_j^{(t)})^2} \right). \end{align*}$$

      (13)

   5. Generate $\sigma_1^{(t+1)}, \dots, \sigma_k^{(t+1)}$ via

      $$\begin{align*} (\sigma_j^2)^{(t+1)} \sim IG \left(\alpha + \frac{n^{(t)}_j}{2}, \beta + \frac{1}{2} \sum_{i : z^{(t)}_i = j} (y_i - \mu^{(t+1)}_j)^2 \right), \end{align*}$$

      (14)

The EM Algorithm¶

The EM algorithm can be seen as a deterministic variant of the Gibbs sampler. We shall look at the EM in more generality later. In the case of fitting finite normal mixtures, the EM algorithm is given by the following:

1. Initialize the parameters $w^{(0)}, \mu^{(0)}, \sigma^{(0)}$.

2. Repeat the following for $t = 0, 1, 2, \dots$ until convergence:

   1. E-step. Compute the responsibility of component $j$ for observation $i$:

      $$\begin{align*} r_{ij}^{(t)} := \frac{w_j^{(t)} \phi\left(y_i, \mu_j^{(t)}, (\sigma_j^{(t)})^2\right)} {\sum_{l=1}^k w_l^{(t)} \phi\left(y_i, \mu_l^{(t)}, (\sigma_l^{(t)})^2\right)}, \quad j = 1, \dots, k. \end{align*}$$

      (15)

   2. Compute the effective counts

      $$\begin{align*} n_j^{(t)} := \sum_{i=1}^n r_{ij}^{(t)}, \quad j = 1, \dots, k. \end{align*}$$

      (16)

   3. M-step. Update the weights:

      $$\begin{align*} w_j^{(t+1)} := \frac{n_j^{(t)}}{n}, \quad j = 1, \dots, k. \end{align*}$$

      (17)

   4. Update the means:

      $$\begin{align*} \mu_j^{(t+1)} := \frac{\sum_{i=1}^n r_{ij}^{(t)}\, y_i}{n_j^{(t)}}, \quad j = 1, \dots, k. \end{align*}$$

      (18)

   5. Update the variances:

      $$\begin{align*} (\sigma_j^2)^{(t+1)} := \frac{\sum_{i=1}^n r_{ij}^{(t)} \left(y_i - \mu_j^{(t+1)}\right)^2}{n_j^{(t)}}, \quad j = 1, \dots, k. \end{align*}$$

      (19)

There are the following close connections between the EM algorithm and the Gibbs sampler:

• The E-step replaces sampling $z_i \sim \text{Categorical}(\ldots)$ with computing its expectation

  $$\begin{align*} r_{ij}^{(t)} = \mathbb{E}\left[I\{z_i = j\} \mid y_i, w^{(t)}, \mu^{(t)}, \sigma^{(t)}\right], \end{align*}$$

  (20)

  which are the same normalised weights used in the Categorical draw.

• The effective count $n_j^{(t)} = \sum_{i=1}^n r_{ij}^{(t)}$ replaces the integer count $\sum_{i=1}^n I\{z_i^{(t)} = j\}$.

• The weight update replaces the Dirichlet draw and takes the maximum likelihood estimate

  $$\begin{align*} w_j^{(t+1)} = \frac{n_j^{(t)}}{n}, \end{align*}$$

  (21)

  which is equivalent to the posterior mean of the Dirichlet as the prior becomes uninformative.

• The mean update replaces the Gaussian conjugate draw and is the responsibility-weighted sample mean

  $$\begin{align*} \mu_j^{(t+1)} = \frac{\sum_{i=1}^n r_{ij}^{(t)}\, y_i}{n_j^{(t)}}, \end{align*}$$

  (22)

  the maximum likelihood estimate with no prior shrinkage toward $m$.

• The variance update replaces the inverse-Gamma draw and is the responsibility-weighted sample variance

  $$\begin{align*} (\sigma_j^2)^{(t+1)} = \frac{\sum_{i=1}^n r_{ij}^{(t)} \left(y_i - \mu_j^{(t+1)}\right)^2}{n_j^{(t)}}, \end{align*}$$

  (23)

  the maximum likelihood estimate with no $\alpha, \beta$ regularisation terms.