---
title: "Learning with Hidden Variables: The EM Algorithm"
module: Learning
moduleNumber: 5
lessonNumber: 4
order: 504
summary: >
  Complete data can be learned by counting; real data usually hide some variables —
  the disease behind the symptoms, the cluster behind the points. This part develops
  the expectation-maximization algorithm, which learns those models by alternating an
  expected completion of the missing data with a re-estimation of the parameters. It
  works the idea through mixtures of Gaussians, Bayesian networks, and hidden Markov
  models, proves the monotone-likelihood guarantee from the evidence lower bound, and
  traces the line from EM to variational inference and the variational autoencoder.
topics: [Learning]
sources:
  - book: AIMA
    ref: "§20.3 Learning with Hidden Variables: The EM Algorithm"
---

This builds on [Learning Probabilistic Models](/artificial-intelligence/learning/probabilistic-learning),
which framed learning as Bayesian inference and derived the maximum-likelihood, MAP,
and Bayesian estimators for the case of **complete data** — where every example fixes
a value for every variable, so each parameter is just a count. Here we drop that
assumption. When some variables are never observed, counting is impossible, and the
learning problem needs the algorithm this part is about.

## Learning with hidden variables: the EM algorithm

Everything so far assumed complete data. Real data are usually **incomplete**: some
variables are never observed. A medical record lists symptoms, diagnosis, and
treatment, but rarely the underlying disease. Such a **hidden** (or **latent**)
variable is worth keeping precisely because it can _dramatically reduce_ the number
of parameters: a disease node with three predisposing causes and three symptoms may
need $78$ parameters, but splicing the disease out and connecting causes directly to
symptoms — now no longer conditionally independent — can balloon that to over $700$.
Hidden variables buy compactness at the cost of a harder learning problem: you
cannot count what you cannot see.

The **expectation–maximization** algorithm, EM, solves this in a general way. Its
one idea: _pretend you know the parameters, use them to infer a probability
distribution over the hidden variables, then refit the parameters as if that
inferred completion were real data_ — and iterate. Each pass alternates an
expectation step (complete the data in expectation) with a maximization step
(re-estimate the parameters).

### Unsupervised clustering with a mixture of Gaussians

The cleanest instance is **unsupervised clustering**: discovering categories in
unlabeled data. Model the data as generated by a **mixture distribution** with $k$ **components**.
A hidden variable $C$ names the component that generated a point,
and the density is

$$
P(\mathbf{x}) = \sum_{i=1}^{k} P(C = i)\, P(\mathbf{x} \mid C = i).
$$

For continuous data the natural component is a multivariate Gaussian, giving a
**mixture of Gaussians** with parameters $w_i = P(C = i)$ (the weight), $\mu_i$ (the
mean), and $\Sigma_i$ (the covariance) for each component $i$. If you _knew_ which
component generated each point you could fit each Gaussian directly; if you _knew_
the parameters you could assign each point to a component. When neither is known,
EM alternates between the two estimates.

Initialize the parameters arbitrarily, then iterate two steps. In the **E-step**,
compute the responsibility $p_{ij} = P(C = i \mid \mathbf{x}_j)$, the probability
that point $\mathbf{x}_j$ came from component $i$, by Bayes' rule from the current
parameters. In the **M-step**, refit each component to _all_ the data, weighting
each point by its responsibility.

$$
% caption: One EM pass on a mixture of Gaussians. The E-step assigns each point a
% soft responsibility to each cluster (shading); the M-step moves each cluster's
% mean and covariance to the responsibility-weighted data. Iterating tightens the
% components onto the true clusters.
\begin{tikzpicture}[>=stealth, font=\small,
  cl/.style={draw, ellipse, thick}]
  \definecolor{acc}{HTML}{2348F2}
  \definecolor{red}{HTML}{C0392B}
  % --- left: E-step ---
  \draw[black] (-0.3,-0.3) rectangle (3.3,3.3);
  \node[font=\footnotesize, anchor=south, black] at (1.5,3.35) {E-step: assign};
  \node[cl, acc, minimum width=13mm, minimum height=10mm] at (1.0,2.1) {};
  \node[cl, red, minimum width=13mm, minimum height=10mm] at (2.2,0.9) {};
  \foreach \p in {(0.7,2.3),(1.2,1.9),(0.9,2.5),(1.3,2.2)} \fill[acc] \p circle (1.6pt);
  \foreach \p in {(2.1,0.7),(2.4,1.0),(2.0,1.1),(2.5,0.7)} \fill[red] \p circle (1.6pt);
  % ambiguous points between clusters, half-shaded
  \fill[black] (1.6,1.5) circle (1.6pt);
  \fill[black] (1.5,1.2) circle (1.6pt);
  % --- arrow ---
  \draw[->, acc, very thick] (3.7,1.5) -- (4.7,1.5) node[midway, above, font=\scriptsize] {M-step};
  % --- right: M-step ---
  \begin{scope}[xshift=5.2cm]
    \draw[black] (-0.3,-0.3) rectangle (3.3,3.3);
    \node[font=\footnotesize, anchor=south, black] at (1.5,3.35) {M-step: update};
    \node[cl, acc, minimum width=11mm, minimum height=9mm] at (1.05,2.15) {};
    \node[cl, red, minimum width=11mm, minimum height=9mm] at (2.2,0.85) {};
    \foreach \p in {(0.7,2.3),(1.2,1.9),(0.9,2.5),(1.3,2.2),(1.0,2.0)} \fill[acc] \p circle (1.6pt);
    \foreach \p in {(2.1,0.7),(2.4,1.0),(2.0,1.1),(2.5,0.7),(2.2,0.9)} \fill[red] \p circle (1.6pt);
  \end{scope}
