A Worked Bivariate Gaussian Mixture Fit

From component probabilities to a fitted mixture

A bivariate Gaussian mixture describes observations in two dimensions as coming from one of several latent components. The fitting calculation has a simple pattern: use the current parameters to calculate how much each component accounts for each observation, then use those fractional contributions to update the parameters. The fractional contributions are called responsibilities. They are not observed labels, and an observation need not be assigned wholly to either component.

This tutorial works through one complete two-component update. It uses full covariance matrices, so the calculation tracks spread in both coordinates and their association. The goal is to make the arithmetic transparent: after each update, the resulting parameters can be checked against the data and the observed-data log-likelihood.

Definition: For observations $x_1,\ldots,x_n\in\mathbb{R}^2$, a two-component bivariate Gaussian mixture has density $$ f(x_i)=\pi_1\phi_2(x_i;\mu_1,\Sigma_1)+\pi_2\phi_2(x_i;\mu_2,\Sigma_2), \qquad \pi_1+\pi_2=1,\quad \pi_k>0, $$ where $\mu_k\in\mathbb{R}^2$, $\Sigma_k$ is positive definite, and $$ \phi_2(x;\mu,\Sigma)=\frac{1}{2\pi|\Sigma|^{1/2}} \exp\!\left\{-\frac12(x-\mu)^\mathsf{T}\Sigma^{-1}(x-\mu)\right\}. $$ The observed-data log-likelihood is $\ell=\sum_{i=1}^n\log f(x_i)$. Given current parameters, the responsibility of component $k$ for observation $i$ is $$ \tau_{ik}=\frac{\pi_k\phi_2(x_i;\mu_k,\Sigma_k)} {\sum_{j=1}^2\pi_j\phi_2(x_i;\mu_j,\Sigma_j)}. $$

For each observation, $\tau_{i1}+\tau_{i2}=1$. The responsibilities are calculated using the parameters before the update; the updated parameters are then calculated from those same responsibilities. This order is essential. In a computation with many observations or highly different component densities, log-sum-exp arithmetic is useful for numerical stability.

Worked Example: A Responsibility at an Overlap Point

Consider equal mixing proportions, identity covariance matrices, and component means $\mu_1=(0,0)^\mathsf{T}$ and $\mu_2=(2,0)^\mathsf{T}$. For the observation $x=(0.5,0)^\mathsf{T}$, calculate its responsibility for component 1.

The squared distances from $x$ to the means are $0.5^2=0.25$ and $(0.5-2)^2=2.25$. Since the proportions and covariance matrices agree, the common Gaussian normalizing factors cancel. The ratio of the component-1 density to the component-2 density is $$ \frac{\phi_2(x;\mu_1,I)}{\phi_2(x;\mu_2,I)} =\exp\!\left\{-\frac{0.25}{2}+\frac{2.25}{2}\right\} =e. $$ Thus $$ \tau_{1,1}=\frac{e}{e+1}\approx0.7311,\qquad \tau_{1,2}=\frac{1}{e+1}\approx0.2689. $$ The observation contributes to both components, with the larger contribution going to component 1. Replacing these values by a hard assignment would discard information used by the mixture update.

Why the weighted updates have this form

Suppose the current responsibilities are fixed, and write $N_k=\sum_{i=1}^n\tau_{ik}$. The part of the expected complete-data log-likelihood involving the parameters of component $k$ is, up to terms not involving those parameters, $$ Q_k=-\frac{N_k}{2}\log|\Sigma_k| -\frac12\sum_{i=1}^n\tau_{ik}(x_i-\mu_k)^\mathsf{T}\Sigma_k^{-1}(x_i-\mu_k). $$ The update maximizes this expression while respecting the mixture-proportion constraint.

