Learning the Numbers in a Probability Model

Every probability model discussed so far arrived with its numbers already in place. A Bayesian network came with its conditional probability tables filled in; a Gaussian came with a mean and a variance. Somebody had to put them there. This article is about how, and it splits cleanly in two: the case where every variable is observed, which is almost trivially easy, and the case where one is not, which requires a genuinely different algorithm.

A. The likelihood, and the three steps

Fix a model with parameters $\theta$. The likelihood of a data set is the probability the model assigns to exactly the observations you got, which for independent examples is a product:

$$P(\mathbf{d} \mid h_\theta) = \prod_{j=1}^{N} P(d_j \mid h_\theta).$$

Maximum-likelihood estimation picks the $\theta$ maximising this. The recipe is always the same three steps: write the likelihood, differentiate its logarithm, set the derivative to zero.

Unwrap 25 candies from a bag with unknown cherry proportion $\theta$ and find 18 cherry. Then $P(\mathbf{d} \mid h_\theta) = \theta^{18}(1-\theta)^{7}$, so

$$L = 18\log\theta + 7\log(1-\theta), \qquad \frac{dL}{d\theta} = \frac{18}{\theta} - \frac{7}{1-\theta} = 0 \;\Rightarrow\; \theta = \frac{18}{25} = 0.72 .$$

A grid over two million values of $\theta$ agrees to six decimals: the log likelihood is $-14.823833$ at $0.72$, against $-14.847959$ at $0.70$ and $-14.849407$ at $0.74$.

The answer is the observed proportion, which is what anyone would have guessed. The value of the derivation is that it proves the guess is the maximiser, and that the same three steps keep working where no guess suggests itself.

Taking the logarithm is worth defending twice. Mathematically it cannot move a maximum, since $\log$ is strictly increasing. Computationally, a product of thousands of probabilities underflows to zero in floating point and a sum of logarithms does not - which is why production code lives in log space.

B. Why complete data makes a network cheap

Give the candies a wrapper too, red or green, with a probability depending on flavour. Three parameters now: $\theta$, $\theta_1 = P(\text{red} \mid \text{cherry})$, $\theta_2 = P(\text{red} \mid \text{lime})$. With counts 12 cherry-red, 6 cherry-green, 2 lime-red, 5 lime-green,

$$L = 12\log(\theta\theta_1) + 6\log(\theta(1-\theta_1)) + 2\log((1-\theta)\theta_2) + 5\log((1-\theta)(1-\theta_2)).$$

Expand the logarithms and every term containing $\theta_1$ contains nothing else. The partial derivatives are independent equations, each the one-parameter problem again:

$$\theta = 0.72, \qquad \theta_1 = \tfrac{12}{18} = 0.666667, \qquad \theta_2 = \tfrac{2}{7} = 0.285714 .$$

Nelder–Mead on the full three-parameter likelihood lands within $1.4 \times 10^{-8}$ of these, at a log likelihood of $-30.468975$.

That decoupling is the whole reason learning a Bayesian network from complete data is fast. Each conditional probability table is fitted by counting within its own conditioning case, with no search and no iteration - a set of independent counting problems rather than one joint optimisation.

The figure opens with the two wrapper sliders on the rounded values 0.67 and 0.29, so its total log likelihood reads $-30.469740$. Press snap to the counted estimates and the sliders move to the maximum-likelihood values, and the total to the $-30.468975$ quoted above.

C. The continuous case, and why least squares

For $N$ draws from a Gaussian, setting the two partial derivatives to zero gives

$$\mu = \frac{\sum_j x_j}{N}, \qquad \sigma = \sqrt{\frac{\sum_j (x_j - \mu)^2}{N}} .$$

Note the denominator. On ten measurements summing to $51.0$ with squared deviations summing to $3.98$, that is $\mu = 5.1$ and $\sigma = 0.630872$, where dividing by $N-1$ gives $0.664998$ - a factor of $\sqrt{N/(N-1)} = 1.054093$ apart. Both are correct answers to different questions: one maximises the likelihood, the other is unbiased. Maximum likelihood was never asked to be unbiased.