\end{tikzpicture}
$$

Writing $n_i = \sum_j p_{ij}$ for the effective number of points assigned to
component $i$, the M-step updates are responsibility-weighted versions of the
Gaussian ML estimators derived earlier:

$$
\mu_i \gets \frac{\sum_j p_{ij}\,\mathbf{x}_j}{n_i},
  \qquad
  \Sigma_i \gets \frac{\sum_j p_{ij}\,(\mathbf{x}_j - \mu_i)(\mathbf{x}_j - \mu_i)^\top}{n_i},
  \qquad
  w_i \gets \frac{n_i}{N}.
$$

The E-step computes the _expected_ values $p_{ij}$ of the hidden **indicator
variables** $Z_{ij}$ ($1$ if point $j$ came from component $i$, else $0$); the
M-step maximizes the log likelihood given those expectations. Applied to data
sampled from a three-component mixture, EM reconstructs a model nearly
indistinguishable from the generator. As pseudocode:

```algorithm
caption: $\textsc{EM-Mixture-of-Gaussians}$ — cluster unlabeled data into $k$ Gaussians
input: data $\mathbf{x}_1, \ldots, \mathbf{x}_N$, number of components $k$
initialize $w_i, \mu_i, \Sigma_i$ arbitrarily for $i = 1$ to $k$
repeat
  for each point $j$ and component $i$ do // E-step
    $p_{ij} \gets \alpha\, w_i\, P(\mathbf{x}_j \mid C = i)$ // responsibility, normalized over $i$
  for each component $i$ do // M-step
    $n_i \gets \sum_j p_{ij}$
    $\mu_i \gets \frac{1}{n_i} \sum_j p_{ij}\, \mathbf{x}_j$
    $\Sigma_i \gets \frac{1}{n_i} \sum_j p_{ij}\, (\mathbf{x}_j - \mu_i)(\mathbf{x}_j - \mu_i)^\top$
    $w_i \gets n_i / N$
until log likelihood $L$ converges
return $\{w_i, \mu_i, \Sigma_i\}$
```

#### One EM iteration, worked in numbers

Consider a single pass on four one-dimensional points, $x = \{1.0, 1.5, 5.0, 6.0\}$,
which plainly form two clusters. Fit a two-component mixture, and initialize
_deliberately badly_ to see EM correct itself: means $\mu_1 = 0.0$, $\mu_2 = 4.0$,
standard deviations $\sigma_1 = \sigma_2 = 2.0$, and equal weights $w_1 = w_2 =
0.5$.

**E-step.** For each point compute the responsibility $p_{1j} = P(C = 1 \mid x_j)$
by Bayes' rule, $p_{1j} = \dfrac{w_1\,\mathcal{N}(x_j; \mu_1, \sigma_1)}{w_1\,
\mathcal{N}(x_j; \mu_1, \sigma_1) + w_2\,\mathcal{N}(x_j; \mu_2, \sigma_2)}$. With
$\mathcal{N}(1.0; 0, 2) = 0.1760$ and $\mathcal{N}(1.0; 4, 2) = 0.0648$, the point
$x = 1.0$ gets $p_{11} = \tfrac{0.5(0.1760)}{0.5(0.1760) + 0.5(0.0648)} = 0.731$.
The full set of responsibilities to component 1:

$$
p_{1,\cdot} = (0.731,\; 0.622,\; 0.047,\; 0.018), \qquad
p_{2,\cdot} = (0.269,\; 0.378,\; 0.953,\; 0.982).
$$

The two low points lean toward component 1, the two high points toward component 2,
with the boundary points fractionally split — exactly the _soft_ assignment that
distinguishes EM from hard $k$-means.

**M-step.** The effective counts are $n_1 = \sum_j p_{1j} = 0.731 + 0.622 + 0.047 +
0.018 = 1.419$ and $n_2 = 4 - n_1 = 2.581$. The new means are
responsibility-weighted averages:

$$
\mu_1 \gets \frac{0.731(1.0) + 0.622(1.5) + 0.047(5.0) + 0.018(6.0)}{1.419} = 1.42,
\qquad
\mu_2 \gets \frac{\ldots}{2.581} = 4.45.
$$

The weighted variances give $\sigma_1 \gets 0.89$ and $\sigma_2 \gets 1.88$, and the
mixing weights become $w_1 \gets n_1/4 = 0.355$, $w_2 \gets 0.645$. In one step the
means have jumped from $(0.0, 4.0)$ toward the true cluster centers near $1.25$ and
$5.5$, and the data log likelihood rose from $-9.40$ to $-7.76$. That increase is
guaranteed on every iteration.