Result: If $N_k>0$, the mixing-proportion and mean updates are $$ \pi_k^{\mathrm{new}}=\frac{N_k}{n}, \qquad \mu_k^{\mathrm{new}}=\frac{\sum_{i=1}^n\tau_{ik}x_i}{N_k}. $$ If the weighted scatter matrix $$ S_k=\frac{1}{N_k}\sum_{i=1}^n\tau_{ik} (x_i-\mu_k^{\mathrm{new}})(x_i-\mu_k^{\mathrm{new}})^\mathsf{T} $$ is positive definite, the unrestricted full-covariance update is $\Sigma_k^{\mathrm{new}}=S_k$.

Here is a derivation of these maximizers. For the proportions, maximize $\sum_kN_k\log\pi_k$ subject to $\sum_k\pi_k=1$. The Lagrange equations give $N_k/\pi_k=\lambda$; summing over $k$ yields $\lambda=n$, since $\sum_kN_k=n$. Hence $\pi_k=N_k/n$. For the mean, expanding the weighted quadratic around the weighted mean gives $$ \sum_i\tau_{ik}(x_i-\mu)(x_i-\mu)^\mathsf{T} = \sum_i\tau_{ik}(x_i-\mu_k^{\mathrm{new}})(x_i-\mu_k^{\mathrm{new}})^\mathsf{T} +N_k(\mu-\mu_k^{\mathrm{new}})(\mu-\mu_k^{\mathrm{new}})^\mathsf{T}. $$ The second term is positive semidefinite, so the quadratic contribution is minimized, and $Q_k$ maximized, at the weighted mean.

For the covariance, substitute that mean and let $S_k$ be the weighted scatter above. The covariance-dependent quantity to minimize is $\log|\Sigma|+\operatorname{tr}(\Sigma^{-1}S_k)$. For positive-definite $S_k$ and $\Sigma$, let $\lambda_1,\lambda_2$ be the eigenvalues of $\Sigma^{-1/2}S_k\Sigma^{-1/2}$. The difference from its value at $\Sigma=S_k$ is $$ \log|\Sigma|+\operatorname{tr}(\Sigma^{-1}S_k)-\log|S_k|-2 =\sum_{j=1}^2(\lambda_j-\log\lambda_j-1)\geq0, $$ because $u-\log u-1\geq0$ for $u>0$, with equality only at $u=1$. Thus $S_k$ is the maximizing covariance update. If the weighted scatter is singular, this unrestricted update is not positive definite; the result does not license using a singular covariance as a valid Gaussian density. Constraints or regularization require their own specified update.

Worked Example: A Weighted Mean and Covariance

Take three bivariate observations $x_1=(0,0)^\mathsf{T}$, $x_2=(2,0)^\mathsf{T}$, and $x_3=(0,2)^\mathsf{T}$. Suppose the responsibilities for one component are $0.8,0.6,0.2$, respectively. Then $N_k=1.6$, and its updated mean is $$ \mu_k^{\mathrm{new}} =\frac{0.8(0,0)^\mathsf{T}+0.6(2,0)^\mathsf{T}+0.2(0,2)^\mathsf{T}}{1.6} =(0.75,0.25)^\mathsf{T}. $$

The weighted second-moment matrix is $$ \frac{1}{1.6}\sum_i\tau_{ik}x_ix_i^\mathsf{T} =\frac{1}{1.6} \begin{pmatrix}2.4&0\\0&0.8\end{pmatrix} = \begin{pmatrix}1.5&0\\0&0.5\end{pmatrix}. $$ Subtracting $\mu_k^{\mathrm{new}}(\mu_k^{\mathrm{new}})^\mathsf{T}$ gives $$ \Sigma_k^{\mathrm{new}} = \begin{pmatrix}1.5&0\\0&0.5\end{pmatrix} - \begin{pmatrix}0.5625&0.1875\\0.1875&0.0625\end{pmatrix} = \begin{pmatrix}0.9375&-0.1875\\-0.1875&0.4375\end{pmatrix}. $$ Its determinant is $0.9375(0.4375)-0.1875^2=0.375>0$, and its first diagonal entry is positive, so this symmetric matrix is positive definite. The update is a valid full covariance matrix; in particular, its nonzero off-diagonal entries should not be silently discarded.