A useful corollary falls out. Fit $y$ as normal about $\theta_1 x + \theta_2$ with fixed variance, and the only parameter-dependent part of the log likelihood is $-\sum_j (y_j - \theta_1 x_j - \theta_2)^2$. Minimising squared error is maximum likelihood under Gaussian noise, which is the reason least squares is the default and not merely a convenient tradition.

D. What it does to an event it has not seen

Unwrap one candy, find it cherry, and the recipe reports $\theta = 1$: lime is impossible, probability exactly zero, a claim no later evidence can multiply its way out of. Maximum likelihood reads "not yet observed" as "cannot happen".

The Bayesian repair is a prior. For $\theta \in [0,1]$ the natural choice is a Beta distribution, whose useful property is conjugacy: a $\text{Beta}(a,b)$ prior with $c$ cherries and $\ell$ limes gives a $\text{Beta}(a+c, b+\ell)$ posterior. Updating is addition, and the posterior stays in the same family, so belief can be updated forever without the representation growing.

$\text{Beta}(a,b)$ behaves like $a+b-2$ imaginary observations already counted. Watch $\text{Beta}(2,2)$ - one imaginary candy of each flavour - handle the pathological case:

The prior does its work early and then gets out of the way: it is the difference between $0.6$ and a claim of certainty at one observation, and worth three and a half thousandths at 250. Note also that maximum likelihood is exactly MAP under a uniform prior, and that adding one to every count - the standard field hack - is the posterior mean under $\text{Beta}(1,1)$. It was a prior all along.

E. Naive Bayes, and one zero

A naive Bayes model assumes the attributes are conditionally independent given the class, so $P(C, x_1, \ldots, x_n) = P(C)\prod_i P(x_i \mid C)$: for $n$ Boolean attributes, $2n+1$ parameters, all fitted by counting.

Take 14 papers labelled theory or applied, described by whether each contains a proof, uses a dataset, and reports a benchmark. Six are theory, eight applied, so the priors are $3/7$ and $4/7$:

A paper with a proof and nothing else scores $5/21$ against $1/448$, a posterior of $320/323 = 0.990712$. Sensible.

Now look at the zero. No theory paper reports a benchmark, so $P(\text{benchmark} = 1 \mid \text{theory}) = 0$. Classify a paper with a proof, no dataset, and a benchmark: the theory score is a product containing that zero, so it is zero exactly. The posterior is $P(\text{theory}) = 0$ regardless of the proof, regardless of the prior. One attribute value absent from the training data holds a veto over the entire prediction - and in text classification, where most words are missing from most classes, that is the normal condition rather than an edge case.

Add-one smoothing turns $0/6$ into $1/8$, and the same paper scores $105/4096$ against $27/1000$: a posterior of $4375/8983 = 0.487031$ against $0.512969$. An honest near-tie instead of a certainty, and the smoothed classifier still labels all 14 training papers correctly.

The independence assumption is of course false - papers with benchmarks tend to have datasets - and wherever attributes are dependent within a class the shared evidence is counted twice, so the probabilities are not to be trusted. The rankings survive anyway, because the decision is an $\arg\max$ and double-counting tends to inflate both scores together. Trust the class, not the number.

F. When the counts cannot be taken

Everything above needed complete data. Hide a variable - which cluster a point came from, which topic wrote a document - and the likelihood of an observed example becomes a sum over what the hidden variable might have been:

$$P(x \mid \theta) = \sum_{z} P(x, Z = z \mid \theta).$$

The logarithm of a sum does not split, so the parameters no longer decouple and there is no closed form to set to zero.

Expectation–maximization substitutes expected counts for the counts it cannot take:

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

The E step computes the posterior over the hidden variables under the current parameters - the fractional membership a point gets is its responsibility. The M step applies the complete-data formulas with every count replaced by a sum of responsibilities. For a mixture of Gaussians that is