$$
% caption: One EM iteration on four 1-D points. Top: initial components (means at
% 0 and 4). The E-step assigns each point soft responsibilities (arrow thickness);
% the M-step pulls each mean to its responsibility-weighted data, so mu1 moves from
% 0 to 1.42 and mu2 from 4 to 4.45, raising log L from -9.40 to -7.76.
\begin{tikzpicture}[>=stealth, font=\footnotesize]
  \definecolor{acc}{HTML}{2348F2}
  \definecolor{red}{HTML}{C0392B}
  % number line
  \draw[black] (0,0) -- (7,0);
  \foreach \x/\lab in {0/0, 1.4/1.4, 4/4, 4.45/4.45} {}
  \foreach \v in {0,1,2,3,4,5,6} \draw[black] (\v,0.06) -- (\v,-0.06) node[anchor=north, font=\scriptsize] {\v};
  % data points
  \foreach \x in {1.0,1.5,5.0,6.0} \fill[black] (\x,0) circle (1.8pt);
  % initial means
  \draw[acc, thick] (0,0.1) -- (0,0.9); \node[acc, anchor=south, font=\scriptsize] at (0,0.9) {mu1 = 0};
  \draw[red, thick] (4,0.1) -- (4,0.9); \node[red, anchor=south, font=\scriptsize] at (4,0.9) {mu2 = 4};
  % updated means (dashed arrows to new position)
  \draw[acc, ->, dashed] (0,0.5) -- (1.4,0.5); \node[acc, anchor=south, font=\scriptsize] at (1.9,0.5) {1.42};
  \draw[red, ->, dashed] (4,-0.55) -- (4.45,-0.55); \node[red, anchor=north, font=\scriptsize] at (5.0,-0.5) {4.45};
  % LL annotation
  \node[anchor=west, font=\scriptsize] at (0,1.7) {log L: -9.40};
  \draw[->, black] (1.4,1.7) -- (2.3,1.7);
  \node[anchor=west, font=\scriptsize] at (2.4,1.7) {-7.76 after one step};
\end{tikzpicture}
$$

Two facts anchor EM's behavior. **It increases the log likelihood of the data at
every iteration** — this can be proved in general — and under mild conditions it
converges to a local maximum. It resembles gradient-based hill-climbing but has no
step-size parameter to tune.

$$
% caption: The log likelihood $L$ of the data as a function of EM iteration. It
% rises monotonically toward the value under the true model (dashed), climbing fast
% early and slowly near convergence — the typical EM trajectory.
\begin{tikzpicture}[>=stealth, font=\small]
  \definecolor{acc}{HTML}{2348F2}
  \draw[->, black] (0,0) -- (7.4,0);
  \node[black, font=\footnotesize, anchor=north] at (3.5,-0.55) {iteration};
  \draw[->, black] (0,0) -- (0,4.0) node[anchor=east, black, font=\footnotesize, rotate=90, yshift=6mm] {log L};
  % true-model level
  \draw[black, dashed] (0,3.5) -- (7,3.5);
  \node[black, anchor=west, font=\scriptsize] at (5.0,3.78) {true value};
  % rising curve, fast then slow
  \draw[acc, very thick] (0,0.6) .. controls (0.6,2.4) and (1.4,2.6) .. (2.4,2.9)
    .. controls (3.6,3.15) and (5,3.4) .. (7,3.48);
  \foreach \x/\y in {0/0.6, 0.8/2.5, 1.6/2.72, 2.4/2.9, 3.4/3.12, 4.6/3.32, 6/3.44}
    \fill[acc] (\x,\y) circle (1.4pt);
\end{tikzpicture}
$$

EM is not foolproof. A component can shrink onto a single point, sending its
variance to zero and its likelihood to infinity; two components can merge onto the
same data. These degenerate local maxima worsen in high dimensions. Placing priors
on the parameters (the MAP version of EM), restarting a collapsing component, and
sensible initialization all help.

### EM for Bayesian networks and HMMs

The same insight extends beyond mixtures. To learn a Bayesian network with hidden
variables, treat the E-step as computing, by ordinary Bayes-net inference, the
_expected counts_ that you would otherwise tabulate from complete data. The canonical
example is the **two-bag candy mixture**: two bags of candy have been tipped together,
so each piece has a $Flavor$ (cherry or lime), a $Wrapper$ (red or green), and either
a $Hole$ or not, but the $Bag$ it came from is hidden. Within a bag the three features
are independent — a naive-Bayes model — but the distribution of each feature depends on
the bag.

$$
% caption: The two-bag candy mixture as a Bayesian network. The hidden Bag at the
% root selects a bag; conditioned on it, Flavor, Wrapper, and Holes are independent
% (naive Bayes). Seven parameters: theta = P(Bag=1), and for each feature a
% probability given each bag. The CPT for Flavor is shown; Wrapper and Holes are
% analogous.
\begin{tikzpicture}[>=stealth, font=\small,
  var/.style={draw, ellipse, minimum width=17mm, minimum height=8mm, font=\footnotesize},
  cpt/.style={draw, black, font=\scriptsize, inner sep=3pt, align=left}]
  \definecolor{acc}{HTML}{2348F2}
  \node[var, draw=acc, text=acc, thick] (bag) at (0,2.0) {Bag};
  \node[var] (fl) at (-3.2,0) {Flavor};
  \node[var] (wr) at (0,0) {Wrapper};
  \node[var] (ho) at (3.2,0) {Holes};
  \draw[->, acc, thick] (bag) -- (fl);
  \draw[->, acc, thick] (bag) -- (wr);
  \draw[->, acc, thick] (bag) -- (ho);
  % Bag prior
  \node[cpt, anchor=west] at (0.9,2.3) {P(Bag=1) = t};
  % Flavor CPT
  \node[cpt, anchor=east] at (-4.5,-0.9)
    {Bag : P(cherry given Bag)\\ 1 : tF1\\ 2 : tF2};