A complete two-component update

Now fit one EM iteration to four observations: $$ (0,0)^\mathsf{T},\quad(0,1)^\mathsf{T},\quad(3,0)^\mathsf{T},\quad(3,1)^\mathsf{T}. $$ Start with $\pi_1=\pi_2=0.5$, $\mu_1=(0,0)^\mathsf{T}$, $\mu_2=(3,0)^\mathsf{T}$, and $\Sigma_1=\Sigma_2=I$. The two component densities have equal normalizing constants. For either observation with first coordinate 0, the ratio of component-1 density to component-2 density is $\exp(4.5)$, since the squared first-coordinate distances are 0 and 9. For either observation with first coordinate 3, the ratio is $\exp(-4.5)$. The second coordinate affects both component densities equally and cancels in the responsibility ratio.

Worked Example: Responsibilities, Means, and Covariances

Let $q=\{1+\exp(-4.5)\}^{-1}\approx0.989013$, so $s=1-q\approx0.010987$. The responsibilities for component 1 are $q$ at each observation with first coordinate 0 and $s$ at each observation with first coordinate 3. Component 2 receives the complementary responsibilities. Each component has total responsibility $2q+2s=2$, so both updated proportions are $2/4=0.5$.

For component 1, the weighted mean has first coordinate $$ \frac{2s(3)}{2}=3s\approx0.032961 $$ and second coordinate $1/2$, because the total responsibility on the two observations with second coordinate 1 is $q+s=1$. Thus $\mu_1^{\mathrm{new}}\approx(0.032961,0.5)^\mathsf{T}$. By symmetry, $\mu_2^{\mathrm{new}}\approx(2.967039,0.5)^\mathsf{T}$.

For component 1, the weighted second moment of the first coordinate is $9s\approx0.098882$. Subtracting the squared mean gives the first-coordinate variance $$ 9s-(3s)^2\approx0.098882-0.001086=0.097796. $$ The second-coordinate variance is $1/2-(1/2)^2=0.25$. The weighted cross-moment is $3s/2$, which equals the product of the updated means, $(3s)(1/2)$; therefore the covariance is zero. Symmetry gives the same updated covariance for component 2: $$ \Sigma_1^{\mathrm{new}}=\Sigma_2^{\mathrm{new}} \approx \begin{pmatrix}0.097796&0\\0&0.25\end{pmatrix}. $$ Both diagonal entries are positive, so these updated covariance matrices are positive definite.

Checking the observed-data likelihood

The likelihood calculation uses the mixture density at each observation, not the complete-data likelihood with a component label chosen for each point. For the starting parameters, at either point with second coordinate 0 the density is $$ f_{\mathrm{old}}(x)=\frac{1}{4\pi}\{1+\exp(-4.5)\}, $$ and at either point with second coordinate 1 it is smaller by a factor $\exp(-0.5)$. Hence $$ \ell_{\mathrm{old}} =2\log\!\left[\frac{1+\exp(-4.5)}{4\pi}\right] +2\log\!\left[\frac{\exp(-0.5)\{1+\exp(-4.5)\}}{4\pi}\right] \approx-11.0799. $$

After the update, write $v=0.097796$, $m=3s$, and $C=1/(\pi\sqrt{v})$. At a point with first coordinate 0 or 3 and second coordinate 0, the mixture density is $d_0=\frac{C}{2}\left[\exp\{-m^2/(2v)\}+\exp\{-(3-m)^2/(2v)\}\right]$. At the corresponding point with second coordinate 1, the density is smaller by $\exp(-0.5)$, because the updated second-coordinate variance is $0.25$ and the displacement from its mean is $0.5$. Therefore $$ \ell_{\mathrm{new}} =4\log(d_0)-1 \approx-3.7240. $$ The increase is approximately $7.3559$. This is a substantial change for this small example; it is not evidence that the update is a global maximum or that every fitted mixture will have well-separated components.