$$r_{jk} = \frac{w_k \mathcal{N}(x_j \mid \mu_k, \sigma_k)}{\sum_{k'} w_{k'} \mathcal{N}(x_j \mid \mu_{k'}, \sigma_{k'})}, \qquad w_k = \frac{n_k}{N}, \quad \mu_k = \frac{\sum_j r_{jk}x_j}{n_k}, \quad \sigma_k = \sqrt{\frac{\sum_j r_{jk}(x_j - \mu_k)^2}{n_k}},$$

with $n_k = \sum_j r_{jk}$. These are the Gaussian formulas from section C, weighted.

G. A measured run, and where k-means loses

Twenty points on a line: eight from a narrow component near zero ($-0.8, -0.7, -0.2, 0.0, 0.1, 0.5, 1.1, 1.3$) and twelve from a much broader one near five ($1.7, 3.0, 3.3, 3.4, 4.0, 4.2, 4.6, 4.8, 6.9, 7.1, 7.2, 9.2$), with the labels withheld.

From a crude start - means 0 and 5, both standard deviations 1 - the log likelihood rises from $-56.616239$ to $-46.633131$ over 82 iterations, never once falling, and settles on

$$w = (0.347288, 0.652712), \quad \mu = (0.128189, 4.581627), \quad \sigma = (0.731077, 2.367312).$$

The larger responsibility per point recovers all twenty withheld labels.

The figure stops at iteration 81, one short of the 82 above, because its stopping tolerance differs; the log likelihood it starts and ends on, $-56.616239$ and $-46.633131$, agrees with the run described here.

The monotone rise is guaranteed: each iteration maximises a lower bound touching the log likelihood at the current parameters, so a decrease means a bug. The global maximum is not guaranteed - over 400 random restarts, 385 reached $-46.633131$ and 15 stopped at $-48.277387$, a worse fixed point no further iteration escapes.

Now compare k-means on the same points. Its optimal split - verified by enumerating all 19 contiguous partitions, so this is not a local minimum - has centres $0.333333$ and $5.245455$ with a within-cluster sum of squares of $47.347273$, and it puts $x = 1.7$ in the narrow cluster. That is one point away from EM's answer and one point away from the truth.

The reason is instructive. $x = 1.7$ lies $1.571811$ from the fitted narrow mean and $2.881627$ from the broad one, so distance alone hands it over, and distance is all k-means has. EM compares densities:

$$0.347288 \times \mathcal{N}(1.7 \mid 0.128189, 0.731077) = 0.018788, \qquad 0.652712 \times \mathcal{N}(1.7 \mid 4.581627, 2.367312) = 0.052436,$$

giving responsibilities $0.263789$ and $0.736211$. The broad component is 3.24 times wider and carries the larger weight, and both count. The point is further from that component and still more likely to have come from it.

Responsibilities elsewhere are mostly decisive but not everywhere: $x = 1.1$ splits $0.677 / 0.323$, and $x = 1.3$ splits $0.555 / 0.445$. EM fits component widths and mixing weights, which k-means does not represent at all.

Key takeaways

• With complete data, maximum likelihood is three steps and the answer is a ratio of counts; the parameters of a network decouple, so fitting it is a set of independent counting problems.
• The maximum-likelihood Gaussian divides by $N$, not $N-1$, and it explains why least squares is the right loss under Gaussian noise.
• Maximum likelihood assigns probability zero to anything unobserved, which destroys a naive Bayes prediction through a single zero factor. A Beta prior is the principled repair, and add-one smoothing is a special case of it.
• Hidden variables put a sum inside the logarithm and break the decoupling. EM replaces counts with expected counts, never decreases the likelihood, and reaches only a local maximum - so restart it.
• Because EM models width and weight rather than distance alone, it recovers structures that hard assignment cannot.

What's next

Every EM instance here fitted a static model. The same two steps fit a model that changes over time - the E step becomes forward–backward smoothing, and the reason it must be smoothing rather than filtering is that estimating a transition needs evidence from after it happened.