\end{tikzpicture}
$$

For that model, with the bag hidden, the expected count of candies from bag $1$ is

$$
\hat N(Bag = 1) = \sum_{j=1}^{N} P(Bag = 1 \mid flavor_j, wrapper_j, holes_j),
$$

and each CPT parameter $\theta_{ijk} = P(X_i = x_{ij} \mid \mathbf{U}_i =
\mathbf{u}_{ik})$ is updated by the ratio of expected counts,

$$
\theta_{ijk} \gets \hat N(X_i = x_{ij}, \mathbf{U}_i = \mathbf{u}_{ik}) \, / \,
  \hat N(\mathbf{U}_i = \mathbf{u}_{ik}).
$$

#### One EM iteration on the candy bags, in numbers

The Gaussian mixture above was carried through with real numbers; the candy net
deserves the same. Generate $N = 1000$ candies from a true model in which the two bags
are equally likely, bag 1 is mostly cherry / red / holed and bag 2 mostly lime / green
/ no-hole:

$$
\theta = 0.5, \qquad \theta_{F1} = \theta_{W1} = \theta_{H1} = 0.8, \qquad
\theta_{F2} = \theta_{W2} = \theta_{H2} = 0.3.
$$

Because $Bag$ is hidden, the data are just counts of the eight observable candy types
$(Flavor, Wrapper, Holes)$, worked here with the sample counts and initialization
AIMA uses:[^candyem]

$$
% caption: The 1000 sampled candies by observable type. Columns are wrapper x holes,
% rows are flavor. Bag is not recorded. Red-wrapped cherries with a hole (273) are the
% most common type, reflecting bag 1's true profile; green no-hole limes (167) reflect
% bag 2.
\begin{tikzpicture}[>=stealth, font=\footnotesize]
  \definecolor{acc}{HTML}{2348F2}
  \definecolor{red}{HTML}{C0392B}
  % header rows
  \node[font=\scriptsize\bfseries, text=black] at (-1.6,2.1) {};
  \node[font=\scriptsize, text=black, anchor=south] at (1.4,2.15) {W = red};
  \node[font=\scriptsize, text=black, anchor=south] at (4.2,2.15) {W = green};
  \foreach \x/\lab in {0.7/{H=1}, 2.1/{H=0}, 3.5/{H=1}, 4.9/{H=0}}
    \node[font=\scriptsize, text=black] at (\x,1.65) {\lab};
  \node[font=\scriptsize, text=black, anchor=east] at (-0.15,1.05) {F = cherry};
  \node[font=\scriptsize, text=black, anchor=east] at (-0.15,0.35) {F = lime};
  % cell borders
  \draw[black] (0,0) rectangle (5.6,1.4);
  \draw[black] (1.4,0) -- (1.4,1.4);
  \draw[black] (2.8,0) -- (2.8,1.4);
  \draw[black] (4.2,0) -- (4.2,1.4);
  \draw[black] (0,0.7) -- (5.6,0.7);
  % cherry row
  \node[text=acc] at (0.7,1.05) {273};
  \node at (2.1,1.05) {93};
  \node at (3.5,1.05) {104};
  \node at (4.9,1.05) {90};
  % lime row
  \node at (0.7,0.35) {79};
  \node at (2.1,0.35) {100};
  \node at (3.5,0.35) {94};
  \node[text=red] at (4.9,0.35) {167};
\end{tikzpicture}
$$

Initialize the parameters deliberately off the truth, $\theta^{(0)} = 0.6$ and every
feature probability $0.6$ for bag 1 and $0.4$ for bag 2.

**E-step.** For each of the eight candy types, compute the responsibility $P(Bag = 1
\mid flavor, wrapper, holes)$ by Bayes' rule on the naive-Bayes model,

$$
P(Bag = 1 \mid f, w, h) =
  \frac{\theta \, P(f \mid B{=}1)\, P(w \mid B{=}1)\, P(h \mid B{=}1)}
       {\sum_{i} P(B{=}i)\, P(f \mid B{=}i)\, P(w \mid B{=}i)\, P(h \mid B{=}i)}.
$$

For the 273 red-wrapped cherry candies with holes, every factor favors bag 1, so

$$
P(B{=}1 \mid cherry, red, hole) =
  \frac{0.6 \cdot 0.6 \cdot 0.6 \cdot 0.6}{0.6 \cdot 0.6^3 + 0.4 \cdot 0.4^3}
  = \frac{0.1296}{0.1552} = 0.8351,
$$

and the eight responsibilities range from $0.835$ down to $0.308$ for the type that
best matches bag 2:

