import plotly.io as pio; pio.renderers.default = "notebook_connected"
Lecture 05: Density Estimation and Gaussian Mixture Models – CS 189, Fall 2026
Accompanying demonstration: fitting continuous data, one Gaussian and then several.
Lecture 4 estimated the parameters of a Bernoulli, a distribution over a single binary variable. This notebook applies the same machinery to continuous data. It fits a Gaussian by maximum likelihood, shows that the resulting variance is biased, then meets data that no single Gaussian describes and fits a mixture with the EM algorithm.
Everything is one-dimensional, matching the slides. The EM implementation in section 6 is written for general $D$, so it also works on multivariate data, but every example here uses $D = 1$.
import numpy as np
import pandas as pd
import plotly.express as px
import plotly.graph_objects as go
from scipy.stats import multivariate_normal
rng = np.random.default_rng(0)
BLUE, GOLD, GREEN, ORANGE = '#002675', '#FDB515', '#028842', '#DF5327'
# A separate stream for building the dataset, so the figures on the slides and the
# numbers here stay identical no matter what else is drawn along the way.
data_rng = np.random.default_rng(189)
1. Fitting a Gaussian by maximum likelihood¶
The feature is the loudness of a one-second audio segment in decibels. The data is simulated, since a recording of the lecture hall cannot be distributed with the notebook, but the shape is representative of a quiet room.
loudness = data_rng.normal(-52, 4.0, size=4000) # background segments, in dB
px.histogram(loudness, nbins=60, histnorm='probability density',
color_discrete_sequence=[BLUE], labels={'value': 'loudness (dB)'},
title='Loudness of 4,000 background segments',
width=850, height=430, template='plotly_white')
The distribution is unimodal and roughly symmetric, which suggests a Gaussian. The maximum likelihood estimates derived in lecture are the sample mean and the sample variance.
mu_ml = loudness.mean()
var_ml = ((loudness - mu_ml) ** 2).mean() # note: divide by N, not N - 1
print(f'mu_ML = {mu_ml:.3f}')
print(f'sigma2_ML = {var_ml:.3f} (sigma = {np.sqrt(var_ml):.3f})')
mu_ML = -51.928 sigma2_ML = 16.090 (sigma = 4.011)
Two numbers, computed in a single pass. They are the sufficient statistics: given $\sum_n x_n$ and $\sum_n x_n^2$, the 4,000 observations are not needed to write down the fitted density.
def gaussian(x, mu, var):
return np.exp(-0.5 * (x - mu) ** 2 / var) / np.sqrt(2 * np.pi * var)
grid = np.linspace(loudness.min() - 3, loudness.max() + 3, 400)
fig = px.histogram(loudness, nbins=60, histnorm='probability density',
color_discrete_sequence=['#c9d3e6'], labels={'value': 'loudness (dB)'})
fig.add_trace(go.Scatter(x=grid, y=gaussian(grid, mu_ml, var_ml),
line=dict(color=BLUE, width=3), name='fitted Gaussian'))
fig.update_layout(title='Maximum likelihood fit', width=850, height=430,
template='plotly_white', showlegend=False)
fig
The optimization surface¶
The closed form conceals a two-dimensional optimization over $\mu$ and $\sigma$. Evaluating the log likelihood on a grid shows the surface whose maximum it locates.
mus = np.linspace(mu_ml - 1.0, mu_ml + 1.0, 120)
sigmas = np.linspace(np.sqrt(var_ml) - 0.8, np.sqrt(var_ml) + 0.8, 120)
MU, SIG = np.meshgrid(mus, sigmas)
n = len(loudness)
# Sum_n (x_n - mu)^2 = Sum x_n^2 - 2 mu Sum x_n + n mu^2, so the two sufficient
# statistics are all that is required; the observations never enter the grid.
s1, s2 = loudness.sum(), (loudness ** 2).sum()
sq = s2 - 2 * MU * s1 + n * MU ** 2
LL = -0.5 * n * np.log(2 * np.pi * SIG ** 2) - sq / (2 * SIG ** 2)
fig = go.Figure(go.Contour(x=mus, y=sigmas, z=LL, ncontours=40,
colorscale='Blues', showscale=False))
fig.add_trace(go.Scatter(x=[mu_ml], y=[np.sqrt(var_ml)], mode='markers',
marker=dict(color=GOLD, size=14, line=dict(color='black', width=1))))
fig.update_xaxes(title='mu'); fig.update_yaxes(title='sigma')
fig.update_layout(title='Log likelihood surface, with the closed-form solution marked',
width=760, height=520, template='plotly_white', showlegend=False)
fig
2. The maximum likelihood variance is biased¶
$\sigma^2_{ML}$ divides by $N$. Lecture established that $\mathbb{E}[\sigma^2_{ML}] = \frac{N-1}{N}\sigma^2$, so on average it underestimates. The deviations are measured from $\mu_{ML}$, which was fitted to the same data and has moved toward those observations.
The claim concerns an expectation over datasets, so it is checked by drawing many datasets.
TRUE_MU, TRUE_VAR = -52.0, 16.0
def average_variance_estimates(N, reps=4000):
s = rng.normal(TRUE_MU, np.sqrt(TRUE_VAR), size=(reps, N))
return s.var(axis=1).mean(), s.var(axis=1, ddof=1).mean()
Ns = np.array([2, 3, 5, 10, 20, 50, 100])
results = np.array([average_variance_estimates(N) for N in Ns])
fig = go.Figure()
fig.add_trace(go.Scatter(x=Ns, y=results[:, 0], mode='lines+markers',
name='sigma2_ML (divide by N)', line=dict(color=ORANGE)))
fig.add_trace(go.Scatter(x=Ns, y=results[:, 1], mode='lines+markers',
name='unbiased (divide by N-1)', line=dict(color=GREEN)))
fig.add_trace(go.Scatter(x=Ns, y=TRUE_VAR * (Ns - 1) / Ns, mode='lines',
name='theory: sigma2 (N-1)/N', line=dict(color=BLUE, dash='dot')))
fig.add_hline(y=TRUE_VAR, line_dash='dash', annotation_text='true sigma2 = 16')
fig.update_xaxes(type='log', title='N (samples per dataset)', tickmode='array',
tickvals=Ns, ticktext=[str(x) for x in Ns], minor_showgrid=False)
fig.update_yaxes(title='average estimate over 4,000 datasets')
fig.update_layout(title='The MLE variance is low by a factor of (N-1)/N',
width=850, height=480, template='plotly_white')
fig
The simulated averages follow the theoretical curve. The bias vanishes as $N$ grows and matters
only for small datasets. It is also the origin of the ddof argument: np.var with the default
ddof=0 is the maximum likelihood estimate, ddof=1 is the unbiased one, and Pandas .var()
defaults to ddof=1.
3. When one Gaussian is the wrong model¶
Everything above assumed the correct family. Here is a full day of loudness measurements, including the segments where somebody was speaking. This is the dataset on the transition slide.
speech = data_rng.normal(-30, 4.5, size=1800) # someone talking
day = np.concatenate([loudness, speech])
data_rng.shuffle(day)
mu_day, var_day = day.mean(), day.var()
print(f'N = {len(day)} mu_ML = {mu_day:.1f} sigma_ML = {np.sqrt(var_day):.1f}')
N = 5800 mu_ML = -45.1 sigma_ML = 11.0
grid2 = np.linspace(day.min() - 3, day.max() + 3, 600)
fig = px.histogram(day, nbins=90, histnorm='probability density',
color_discrete_sequence=['#c3ccdd'], labels={'value': 'loudness (dB)'})
fig.add_trace(go.Scatter(x=grid2, y=gaussian(grid2, mu_day, var_day),
line=dict(color=ORANGE, width=4), name='best single Gaussian'))
fig.add_vline(x=mu_day, line_dash='dot', line_color=ORANGE)
fig.update_layout(title='The best single Gaussian puts its peak where the data is sparsest',
width=880, height=470, template='plotly_white', showlegend=False)
fig
The fitted density places its mode between the two modes of the data, at a loudness that is rarely observed. This is not a failure of the optimization: it is the best available single Gaussian, and the family was the wrong choice.
One remedy abandons the parametric family and lets the data determine the shape.
from scipy.stats import gaussian_kde
fig = px.histogram(day, nbins=90, histnorm='probability density',
color_discrete_sequence=['#c3ccdd'], labels={'value': 'loudness (dB)'})
for bw, color, dash in [(0.05, GOLD, 'dot'), (0.25, GREEN, 'solid'), (1.0, ORANGE, 'dash')]:
fig.add_trace(go.Scatter(x=grid2, y=gaussian_kde(day, bw_method=bw)(grid2),
line=dict(color=color, dash=dash), name=f'bandwidth {bw}'))
fig.update_layout(title='Kernel density estimation at three bandwidths',
width=880, height=470, template='plotly_white')
fig
A small bandwidth tracks individual observations; a large one merges the two modes. The bandwidth plays the role that bin width plays in a histogram. The cost is that the estimate retains all 5,800 observations, where the Gaussian needed two numbers.
The other remedy keeps a parametric model but enlarges the family, which is the Gaussian mixture.
4. A mixture of two Gaussians¶
A mixture assigns each component a weight $\pi_k$, a mean $\mu_k$, and a variance $\sigma_k^2$:
$$p(x) = \sum_{k=1}^{K} \pi_k\, \mathcal{N}(x \mid \mu_k, \sigma_k^2)$$
Fitting one with scikit-learn takes a single call. The data must be shaped $(N, D)$ even when $D = 1$.
from sklearn.mixture import GaussianMixture
X = day.reshape(-1, 1) # (N, 1): one feature
gmm = GaussianMixture(n_components=2, random_state=0).fit(X)
for k in np.argsort(gmm.means_.ravel()):
print(f'component {k}: pi = {gmm.weights_[k]:.3f} '
f'mu = {gmm.means_[k, 0]:7.2f} sigma = {np.sqrt(gmm.covariances_[k, 0, 0]):.2f}')
component 1: pi = 0.690 mu = -51.92 sigma = 4.02 component 0: pi = 0.310 mu = -29.86 sigma = 4.56
Compare those against the values the data was generated from: weights 4000/5800 and 1800/5800, means -52 and -30, standard deviations 4.0 and 4.5.
mix = np.exp(gmm.score_samples(grid2.reshape(-1, 1)))
fig = px.histogram(day, nbins=90, histnorm='probability density',
color_discrete_sequence=['#c3ccdd'], labels={'value': 'loudness (dB)'})
fig.add_trace(go.Scatter(x=grid2, y=gaussian(grid2, mu_day, var_day),
line=dict(color=ORANGE, width=3, dash='dash'), name='single Gaussian'))
fig.add_trace(go.Scatter(x=grid2, y=mix, line=dict(color=BLUE, width=4), name='mixture of 2'))
for k, c in zip(range(2), [GREEN, GOLD]):
comp = gmm.weights_[k] * gaussian(grid2, gmm.means_[k, 0], gmm.covariances_[k, 0, 0])
fig.add_trace(go.Scatter(x=grid2, y=comp, line=dict(color=c, dash='dot'),
name=f'component {k}'))
fig.update_layout(title='One Gaussian against a mixture of two',
width=880, height=470, template='plotly_white')
fig
Soft assignments¶
predict_proba gives the responsibility of each component for each point. In one dimension we can
plot it against $x$ directly, which shows exactly where the model is uncertain.
resp = gmm.predict_proba(grid2.reshape(-1, 1))
order = np.argsort(gmm.means_.ravel())
fig = go.Figure()
for j, k in enumerate(order):
fig.add_trace(go.Scatter(x=grid2, y=resp[:, k], line=dict(color=[GREEN, GOLD][j], width=3),
name=f'p(z = {j} | x)'))
fig.add_hline(y=0.5, line_dash='dot', line_color='#888')
fig.update_xaxes(title='loudness (dB)')
fig.update_yaxes(title='responsibility', range=[-0.03, 1.03])
fig.update_layout(title='Responsibilities: the crossover is where a hard assignment would be a guess',
width=880, height=440, template='plotly_white')
fig
The curves are near 0 or 1 almost everywhere and swap over a narrow band around -40 dB. Points in that band are genuinely ambiguous, and a hard assignment would discard exactly that information. This is the difference between k-means and a mixture model.
5. A GMM is a generative model¶
The mixture factorizes as $p(x, z) = p(z)\,p(x \mid z)$, so it can be sampled by ancestor sampling: draw the component first, then draw the point from that component.
pi_true = np.array([0.2, 0.5, 0.3])
mu_true = np.array([-1.0, 2.0, 5.0])
var_true = np.array([0.2, 0.5, 0.1])
N = 3000
z = rng.choice(3, size=N, p=pi_true) # first the latent component
x = rng.normal(mu_true[z], np.sqrt(var_true[z])) # then the observation
pd.Series(z).value_counts(normalize=True).sort_index().round(3)
0 0.204 1 0.491 2 0.304 Name: proportion, dtype: float64
grid3 = np.linspace(x.min() - 1, x.max() + 1, 600)
density = sum(pi_true[k] * gaussian(grid3, mu_true[k], var_true[k]) for k in range(3))
fig = px.histogram(x, nbins=90, histnorm='probability density',
color_discrete_sequence=['#c9d3e6'], labels={'value': 'x'})
fig.add_trace(go.Scatter(x=grid3, y=density, line=dict(color=BLUE, width=3), name='mixture'))
for k, c in enumerate([GOLD, GREEN, ORANGE]):
fig.add_trace(go.Scatter(x=grid3, y=pi_true[k] * gaussian(grid3, mu_true[k], var_true[k]),
line=dict(color=c, dash='dash'), name=f'component {k}'))
fig.update_layout(title='Samples drawn by ancestor sampling, with the density that generated them',
width=880, height=470, template='plotly_white')
fig
The dashed curves are the weighted components and the solid curve is their sum.
Note what we needed in order to sample: the component $z$ of every point. When fitting, that is precisely what we do not have.
6. Implementing EM¶
We now fit a mixture without knowing $z$. The two steps from lecture translate directly into code.
The functions below are written for general $D$, taking data of shape $(N, D)$ and covariance matrices of shape $(K, D, D)$. Every example in this notebook uses $D = 1$, which is what the slides derive, but the same code fits multivariate mixtures unchanged.
E-step. With the parameters fixed, compute the responsibility of each component for each point, which is the posterior over $z$ obtained by Bayes' theorem:
$$\gamma_{nk} = \frac{\pi_k\,\mathcal{N}(x_n \mid \mu_k, \sigma_k^2)} {\sum_{k'} \pi_{k'}\,\mathcal{N}(x_n \mid \mu_{k'}, \sigma_{k'}^2)}$$
# Posterior over the latent component for every point: returns an (N, K) array.
# Written for general D; with D = 1 each Sigma[k] is a 1x1 matrix holding sigma_k^2.
def E_step(x, mu, Sigma, pi):
N, D = x.shape
K = len(pi)
gamma = np.zeros((N, K))
for k in range(K):
gamma[:, k] = pi[k] * multivariate_normal(mu[k], Sigma[k]).pdf(x)
return gamma / gamma.sum(axis=1, keepdims=True)
M-step. With the responsibilities fixed, maximize the expected complete-data log likelihood. Each parameter has a closed form, derived in lecture. In one dimension the last line reads $\sigma_k^2 = \frac{1}{N_k}\sum_n \gamma_{nk}(x_n - \mu_k)^2$, the same weighted average of squared deviations we derived for a single Gaussian.
$$N_k = \sum_n \gamma_{nk} \qquad \pi_k = \frac{N_k}{N} \qquad \mu_k = \frac{1}{N_k}\sum_n \gamma_{nk} x_n$$
# Closed-form parameter updates. `reg` keeps a component from collapsing onto one point.
def M_step(x, gamma, reg=1e-6):
N, D = x.shape
K = gamma.shape[1]
mu = np.zeros((K, D))
Sigma = np.zeros((K, D, D))
pi = np.zeros(K)
for k in range(K):
N_k = gamma[:, k].sum()
mu[k] = gamma[:, k] @ x / N_k
centred = x - mu[k]
Sigma[k] = (gamma[:, k] * centred.T) @ centred / N_k + reg * np.eye(D)
pi[k] = N_k / N
return mu, Sigma, pi
The reg term addresses the singularity from lecture. Without it a component can shrink onto a
single point, its variance collapsing toward zero and the likelihood diverging.
Initialization uses k-means, as sklearn does.
from sklearn.cluster import KMeans
def initialize(x, K, seed=0):
D = x.shape[1]
centers = KMeans(n_clusters=K, n_init=10, random_state=seed).fit(x).cluster_centers_
Sigma = np.array([np.atleast_2d(np.cov(x.T)) for _ in range(K)])
return centers, Sigma, np.ones(K) / K
def log_likelihood(x, mu, Sigma, pi):
per_component = np.column_stack(
[pi[k] * multivariate_normal(mu[k], Sigma[k]).pdf(x) for k in range(len(pi))])
return np.log(per_component.sum(axis=1)).sum()
def em(x, K, iters=50, seed=0):
mu, Sigma, pi = initialize(x, K, seed)
history = [log_likelihood(x, mu, Sigma, pi)]
for _ in range(iters):
gamma = E_step(x, mu, Sigma, pi)
mu, Sigma, pi = M_step(x, gamma)
history.append(log_likelihood(x, mu, Sigma, pi))
return mu, Sigma, pi, history
mu_em, Sigma_em, pi_em, history = em(X, K=2) # X is day.reshape(-1, 1)
for k in np.argsort(mu_em.ravel()):
print(f'component {k}: pi = {pi_em[k]:.3f} mu = {mu_em[k,0]:7.2f} '
f'sigma = {np.sqrt(Sigma_em[k,0,0]):.2f}')
print(f'\nlog likelihood: {history[0]:.1f} -> {history[-1]:.1f}')
component 1: pi = 0.690 mu = -51.92 sigma = 4.02 component 0: pi = 0.310 mu = -29.86 sigma = 4.56 log likelihood: -22775.7 -> -20041.1
Those are the same parameters scikit-learn found, recovered from an implementation that is about twenty lines long.
Each iteration increases the log likelihood¶
This is the property that guarantees convergence, and it is worth seeing rather than taking on faith.
fig = px.line(y=history, markers=True, labels={'x': 'iteration', 'y': 'log likelihood'},
title='EM increases the log likelihood at every step',
width=820, height=430, template='plotly_white')
fig.update_traces(line_color=BLUE)
fig
print('non-decreasing? ', bool(np.all(np.diff(history) >= -1e-6)))
print('largest single-step decrease:', np.diff(history).min())
print('\nours vs sklearn (means, sorted):',
np.round(np.sort(mu_em.ravel()), 3), np.round(np.sort(gmm.means_.ravel()), 3))
non-decreasing? True largest single-step decrease: -7.275957614183426e-12 ours vs sklearn (means, sorted): [-51.922 -29.864] [-51.921 -29.86 ]
Local optima¶
EM converges to a local maximum, so the starting point matters. Initializing at random rather than from k-means makes that visible.
def em_random_init(x, K, iters=50, seed=0):
g = np.random.default_rng(seed)
mu = x[g.choice(len(x), size=K, replace=False)] + g.normal(0, 8, size=(K, x.shape[1]))
Sigma = np.array([np.atleast_2d(np.cov(x.T)) for _ in range(K)])
pi = np.ones(K) / K
for _ in range(iters):
mu, Sigma, pi = M_step(x, E_step(x, mu, Sigma, pi))
return log_likelihood(x, mu, Sigma, pi)
random_finals = [em_random_init(X, 3, seed=s) for s in range(8)]
kmeans_finals = [em(X, K=3, iters=50, seed=s)[3][-1] for s in range(8)]
pd.DataFrame({'random init': np.round(random_finals, 1),
'k-means init': np.round(kmeans_finals, 1)})
| random init | k-means init | |
|---|---|---|
| 0 | -20043.4 | -20041.2 |
| 1 | -20043.6 | -20041.2 |
| 2 | -20041.0 | -20041.2 |
| 3 | -20043.4 | -20041.2 |
| 4 | -20041.2 | -20041.2 |
| 5 | -20041.2 | -20041.2 |
| 6 | -20043.3 | -20041.2 |
| 7 | -20043.6 | -20041.2 |
for name, vals in [('random init', random_finals), ('k-means init', kmeans_finals)]:
print(f'{name:13s} best {max(vals):10.1f} worst {min(vals):10.1f} '
f'spread {max(vals) - min(vals):8.1f}')
random init best -20041.0 worst -20043.6 spread 2.6 k-means init best -20041.2 worst -20041.2 spread 0.0
Fitting three components to data that really has two leaves genuine ambiguity about how to split it, and the random starts settle into two distinct optima a couple of nats apart. The k-means starts reach the same solution every time.
Two honest caveats. The gap here is small, because well-separated components in one dimension
make EM fairly robust; local optima bite much harder with more components and in higher
dimensions. And k-means initialization is not automatically better, only consistent: one random
start here found a marginally higher likelihood. What it buys is reproducibility, which is why it
is the default in sklearn.mixture.GaussianMixture, with n_init available for the cases where
consistency is not enough.
7. k-means is EM in a limit¶
Lecture claimed that constraining every component to the same fixed variance $\sigma_k^2 = \varepsilon$ and letting $\varepsilon \to 0$ turns the E-step into nearest-center assignment. We can watch the responsibilities harden as $\varepsilon$ shrinks.
# At small epsilon every density underflows to zero, so the responsibilities have to be
# computed in log space: subtract the row maximum before exponentiating.
def responsibilities_logspace(x, mu, Sigma, pi):
logs = np.column_stack([np.log(pi[k]) + multivariate_normal(mu[k], Sigma[k]).logpdf(x)
for k in range(len(pi))])
logs -= logs.max(axis=1, keepdims=True)
r = np.exp(logs)
return r / r.sum(axis=1, keepdims=True)
centers = KMeans(n_clusters=2, n_init=10, random_state=0).fit(X).cluster_centers_
pi_eq = np.ones(2) / 2
nearest = np.argmin(((X[:, None, :] - centers[None]) ** 2).sum(-1), axis=1)
rows = []
for eps in (100.0, 25.0, 5.0, 1.0, 0.1, 0.01):
Sig = np.array([eps * np.eye(1)] * 2)
g = responsibilities_logspace(X, centers, Sig, pi_eq)
rows.append({'epsilon': eps,
'mean largest responsibility': g.max(axis=1).mean(),
'fraction above 0.99': (g.max(axis=1) > 0.99).mean(),
'matches nearest center': (g.argmax(axis=1) == nearest).mean()})
pd.DataFrame(rows)
| epsilon | mean largest responsibility | fraction above 0.99 | matches nearest center | |
|---|---|---|---|---|
| 0 | 100.00 | 0.894111 | 0.009655 | 1.0 |
| 1 | 25.00 | 0.992857 | 0.915517 | 1.0 |
| 2 | 5.00 | 0.999410 | 0.994655 | 1.0 |
| 3 | 1.00 | 0.999923 | 0.999483 | 1.0 |
| 4 | 0.10 | 1.000000 | 1.000000 | 1.0 |
| 5 | 0.01 | 1.000000 | 1.000000 | 1.0 |
As $\varepsilon$ falls, the largest responsibility approaches 1 for every point: the soft assignment becomes a hard one.
The last column is 1.0 throughout, which is worth noticing. With equal weights and equal variances the most probable component is always the nearest center, whatever $\varepsilon$ is. What $\varepsilon$ controls is not which component wins but how confidently it wins.
k-means is therefore not a separate algorithm that resembles EM. It is EM for a mixture of equal, infinitely narrow Gaussians, which is also why it cannot represent components of different widths and a full mixture can.
Summary¶
| Gaussian MLE | the sample mean and sample variance, depending on the data only through $\sum_n x_n$ and $\sum_n x_n^2$ |
| Bias | $\sigma^2_{ML}$ is low by $(N-1)/N$, which is where ddof comes from |
| Model choice | maximum likelihood optimizes within a family and cannot rescue the wrong family |
| Mixtures | a weighted sum of Gaussians, with a latent $z$ naming the component |
| EM | alternate the posterior over $z$ with closed-form parameter updates; the log likelihood increases every step |
| k-means | EM for equal-variance components in the zero-variance limit |
Lecture 6 applies the Gaussian again, this time as a noise model, and shows that least squares is maximum likelihood.