The last lecture of Phase A, and the one that ties it together. Maximum likelihood is a single principle from which squared error, cross-entropy, ridge and LASSO all fall out — and Bayes' theorem is the reason a 99%-accurate test can be wrong most of the time.
For events \(A,B\) with \(P(B)>0\),
\[ P(A\mid B) = \frac{P(A\cap B)}{P(B)}. \]Read: having learned that \(B\) happened, restrict attention to \(B\) and renormalise. The denominator is what makes the conditional probabilities of a partition sum to 1 again.
Independence is the statement that conditioning changes nothing: \(P(A\mid B)=P(A)\), equivalently \(P(A\cap B)=P(A)P(B)\).
The lesson generalises well beyond medicine, and it is the single most useful piece of probability for an ML practitioner to hold on to. A classifier's accuracy on the positive class tells you almost nothing about the reliability of a positive prediction, when the positive class is rare.
Fraud detection, rare-disease screening, anomaly detection, safety filters — all live in this regime. Precision \(P(\text{actually positive}\mid\text{predicted positive})\) depends on the base rate; recall \(P(\text{predicted positive}\mid\text{actually positive})\) does not. That asymmetry is why you should report both, why the precision-recall curve is preferred to the ROC curve for imbalanced problems, and why a model with 99.9% accuracy on a 0.1%-prevalence problem may be worthless — predicting "negative" always would score 99.9%.
Four, with their means, variances, and what they model.
\[\begin{array}{llll} \textbf{Bernoulli}(p) & P(X{=}1)=p & \E=p & \Var=p(1-p) \\[3pt] \textbf{Categorical}(\boldsymbol\pi) & P(X{=}k)=\pi_k & - & - \\[3pt] \textbf{Gaussian}(\mu,\sigma^2) & f(x)=\frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{(x-\mu)^2}{2\sigma^2}} & \E=\mu & \Var=\sigma^2\\[3pt] \textbf{Multivariate Gaussian} & f(\vx)=\frac{1}{\sqrt{(2\pi)^d|\mS|}}e^{-\frac12(\vx-\vm)^\top\mS^{-1}(\vx-\vm)} & \E=\vm & \Cov=\mS \end{array}\]
Look at the exponent of the multivariate Gaussian: \((\vx-\vm)^\top\mS^{-1}(\vx-\vm)\) is a quadratic form — exactly Lecture 11's object. Its level sets are the ellipsoids we drew there, with axes along the eigenvectors of \(\mS\) and half-lengths proportional to \(\sqrt{\lambda_i}\). The density is constant on those ellipsoids. Everything you learned about quadratic forms is a fact about Gaussians.
Given data \(\mathcal{D}=\{x_1,\dots,x_n\}\) modelled as independent draws from \(p(x\mid\vt)\), the likelihood is
\[ L(\vt) = p(\mathcal{D}\mid\vt) = \prod_{i=1}^{n} p(x_i\mid\vt), \]and the maximum likelihood estimate is \(\hat{\vt}_{\text{MLE}}=\argmax_{\vt}L(\vt)\).
Crucially, \(L\) is a function of \(\vt\), not of the data. The data are fixed and observed; we are asking which parameter setting makes what we saw most probable. This is the single most common conceptual slip with likelihoods, and \(L(\vt)\) is not a probability distribution over \(\vt\) — it does not integrate to 1.
By Lecture 2's argmax theorem we may take logs, converting the product into a sum. The recipe is always the same four steps:
Data \(x_1,\dots,x_n\) drawn independently from \(\mathcal{N}(\mu,\sigma^2)\), with both parameters unknown. Find \(\hat\mu\) and \(\hat\sigma^2\).
Hint: write the log-likelihood, then differentiate with respect to each parameter in turn.
Log-likelihood.
\[ \ell(\mu,\sigma^2) = \sum_{i=1}^n \log\!\left[\frac{1}{\sqrt{2\pi\sigma^2}} e^{-\frac{(x_i-\mu)^2}{2\sigma^2}}\right] = -\frac{n}{2}\log(2\pi) - \frac{n}{2}\log\sigma^2 - \frac{1}{2\sigma^2}\sum_i (x_i-\mu)^2. \]For \(\mu\): only the last term involves it.
\[\frac{\partial\ell}{\partial\mu} = \frac{1}{\sigma^2}\sum_i (x_i-\mu) = 0 \;\Longrightarrow\; \sum_i x_i = n\mu \;\Longrightarrow\; \boxed{\hat\mu = \bar x}.\]The sample mean. Note this is independent of \(\sigma^2\), which is why we can solve for it first. It is also Lecture 13's result that the mean minimises squared deviation — the same computation.
For \(\sigma^2\): treat \(v=\sigma^2\) as the variable (cleaner than \(\sigma\)).
\[\frac{\partial\ell}{\partial v} = -\frac{n}{2v} + \frac{1}{2v^2}\sum_i (x_i-\mu)^2 = 0.\]Multiply by \(2v^2>0\): \(\;-nv + \sum_i(x_i-\mu)^2 = 0\), so
\[\boxed{\hat\sigma^2 = \frac{1}{n}\sum_{i=1}^n (x_i-\hat\mu)^2}.\]The catch. Note the \(n\), not \(n-1\). By Lecture 13's exercise, \(\E[\hat\sigma^2_{\text{MLE}}]=\frac{n-1}{n}\sigma^2\) — the maximum likelihood estimate of the variance is biased low, systematically underestimating the spread. It is consistent (the bias vanishes as \(n\to\infty\)) but not unbiased.
This is worth sitting with: maximum likelihood is a principle, not a guarantee. It is asymptotically efficient under regularity conditions, but at finite \(n\) it can be biased, and here it is biased in the dangerous direction — towards overconfidence. The same phenomenon is why a model's training loss underestimates its test loss.
Now the payoff. Two standard losses, both obtained by writing down a noise model and taking the negative log-likelihood.
MLE treats \(\vt\) as a fixed unknown. The Bayesian alternative treats it as random, with a prior, and maximises the posterior. Bayes' theorem gives
\[p(\vt\mid\mathcal{D}) = \frac{p(\mathcal{D}\mid\vt)\,p(\vt)}{p(\mathcal{D})} \;\propto\; \underbrace{p(\mathcal{D}\mid\vt)}_{\text{likelihood}}\ \underbrace{p(\vt)}_{\text{prior}},\]
where the evidence \(p(\mathcal{D})\) is dropped because it does not depend on \(\vt\) and so cannot move the argmax.
Three consequences worth carrying into Phase B.
That is Phase A. It is worth pausing to see the shape of it, because the pieces connect more tightly than they appeared to at the time.
Phase B builds classical machine learning on this. You will find that most of it is already half-derived: linear regression is Lecture 10, logistic regression is Lectures 2, 5 and 14, PCA is Lecture 11, gradient descent is Lectures 4, 6 and 12, and the bias-variance decomposition is Lecture 13. That is the intended effect. Foundations built properly make the material above them feel inevitable rather than arbitrary.
A factory has three machines producing 50%, 30% and 20% of output, with defect rates 1%, 2% and 3% respectively.
1. Law of total probability.
\[P(D) = 0.5(0.01)+0.3(0.02)+0.2(0.03) = 0.005+0.006+0.006 = 0.017 = 1.7\%.\]2. Bayes.
\[P(M_3\mid D) = \frac{P(D\mid M_3)P(M_3)}{P(D)} = \frac{0.03\times0.20}{0.017} = \frac{0.006}{0.017} \approx 0.3529.\]About 35.3%.
3. The three posteriors are \(\frac{0.005}{0.017}=29.4\%\), \(\frac{0.006}{0.017}=35.3\%\), \(\frac{0.006}{0.017}=35.3\%\). Machines 2 and 3 are tied as the most likely source, at 35.3% each.
And no, this is not the same as the highest defect rate: machine 3 is three times as defect-prone as machine 1, yet contributes only slightly more defects, because machine 1 produces two and a half times as much. Posterior = likelihood × prior, and the prior here is production volume. This is the base-rate lesson again, in an industrial setting: to reduce total defects you should weigh volume alongside rate.
Derive the MLE for the rate \(\lambda\) of a Poisson distribution, \(P(X=k)=\frac{\lambda^k e^{-\lambda}}{k!}\), from data \(k_1,\dots,k_n\).
Then state the corresponding loss function for a regression model that predicts counts, and say when you would use it in preference to squared error.
Log-likelihood.
\[\ell(\lambda) = \sum_{i=1}^n\Big[k_i\log\lambda - \lambda - \log(k_i!)\Big] = \left(\sum_i k_i\right)\log\lambda - n\lambda - \sum_i\log(k_i!).\]The last term does not involve \(\lambda\) and can be ignored for optimisation.
\[\ell'(\lambda) = \frac{\sum_i k_i}{\lambda} - n = 0 \;\Longrightarrow\; \boxed{\hat\lambda = \frac1n\sum_i k_i = \bar k}.\]The sample mean, again. And \(\ell''(\lambda) = -\frac{\sum k_i}{\lambda^2}<0\), so it is the unique maximum. ∎
Poisson regression. Model \(y_i\sim\text{Poisson}(\lambda_i)\) with \(\lambda_i=\exp(\vw^\top\vx_i)\) — the exponential link, which guarantees \(\lambda_i>0\) as a rate must be (a domain constraint, Lecture 1). The negative log-likelihood is
\[ L(\vw) = \sum_i\Big[\exp(\vw^\top\vx_i) - y_i\,\vw^\top\vx_i\Big], \]dropping the \(\log(y_i!)\) constant. Its gradient is \(\sum_i(\hat\lambda_i - y_i)\vx_i = \mX^\top(\hat{\boldsymbol\lambda}-\vy)\) — the same design-transposed-times-error shape as least squares and logistic regression, for the same exponential-family reason.
When to prefer it. Whenever the target is a non-negative integer count: website visits, defect counts, photon arrivals, word frequencies. Two concrete advantages over squared error:
Show that for the multivariate Gaussian \(\mathcal{N}(\vm,\mS)\), maximising the log-likelihood over \(\vm\) with \(\mS\) known gives \(\hat{\vm}=\bar{\vx}\), the sample mean vector.
Then explain what the level sets of the density look like in terms of the eigenvalues of \(\mS\), and what the quantity \((\vx-\vm)^\top\mS^{-1}(\vx-\vm)\) is called and used for.
Log-likelihood. Dropping terms without \(\vm\):
\[\ell(\vm) = -\frac12\sum_{i=1}^n (\vx_i-\vm)^\top\mS^{-1}(\vx_i-\vm) + \text{const}.\]Differentiate using Lecture 12's rule \(\grad_{\vz}(\vz^\top\mA\vz)=2\mA\vz\) for symmetric \(\mA\) (and \(\mS^{-1}\) is symmetric because \(\mS\) is), together with the chain rule for the inner \(\vz=\vx_i-\vm\), whose Jacobian with respect to \(\vm\) is \(-\mI\):
\[\grad_{\vm}\ell = \sum_{i=1}^n \mS^{-1}(\vx_i-\vm) = \mS^{-1}\left(\sum_i\vx_i - n\vm\right) = \vzero.\]Since \(\mS^{-1}\) is invertible we may multiply by \(\mS\): \(\sum_i\vx_i = n\vm\), so \(\hat{\vm}=\bar{\vx}\). ∎
The Hessian is \(-n\mS^{-1}\), which is negative definite because \(\mS\succ0\) implies \(\mS^{-1}\succ0\) (Lecture 11: inverting a positive definite matrix inverts its eigenvalues, keeping them positive). So it is the unique global maximum.
Level sets. The density depends on \(\vx\) only through the quadratic form, so level sets of the density are level sets of that form: ellipsoids centred at \(\vm\). By Lecture 11, if \(\mS=\mQ\Lambda\mQ^\top\) then \(\mS^{-1}=\mQ\Lambda^{-1}\mQ^\top\), so the form is \(\sum_i z_i^2/\lambda_i\) in the eigenbasis. The ellipsoid axes point along the eigenvectors \(\vq_i\), with half-lengths proportional to \(\sqrt{\lambda_i}\) — long in directions of high variance. An isotropic \(\mS=\sigma^2\mI\) gives spheres; strongly correlated variables give long thin cigars.
The quantity. \(d_M(\vx)^2 = (\vx-\vm)^\top\mS^{-1}(\vx-\vm)\) is the squared Mahalanobis distance. It is Euclidean distance measured in units of the local standard deviation, in every direction at once — equivalently, ordinary distance after whitening by \(\mS^{-1/2}\). Uses you will meet: multivariate outlier detection (under a Gaussian, \(d_M^2\sim\chi^2_d\), giving a principled threshold), the discriminant in Gaussian classifiers and LDA, and the natural-gradient and trust-region methods that measure parameter changes in the metric induced by the Fisher information rather than in raw Euclidean distance.
Coding + analysis, the Phase A capstone. Build a Bayesian linear regression by hand and compare it with MLE and MAP.
import numpy as np
rng = np.random.default_rng(3)
n, deg, sigma, tau = 15, 12, 0.25, 1.0
xs = np.sort(rng.uniform(-1, 1, n))
f = lambda x: 1.0 - 2.0*x + 0.5*x**2 + 3.0*x**3
y = f(xs) + sigma*rng.normal(size=n)
P = lambda x: np.vander(np.atleast_1d(x), deg+1, increasing=True)
X = P(xs)
w_mle, *_ = np.linalg.lstsq(X, y, rcond=None)
grid = np.linspace(-1.15, 1.15, 400); G = P(grid)
print("||w_mle|| =", np.linalg.norm(w_mle).round(2))
print("max |fit| on the grid =", np.abs(G @ w_mle).max().round(2),
" vs max |data| =", np.abs(y).max().round(2))
for lam in [1e-8, 1e-4, 1e-2, 1e0, 1e2]:
w = np.linalg.solve(X.T@X + lam*np.eye(deg+1), X.T@y)
train = np.mean((X@w - y)**2)
true = np.mean((G@w - f(grid))**2) # error against the TRUE function
print(f"lambda={lam:8.0e} ||w||={np.linalg.norm(w):9.2f}"
f" train MSE={train:.5f} error vs truth={true:.5f}")
S_N = np.linalg.inv(X.T@X/sigma**2 + np.eye(deg+1)/tau**2)
m_N = S_N @ X.T @ y / sigma**2
L = np.linalg.cholesky(S_N + 1e-12*np.eye(deg+1))
samples = m_N + (L @ rng.normal(size=(deg+1, 30))).T
band = samples @ G.T # 30 x len(grid)
print("\\nposterior sd of the prediction, inside vs outside the data:")
i0 = np.argmin(np.abs(grid))
print(f" x= 0.00 sd={band[:, i0].std():.4f}")
print(f" x= 1.15 sd={band[:, -1].std():.4f} <- outside the data range")2. With this seed \(\|\vw_{\text{MLE}}\|\approx 3.5\times10^4\), and the fitted curve reaches a magnitude of about \(9\times10^4\) on the plotting grid — against data whose largest value is \(1.89\). It passes essentially through every data point while oscillating violently between them and exploding outside the data range. Twelve degrees of freedom fitted to fifteen noisy points is near-interpolation: the model has fitted the noise. This is the textbook picture of overfitting, and also the ill-conditioning of Lecture 10 — the Vandermonde matrix here has \(\operatorname{cond}(\mX)\approx 2.8\times10^6\).
3. Error against the truth is U-shaped in \(\lambda\) and spans an enormous range: with this seed it runs \(3.9\times10^{4}\) at \(\lambda=10^{-8}\), down through \(129\) at \(10^{-4}\) and \(1.52\) at \(10^{-2}\), to a minimum of about \(0.87\) at \(\lambda=1\), then back up to \(1.71\) at \(\lambda=100\). Meanwhile training MSE rises monotonically the whole way, from \(0.011\) to \(1.18\) — training error alone would have told you to pick the worst model on the list. Too small and you overfit; too large and you underfit towards the zero function. Choosing it without peeking means cross-validation: hold out folds, fit on the rest, and pick the \(\lambda\) minimising held-out error. Never tune on the test set, and never tune against the true function — you do not have it. (The Bayesian alternative is to maximise the marginal likelihood, or "evidence", which integrates the parameters out rather than holding data out.)
4. The posterior samples agree tightly with each other where there is data and fan out dramatically outside the data range and in any gaps. Concretely: the predictive standard deviation is about \(0.09\) at \(x=0\), in the middle of the data, and about \(9.0\) at \(x=1.15\), just beyond it — a hundredfold increase over a short distance.
That is the thing a point estimate cannot tell you. \(\vw_{\text{MAP}}\) gives one curve with no indication of confidence; the posterior says here is what I know, and here is where I am guessing. This distinction matters enormously in any setting where acting on a confident-but-wrong prediction is costly — medicine, control, active learning, Bayesian optimisation — and it is why calibration and uncertainty quantification are active research areas rather than solved problems.
It is worth appreciating what made this tractable: the posterior is Gaussian in closed form only because a Gaussian prior is conjugate to a Gaussian likelihood. For a neural network no such closed form exists, and everything in approximate Bayesian deep learning — Laplace approximation, variational inference, deep ensembles, MC dropout — is an attempt to recover a usable version of the fan you just plotted.
Phase A complete. Take the index page's advice seriously before moving on: if there are exercises here you could not do, tell me which ones. That is more useful to me than a claim that everything made sense, and it decides whether Phase B starts next or whether we revisit something first.