$$
% caption: E-step responsibilities P(Bag=1) for the eight candy types under the
% initial parameters, with each type's count. A cherry/red/hole candy is 0.835 likely
% to be bag 1; a lime/green/no-hole candy only 0.308. The types matching neither
% profile sit at 0.5. These are soft, not hard, assignments.
\begin{tikzpicture}[>=stealth, font=\footnotesize,
  row/.style={font=\scriptsize, anchor=west},
  val/.style={font=\scriptsize, anchor=east, text=acc}]
  \definecolor{acc}{HTML}{2348F2}
  \definecolor{red}{HTML}{C0392B}
  \node[font=\scriptsize\bfseries, text=black, anchor=west] at (0,3.5) {candy type (F, W, H)};
  \node[font=\scriptsize\bfseries, text=black, anchor=east] at (7.6,3.5) {count};
  \node[font=\scriptsize\bfseries, text=black, anchor=east] at (9.8,3.5) {P(Bag=1)};
  \foreach \i/\lab/\cnt/\p in {
    0/{cherry, red, hole}/273/0.835,
    1/{cherry, green, hole}/104/0.692,
    2/{cherry, red, no-hole}/93/0.692,
    3/{lime, red, hole}/79/0.692,
    4/{cherry, green, no-hole}/90/0.500,
    5/{lime, red, no-hole}/100/0.500,
    6/{lime, green, hole}/94/0.500,
    7/{lime, green, no-hole}/167/0.308} {
    \node[row] at (0, 3.0 - \i*0.38) {\lab};
    \node[font=\scriptsize, anchor=east] at (7.6, 3.0 - \i*0.38) {\cnt};
    \node[val] at (9.8, 3.0 - \i*0.38) {\p};
    % small bar for the probability
    \draw[acc!30, line width=3pt] (10.1, 3.0 - \i*0.38) -- ++(\p*2.2, 0);
  }
\end{tikzpicture}
$$

**M-step.** Sum the responsibilities over all candies to get the expected count of
bag-1 candies, then divide by $N$ for the new prior. Weighting each type's count by its
responsibility (the $273$ candies contribute $273 \times 0.835 = 227.97$, and so on
across all eight), the expected bag-1 total is $\hat N(Bag{=}1) = 612.4$, so

$$
\theta^{(1)} = \hat N(Bag{=}1)/N = 612.4 / 1000 = 0.6124.
$$

The feature parameters follow the same ratio-of-expected-counts recipe. For $\theta_{F1}
= P(cherry \mid Bag{=}1)$, sum the responsibilities of the four _cherry_ types for the
numerator and all eight for the denominator, $\theta_{F1}^{(1)} = \hat N(cherry,
Bag{=}1) / \hat N(Bag{=}1)$. Doing this for every parameter gives

$$
\theta^{(1)} = 0.6124, \quad
\theta_{F1}^{(1)} = 0.6684, \quad
\theta_{W1}^{(1)} = 0.6483, \quad
\theta_{H1}^{(1)} = 0.6558,
$$

$$
\theta_{F2}^{(1)} = 0.3887, \quad
\theta_{W2}^{(1)} = 0.3817, \quad
\theta_{H2}^{(1)} = 0.3827.
$$

Every bag-1 parameter has moved up from $0.6$ toward the true $0.8$, and every bag-2
parameter down from $0.4$ toward the true $0.3$ — one iteration already separates the
bags in the right direction. The log likelihood of the
1000 candies rises from about $-2044$ under $\theta^{(0)}$ to about $-2021$ after this
single step, an improvement in the likelihood itself by a factor of roughly $e^{23}
\approx 10^{10}$.

$$
% caption: The data log likelihood over EM iterations on the candy bags. One step
% lifts it from about -2044 to -2021; by the tenth iteration the learned model (about
% -1982) fits better than the generating model, after which progress slows to a crawl.
\begin{tikzpicture}[>=stealth, font=\footnotesize]
  \definecolor{acc}{HTML}{2348F2}
  \definecolor{red}{HTML}{C0392B}
  \draw[->, black] (0,0) -- (7.6,0);
  \node[black, font=\footnotesize, anchor=north] at (3.6,-0.5) {iteration};
  \draw[->, black] (0,0) -- (0,4.3);
  \node[black, font=\footnotesize, rotate=90, anchor=south] at (-1.05,2.1) {log L};
  \foreach \x/\l in {0/0, 2/2, 4/6, 6/10} \draw[black] (\x,0.05) -- (\x,-0.05) node[anchor=north, font=\scriptsize] {\l};
  \foreach \y/\l in {0.3/{-2044}, 2.0/{-2021}, 3.6/{-1982}} \draw[black] (0.05,\y) -- (-0.05,\y) node[anchor=east, font=\scriptsize] {\l};
  % rising curve
  \draw[acc, very thick] (0,0.3) .. controls (0.7,1.7) and (1.3,1.95) .. (2,2.0)
    .. controls (3.5,2.15) and (4,3.45) .. (6,3.6) -- (7,3.6);
  \foreach \x/\y in {0/0.3, 1/1.75, 2/2.0, 3/2.5, 4/3.2, 5/3.5, 6/3.6}
    \fill[acc] (\x,\y) circle (1.5pt);
  \draw[black, dashed] (0,3.5) -- (7,3.5);
  \node[black, anchor=east, font=\scriptsize] at (6.9,3.22) {true-model L};
  \node[acc, anchor=west, font=\scriptsize] at (0.15,0.55) {start};
\end{tikzpicture}
$$