Why the likelihood does not decrease

The EM monotonicity statement concerns an update that maximizes the expected complete-data log-likelihood, or at least does not decrease it, over the allowed parameter space. To see why it applies here, let $z$ denote a component label and let $p_\theta(x_i,z=k)=\pi_k\phi_2(x_i;\mu_k,\Sigma_k)$. Assume the current and updated mixture densities are positive at the observed data, the current responsibilities are well-defined, and the proposed parameters are admissible.

Result: For an EM update from $\theta$ to $\theta'$, $$ \ell(\theta')-\ell(\theta) \geq Q(\theta'\mid\theta)-Q(\theta\mid\theta), $$ where $Q(\theta'\mid\theta)=\sum_{i,k}\tau_{ik}(\theta)\log p_{\theta'}(x_i,z=k)$. Thus an update that does not decrease $Q$ cannot decrease the observed-data log-likelihood.

For each observation, use the current posterior probabilities $\tau_{ik}(\theta)=p_\theta(x_i,z=k)/p_\theta(x_i)$. Then $$ \frac{p_{\theta'}(x_i)}{p_\theta(x_i)} =\sum_k\tau_{ik}(\theta) \frac{p_{\theta'}(x_i,z=k)}{p_\theta(x_i,z=k)}. $$ Apply the concavity of the logarithm to this finite weighted sum: $$ \log\frac{p_{\theta'}(x_i)}{p_\theta(x_i)} \geq \sum_k\tau_{ik}(\theta) \log\frac{p_{\theta'}(x_i,z=k)}{p_\theta(x_i,z=k)}. $$ Summing over $i$ gives the stated inequality, because the right-hand side is $Q(\theta'\mid\theta)-Q(\theta\mid\theta)$. The weighted updates above maximize $Q$ when they are admissible, so their observed-data log-likelihood cannot be lower than the starting value. This monotonicity argument does not assert that each iteration finds the global likelihood maximum.

Common mistakes and a numerical check

Using hard assignments. Choosing the larger responsibility and setting the other to zero changes the weighted update. Use both posterior probabilities unless the model or method explicitly calls for hard labels.
Using the wrong denominator. Divide a weighted mean or covariance by $N_k=\sum_i\tau_{ik}$, not by the sample size $n$. The mixing proportion, in contrast, is $N_k/n$.
Confusing covariance with second moments. The weighted covariance is centered at the updated mean. A weighted second-moment matrix alone is not the covariance.
Rounding before aggregating. Retain unrounded responsibilities and log densities through the sum, then round the reported result. For example, if four equal per-observation log densities give an unrounded total of about $-10.554636$, report $-10.5546$, not a total formed by multiplying the prematurely rounded value $-2.6387$. Against an old log-likelihood of $-11.6163849$, the increase is about $1.0617485$, which rounds to $1.0617$.

A numerical update should be checked at several levels: responsibilities sum to one within each observation; component responsibility totals sum to $n$; updated mixing proportions sum to one; covariance matrices are symmetric and, when an ordinary Gaussian density is required, positive definite; and the observed-data likelihood is evaluated as a sum of log mixture densities.

Check Your Understanding

Use the definitions and calculations above to answer these questions.

Why does a responsibility use the mixing proportion as well as the component density?
If a component has total responsibility $N_k=1.8$ in a sample of six observations, what is its updated mixing proportion?
What condition on the weighted scatter matrix is needed for the unrestricted full-covariance update to be positive definite?
In the four-point example, why does the second-coordinate variance update to $0.25$ for both components?
What does the EM monotonicity inequality guarantee, and what does it not guarantee about the likelihood maximum?