Two features of the candy run mirror the general theory. The assignments are _soft_ —
the ambiguous types sit at exactly $0.5$ and contribute to both bags — which is why EM,
unlike a hard bag-by-bag count, can make progress with the bag entirely unobserved.
And each parameter update reads off _only_ a local posterior: $\theta_{F1}$ needs the
posterior over $Bag$ for each candy and nothing about $Wrapper$ or $Holes$ beyond what
that posterior already summarizes. This locality is what turns the update into a
by-product of ordinary Bayes-net inference.

The general lesson: **the parameter updates for Bayesian-network learning with
hidden variables are read directly off the results of inference on each example**,
and only _local_ posteriors — over each variable and its parents — are needed.
Exact inference algorithms such as variable elimination produce these as a
by-product, with no learning-specific computation.

The final case is the **hidden Markov model**. An HMM is a dynamic Bayesian network
with one discrete state variable, and each data point is an observation _sequence_.
The complication over general Bayes nets is that the transition probabilities
$\theta_{ijt} = P(X_{t+1} = j \mid X_t = i)$ are shared across all time steps
($\theta_{ijt} = \theta_{ij}$ for all $t$), so the update sums expected transition
counts over time,

$$
\theta_{ij} \gets \frac{\sum_t \hat N(X_{t+1} = j,\, X_t = i)}{\sum_t \hat N(X_t = i)}.
$$

The expected counts come from the **forward–backward** algorithm, run as
_smoothing_ rather than filtering — you must consider later evidence to estimate
that a transition occurred. This instance of EM is the **Baum–Welch** algorithm.

### The general form

Strip the examples away and one equation is left. Let $\mathbf{x}$ be all observed
values, $\mathbf{Z}$ all hidden variables, and $\boldsymbol{\theta}$ all parameters.
EM iterates

$$
\boldsymbol{\theta}^{(i+1)} = \arg\max_{\boldsymbol{\theta}}
  \sum_{\mathbf{z}} P(\mathbf{Z} = \mathbf{z} \mid \mathbf{x}, \boldsymbol{\theta}^{(i)})\,
  L(\mathbf{x}, \mathbf{Z} = \mathbf{z} \mid \boldsymbol{\theta}).
$$

The E-step is the summation: the expected complete-data log likelihood under the
posterior $P(\mathbf{Z} \mid \mathbf{x}, \boldsymbol{\theta}^{(i)})$ over the hidden
variables. The M-step is the maximization of that expected log likelihood over the
parameters. Every instance is this equation with $\mathbf{Z}$ specialized: the
component indicator $Z_{ij}$ for mixtures, an unobserved variable's value for Bayes
nets, the sequence state for HMMs.

```algorithm
caption: $\textsc{Expectation-Maximization}$ — parameter learning with hidden variables
input: observed data $\mathbf{x}$, model with hidden variables $\mathbf{Z}$
initialize parameters $\boldsymbol{\theta}^{(0)}$
$i \gets 0$
repeat
  compute the posterior $P(\mathbf{Z} \mid \mathbf{x}, \boldsymbol{\theta}^{(i)})$ // E-step
  $Q(\boldsymbol{\theta}) \gets \sum_{\mathbf{z}} P(\mathbf{z} \mid \mathbf{x}, \boldsymbol{\theta}^{(i)})\, L(\mathbf{x}, \mathbf{z} \mid \boldsymbol{\theta})$
  $\boldsymbol{\theta}^{(i+1)} \gets \arg\max_{\boldsymbol{\theta}} Q(\boldsymbol{\theta})$ // M-step
  $i \gets i + 1$
until $L(\mathbf{x} \mid \boldsymbol{\theta}^{(i)})$ converges
return $\boldsymbol{\theta}^{(i)}$
```

### Why EM climbs: the evidence lower bound

The monotone-increase guarantee follows from one inequality. The
quantity EM actually wants to raise is the log likelihood of the observed data,
$\log P(\mathbf{x} \mid \boldsymbol{\theta}) = \log \sum_{\mathbf{z}} P(\mathbf{x},
\mathbf{z} \mid \boldsymbol{\theta})$, a log-of-a-sum that is hard to optimize
directly because the hidden $\mathbf{z}$ is buried inside. Introduce any
distribution $q(\mathbf{z})$ over the hidden variables and rewrite, then apply
**Jensen's inequality** (the log of an average is at least the average of the logs,
since $\log$ is concave):

$$
\log P(\mathbf{x} \mid \boldsymbol{\theta})
  = \log \sum_{\mathbf{z}} q(\mathbf{z})\, \frac{P(\mathbf{x}, \mathbf{z} \mid \boldsymbol{\theta})}{q(\mathbf{z})}
  \;\ge\; \sum_{\mathbf{z}} q(\mathbf{z}) \log \frac{P(\mathbf{x}, \mathbf{z} \mid \boldsymbol{\theta})}{q(\mathbf{z})}
  \;=\; \mathcal{L}(q, \boldsymbol{\theta}).
$$

The right-hand side $\mathcal{L}(q, \boldsymbol{\theta})$ is the **evidence lower
bound** (ELBO). The gap between it and the true log likelihood coincides with the
Kullback–Leibler divergence $\text{KL}\big(q(\mathbf{z}) \,\|\, P(\mathbf{z} \mid
\mathbf{x}, \boldsymbol{\theta})\big) \ge 0$, so the bound is tight — equality holds
— precisely when $q(\mathbf{z}) = P(\mathbf{z} \mid \mathbf{x}, \boldsymbol{\theta})$,
the posterior over the hidden variables.

EM is coordinate ascent on this bound. The **E-step** maximizes $\mathcal{L}$ over
$q$ with $\boldsymbol{\theta}$ fixed, and the maximizer is the posterior
$q(\mathbf{z}) = P(\mathbf{z} \mid \mathbf{x}, \boldsymbol{\theta}^{(i)})$ — which is
why the E-step computes exactly that posterior, and after it the bound touches the
log likelihood. The **M-step** then maximizes $\mathcal{L}$ over
$\boldsymbol{\theta}$ with $q$ fixed, which (dropping the $\boldsymbol{\theta}$-free
entropy of $q$) is the expected complete-data log likelihood $Q(\boldsymbol{\theta})$
maximized earlier. Because each step never lowers $\mathcal{L}$ and the E-step makes
$\mathcal{L}$ equal the log likelihood, the log likelihood itself cannot decrease:

$$
\log P(\mathbf{x} \mid \boldsymbol{\theta}^{(i+1)})
  \;\ge\; \mathcal{L}(q^{(i)}, \boldsymbol{\theta}^{(i+1)})
  \;\ge\; \mathcal{L}(q^{(i)}, \boldsymbol{\theta}^{(i)})
  \;=\; \log P(\mathbf{x} \mid \boldsymbol{\theta}^{(i)}).
$$

That chain is the whole convergence proof. It also explains the failure modes: EM
maximizes a bound that only touches the true objective at the current parameters, so
it can settle at any local maximum, and where the likelihood is unbounded (a Gaussian
component collapsing onto one point) EM will climb toward infinity.

$$
% caption: EM as coordinate ascent on the evidence lower bound. The E-step raises the
% bound L until it touches the log likelihood at the current theta (KL = 0); the
% M-step slides theta to the bound's maximum, raising the true log likelihood; the
% next E-step re-tightens the bound at the new theta.
\begin{tikzpicture}[>=stealth, font=\footnotesize]
  \definecolor{acc}{HTML}{2348F2}
  \definecolor{red}{HTML}{C0392B}
  \draw[->, black] (0,0) -- (7.4,0) node[right, font=\scriptsize] {theta};
  \draw[->, black] (0,0) -- (0,4.0) node[above, font=\scriptsize] {value};
  % true log-likelihood: broad concave-ish hump
  \draw[black, very thick] (0.3,1.2) .. controls (2.5,3.6) and (4.5,3.7) .. (7.0,2.6);
  \node[anchor=west, font=\scriptsize] at (5.4,3.4) {log P(x $\mid$ theta)};
  % lower bound at theta_i, touches at x=1.5
  \draw[acc, thick] (0.4,0.5) .. controls (1.5,2.2) and (2.2,2.2) .. (3.0,0.9);
  \node[acc, anchor=west, font=\scriptsize] at (4.4,0.7) {bound at i};
  \fill[black] (1.5,2.16) circle (1.6pt);
  \draw[black, dotted] (1.5,0) -- (1.5,2.16); \node[anchor=north, font=\scriptsize] at (1.5,0) {i};
  % M-step moves to peak of bound at ~1.85; new theta
  \draw[red, dotted] (1.85,0) -- (1.85,2.22);
  \draw[->, red, thick] (1.5,1.1) -- (1.85,1.1); \node[red, anchor=south, font=\scriptsize] at (1.68,1.1) {M};
  % next bound at theta_{i+1}, touches higher
  \draw[acc, thick, dashed] (2.2,0.6) .. controls (3.4,3.0) and (4.2,3.0) .. (5.2,1.1);
  \fill[black] (3.6,2.98) circle (1.6pt);
  \node[anchor=north, font=\scriptsize] at (3.7,0) {i+1};
  \draw[black, dotted] (3.6,0) -- (3.6,2.98);
\end{tikzpicture}
$$

Once the general form is understood, variants follow. When the exact E-step is
intractable — large Bayes nets — an _approximate_ E-step still yields effective
learning: with an MCMC sampler, each configuration of hidden and observed variables
visited is treated as a complete observation, and parameters update after each
transition. Variational and loopy methods serve the same role for very large
networks. The same latent-variable-and-EM story reappears in the deep-learning
treatment of
[latent-variable models](/deep-learning/generative-models/variational-autoencoders),
where a neural network parameterizes the components and the E-step becomes an
_approximate_ inference network rather than exact Bayes' rule — the variational
autoencoder is EM with the posterior itself learned.

## From EM to modern probabilistic learning

The candy bags and mixtures of Gaussians in this lesson are the classical core of
statistical learning; the same machinery scales into much of modern machine
learning, and the public literature fills in the pieces AIMA only sketches.

**EM's own history and generality.** The algorithm was named and given its general
form by Dempster, Laird, and Rubin (1977), who proved the monotone-likelihood result
via the lower-bound argument above; the coordinate-ascent-on-the-ELBO view was made
explicit by Neal and Hinton (1998), which is what licenses _partial_ or _incremental_
E-steps that update only some responsibilities per pass and still converge.[^em-dlr]
The Baum–Welch algorithm for HMMs actually predates the general EM formulation (Baum
et al., 1970), a special case discovered first.[^em-baum]

**Conjugate priors and the exponential family.** The beta–binomial and
Dirichlet–multinomial conjugacies used here are instances of a general fact: every
distribution in the **exponential family** has a conjugate prior, and Bayesian
updating reduces to adding sufficient statistics to the prior's hyperparameters. This
structure, laid out in Bishop's _Pattern Recognition and Machine Learning_ (2006), is
what makes tractable Bayesian inference possible at all, and it is the backbone of
**latent Dirichlet allocation** (Blei, Ng, and Jordan, 2003), the topic model that
applies exactly the mixture-and-EM idea of this lesson to documents, with words drawn
from latent topics.[^lda]

**When EM is intractable: variational inference and MCMC.** For models where the
exact posterior in the E-step cannot be computed, two families of approximation
dominate. **Variational inference** replaces the true posterior with the closest
member of a tractable family by maximizing the same ELBO derived above — turning
inference into optimization (Jordan et al., 1999; Blei, Kucukelbir, and McAuliffe,
2017).[^vi] Markov-chain Monte-Carlo instead draws samples from the posterior; Gelman
et al.'s _Bayesian Data Analysis_ is the standard reference. The **variational
autoencoder** (Kingma and Welling, 2013) is the meeting point of both threads and
this lesson: it is EM in which a neural network amortizes the E-step, learning a
single inference network $q_\phi(\mathbf{z} \mid \mathbf{x})$ that approximates the
posterior for every data point at once, trained by maximizing the ELBO with the
_reparameterization trick_.[^vae] Read that way, the latent-variable models behind
modern generative AI are direct descendants of the mixture-of-Gaussians EM
worked here.

**Learning Bayes-net structure, in practice.** The score-and-search sketch of this
lesson corresponds to real algorithms: the **BIC/MDL score** penalizes structure by
a term proportional to $\tfrac{1}{2}(\log N)\,|\boldsymbol{\theta}|$, and the
constraint-based **PC algorithm** (Spirtes, Glymour, and Scheines, 1993) recovers
structure from conditional-independence tests, a starting point for the modern field
of causal discovery.[^struct]

EM is the standard method for learning when data are incomplete. Whenever a model
has variables you cannot observe — clusters without labels, diseases behind
symptoms, states behind an emission sequence — the procedure is the same: complete the
data in expectation, maximize as if it were real, repeat.
[^candyem]: **AIMA**, §20.3.1 — Unsupervised clustering: learning mixtures of Gaussians, and the two-bag naive-Bayes candy mixture worked through one EM iteration. The 1000-sample counts, the initialization $\theta^{(0)} = 0.6$, the E-step responsibility of the 273 red-holed cherries ($\approx 0.835$), the resulting $\theta^{(1)} = 0.6124$ and the other updated parameters, and the log-likelihood rise from $\approx -2044$ to $\approx -2021$ (and to $\approx -1982$ by the tenth iteration) are all from AIMA's Figure 20.13 example.
[^em-dlr]: Dempster, A. P., Laird, N. M., and Rubin, D. B. (1977), "Maximum Likelihood from Incomplete Data via the EM Algorithm," _Journal of the Royal Statistical Society, Series B_ 39(1): 1–38 — the general formulation and monotone-likelihood proof; Neal, R., and Hinton, G. (1998), "A View of the EM Algorithm that Justifies Incremental, Sparse, and Other Variants," gives the ELBO coordinate-ascent view.
[^em-baum]: Baum, L. E., Petrie, T., Soules, G., and Weiss, N. (1970), "A Maximization Technique Occurring in the Statistical Analysis of Probabilistic Functions of Markov Chains," _Annals of Mathematical Statistics_ 41(1): 164–171 — the forward–backward (Baum–Welch) procedure for HMMs, a special case of EM found before the general algorithm.
[^lda]: Bishop, C. M. (2006), _Pattern Recognition and Machine Learning_, Springer — exponential families and conjugate priors; Blei, D., Ng, A., and Jordan, M. (2003), "Latent Dirichlet Allocation," _Journal of Machine Learning Research_ 3: 993–1022 — a mixture model over documents fit by variational EM.
[^vi]: Jordan, M. I., Ghahramani, Z., Jaakkola, T. S., and Saul, L. K. (1999), "An Introduction to Variational Methods for Graphical Models," _Machine Learning_ 37: 183–233; Blei, D., Kucukelbir, A., and McAuliffe, J. (2017), "Variational Inference: A Review for Statisticians," _JASA_ 112: 859–877 — approximate posterior inference by maximizing the ELBO over a tractable family.
[^vae]: Kingma, D. P., and Welling, M. (2013), "Auto-Encoding Variational Bayes," _ICLR 2014_ (arXiv:1312.6114) — the variational autoencoder, amortizing the E-step with a learned inference network and the reparameterization trick.
[^struct]: Spirtes, P., Glymour, C., and Scheines, R. (1993), _Causation, Prediction, and Search_, Springer — the constraint-based PC algorithm for recovering network structure from conditional-independence tests; the BIC score is from Schwarz, G. (1978), "Estimating the Dimension of a Model," _Annals of Statistics_ 6: 461–464.