Nyimbi Odero

Series 01 · a numbered run · kept in public

X0.3.2

Chapter: Phase 0 — The Bridge · Module 0.3 — Probability and Estimation

Lesson 0.3.2 — Maximum Likelihood and Maximum A Posteriori Estimation: The Two Point Estimates and When to Use Which

Week 4 · Prerequisites: 0.3.1, 0.2.1, 0.1.1, 0.1.2.

Prior results used. This lesson uses 0.3.1’s Bayes’ rule, the Gaussian pdf, and the conjugate update (precision form and covariance form); 0.2.1’s SE(3) composition and inverse (for the pose-estimation example); and 0.1.2’s autograd (for the gradient-based optimization of the log-likelihood). Each is restated at the point of use.

Hands forward. This lesson hands forward one capability: the learner can formulate a point-estimation problem as a maximum-likelihood or maximum-a-posteriori optimization, derive the closed-form solution for the Gaussian case, verify the solution by hand, and diagnose the four canonical estimation bugs (the missing log, the biased variance, the flat prior where a Gaussian is intended, and the MAP-as-posterior confusion) from their symptoms — the capability that Lesson 0.3.3 (filtering vs. smoothing) and 0.3.4 (the Kalman filter) assume.

Running case. The same 6-DOF MuJoCo arm from 0.2.1–0.3.1, with the noisy joint-2 encoder now read N = 5 times for the same unknown angle q₂. The ML estimate is the sample mean; the MAP estimate is the precision-weighted mean with the prior; and the lesson’s central result is that ML is MAP with a flat prior, and MAP is the conjugate update’s mode (which, for the Gaussian, is the same as its mean).


§L.0 — Lesson map, objectives, and misconceptions

§L.0.1 — Lesson map

PartContent
AMaximum likelihood: definition, log-likelihood, iid case, Gaussian closed form
BMaximum a posteriori: definition, prior, Gaussian closed form, the ML–MAP relationship
CNumerical stability: the log-sum-exp trick, why a product of densities leaves float64’s range, why the log is the optimization variable
DBias and regularization: the MLE variance’s bias, the MAP as regularization, the bias–variance trade-off
EFailure taxonomy, lab synthesis, cheat sheet, exercises, exit test

§L.0.2 — Objectives

By the end of this lesson you will be able to do eight things:

  1. State the maximum-likelihood (ML) estimator and the maximum-a-posteriori (MAP) estimator, and the one term that distinguishes them (the prior).
  2. Derive the ML estimator of a Gaussian’s mean and variance, and the MAP estimator of the mean under a Gaussian prior, by setting the derivative to zero and solving the resulting equation.
  3. Compute by hand the ML and MAP estimates for a concrete dataset, and verify the ML estimate equals the MAP estimate when the prior is flat.
  4. State the log-likelihood and its role: the log is a monotone transformation (so the argmax is unchanged), and the sum of logs is numerically stable where the product of pdfs is not.
  5. Apply the log-sum-exp trick to compute a log-likelihood whose naive product would leave float64’s range, and state the mechanism (the max subtraction).
  6. State the MLE variance’s bias (E[σ̂²_ML] = (N−1)/N σ²), and state the unbiased estimator (1/(N−1) denominator).
  7. Distinguish MAP from the posterior (MAP is a point estimate, the mode; the posterior is the full distribution), and state when the two coincide (Gaussian/symmetric posterior).
  8. Diagnose four estimation bugs — missing log (numerical range loss), biased variance (wrong normalization), flat prior where Gaussian intended (under-regularized estimate), and MAP-as-posterior (overconfident decision) — from their symptoms.

§L.0.3 — Misconception inventory

Wrong model the learner probably holdsWhy it is attractiveWhere displacedThe correct model in one sentence
M1 — “ML and MAP are the same thing.”Both are point estimates of the same parameter, and for large data the MAP estimate is close to the ML estimate, so the distinction feels pedantic.§L.2 (Part B)ML maximizes the likelihood `p(D
M2 — “The MLE variance (1/N)Σ(y−ȳ)² is unbiased.”It is the textbook MLE, and “the MLE is asymptotically unbiased” is a true statement that sounds like “it is unbiased.”§L.4 (Part D)E[(1/N)Σ(y−ȳ)²] = (N−1)/N σ², which is less than σ² for every finite N; the unbiased estimator divides by N−1, and the difference matters for small N (e.g., N = 5 gives a 20% underestimate).
M3 — “Maximizing the likelihood is the same as maximizing the log-likelihood, so the log is optional.”The log is monotone, so the argmax is the same; the log is a convenience, not a necessity.§L.3 (Part C)The argmax is the same (the log is monotone), but the numerical computation is not — the product of N densities leaves float64’s range for large N (underflowing or overflowing, depending on the model’s mean log-density), while the sum of logs stays well-scaled; the log is required for a numerically stable implementation, not just a convenience.
M4 — “MAP is Bayesian.”MAP uses the posterior `p(θD)`, which is a Bayesian object; the estimator is derived from Bayes’ rule.§L.2 (Part B)
M5 — “A flat prior does not matter.”A uniform prior has constant density, so it drops out of the argmax; the estimate is “just the data’s.”§L.2 (Part B)A proper flat prior (a uniform density on a bounded support) is a legitimate prior that gives ML; an improper flat prior (constant density on an unbounded support) is not a probability distribution, and the “posterior” it produces may not be integrable — the estimate is well-defined but the “posterior” is not a probability, which breaks any downstream Bayesian computation.

Each misconception is displaced in the section named.

Everything is self-contained. Where a prerequisite lesson is used, the result is restated at the point of use.


§L.0.5 — Retrieval warm-up

This is the tenth lesson. The warm-up draws from 0.2.1, 0.2.2, and 0.1.1–0.1.2 — 0.3.1 is one step back and is excluded; 0.2.3 and 0.2.4 are relevant too, but the three below suffice for spacing.

Warm-up 1 (from 0.2.1, the inverse). For T = (R, t), write T^{-1} and state the frame its translation is in. Answer: T^{-1} = (Rᵀ, −Rᵀt), and the translation −Rᵀt is in the child frame. How this is used: §L.2’s pose-estimation example uses T^{-1} to compute the pose of the camera relative to the object, and the direction of the inverse matters for which frame’s parameters are being estimated.

Warm-up 2 (from 0.2.2, the exp map’s singularity). log(R) for R = I is ω = 0, but the axis is undefined; log(R) for R a half-turn is ω = πu with the axis ambiguous. What are the two singularities, and what is the difference? Answer: θ = 0 is a rank singularity (the left Jacobian’s determinant is zero, but the value is defined); θ = π is a value singularity (the log’s output is ambiguous, and the gradient is unbounded). How this is used: §L.2’s MAP on a rotation needs the log as its parameterization, and the π singularity is the reason a two-chart or a π-clip is used in the estimate’s optimization.

Warm-up 3 (from 0.1.2, autograd). A loss L(θ) = Σ log p(y_i|θ) + log p(θ) is a sum of differentiable functions. Does .backward() flow through all of them? Answer: Yes — each term is a graph of differentiable operations, and the sum’s backward accumulates .grad on θ from every term. How this is used: §L.2’s gradient-based MAP is the sum’s gradient, which is exactly what the log-posterior provides.

Warm-up 4 (from 0.1.1, view vs. copy and in-place). A learner accumulates a log-likelihood as log_L += term in a loop. Is log_L a tensor or a Python float, and does the in-place += matter? Answer: If log_L starts as a Python float and term is a float (from .item()), log_L is a float; if term is a tensor, log_L is a tensor and += is in-place on the tensor. The in-place += on a tensor that is a leaf with requires_grad=True is an error; on a non-leaf, it is a valid accumulation. The safe pattern is log_L = log_L + term (rebinding), not log_L += term. How this is used: §L.3’s numerical-stability discussion uses the sum of logs, and the accumulation’s type (float vs. tensor) determines whether the gradient flows.


Part A — Maximum Likelihood Estimation

§L.1 — Objects: the estimator, the likelihood, and the log-likelihood

Given data D = {y_1, y_2, …, y_N} drawn from a distribution p(y|θ) parameterized by θ, the likelihood is L(θ) = p(D|θ), considered as a function of θ with D fixed. The maximum-likelihood (ML) estimator is:

θ^ML=arg⁡max⁡θ p(D∣θ)=arg⁡max⁡θ log⁡p(D∣θ)\hat{\theta}_{\text{ML}} = \arg\max_\theta \, p(D \mid \theta) = \arg\max_\theta \, \log p(D \mid \theta)

The log-likelihood ℓ(θ) = log p(D|θ) is the log of the likelihood, and the argmax is unchanged because log is monotone. The log is the optimization variable because:

  • The product p(y_1|θ) · p(y_2|θ) · … · p(y_N|θ) is a product of N factors, and it leaves float64’s range for large N — downward when the mean per-sample log-density is negative, upward when it is positive (Part C).
  • The log converts the product to a sum, which is well-scaled and differentiable term-by-term.
  • The log’s derivative is simpler: d/dθ log p(y|θ) = (1/p(y|θ)) · d/dθ p(y|θ) (the score function), and for the Gaussian, the score is linear in y.

For iid data (independent and identically distributed), the likelihood factors:

p(D∣θ)=∏i=1Np(yi∣θ),ℓ(θ)=∑i=1Nlog⁡p(yi∣θ)p(D \mid \theta) = \prod_{i=1}^{N} p(y_i \mid \theta), \qquad \ell(\theta) = \sum_{i=1}^{N} \log p(y_i \mid \theta)

and the ML estimator is the θ that maximizes the sum of the per-sample log-likelihoods.

§L.1.1 — Result: for iid Gaussian data with known variance, the ML mean is the sample mean; the ML variance is the biased sample variance

Claim. For y_1, …, y_N ~ iid N(μ, σ²), the ML estimate of μ is μ̂_ML = (1/N) Σ y_i = ȳ (the sample mean), and the ML estimate of σ² is σ̂²_ML = (1/N) Σ (y_i − ȳ)² (the biased sample variance, dividing by N, not N−1).

Proof sketch. The log-likelihood is ℓ(μ, σ²) = −(N/2) log(2πσ²) − (1/(2σ²)) Σ (y_i − μ)². Setting ∂ℓ/∂μ = 0 gives (1/σ²) Σ (y_i − μ) = 0 ⇒ Σ y_i = N μ ⇒ μ̂ = ȳ. Setting ∂ℓ/∂σ² = 0 gives −N/(2σ²) + (1/(2σ⁴)) Σ (y_i − μ̂)² = 0 ⇒ σ̂² = (1/N) Σ (y_i − μ̂)². The second derivative is negative at the solution (the log-likelihood is concave), so the stationary point is the maximum.

Corollary. The ML mean is unbiased (E[ȳ] = μ), and its variance is σ²/N (the 1/√n reduction from 0.3.1’s sum closure). The ML variance is biased: E[σ̂²_ML] = (N−1)/N σ² < σ², because ȳ is estimated from the same data and “uses up” one degree of freedom. The unbiased estimator divides by N − 1.

Corollary. The ML estimate’s asymptotic properties (consistency, asymptotic normality, efficiency) hold as N → ∞, but the finite-N properties (bias, variance) are what matter for a robot with a small number of sensor readings, and the bias is non-zero at every finite N.

§L.1.2 — Worked-example ladder

WE-L.1.a — Fully worked: the ML estimate for five joint-2 readings, by hand

Input. The running case’s joint-2 encoder readings: y = [0.5, 0.7, 0.6, 0.8, 0.4], with N = 5, assumed y_i ~ N(q₂, σ²).

Trace (ML mean).

  1. Σ y_i = 0.5 + 0.7 + 0.6 + 0.8 + 0.4 = 3.0.
  2. ȳ = 3.0 / 5 = 0.6.
  3. μ̂_ML = 0.6.

Trace (ML variance, biased).

  1. Deviations: [0.5−0.6, 0.7−0.6, 0.6−0.6, 0.8−0.6, 0.4−0.6] = [-0.1, 0.1, 0, 0.2, -0.2].
  2. Squared: [0.01, 0.01, 0, 0.04, 0.04].
  3. Σ = 0.10.
  4. σ̂²_ML = 0.10 / 5 = 0.02.

Trace (ML variance, unbiased).

  1. σ̂²_unbiased = 0.10 / 4 = 0.025.

Output and interpretation. The ML mean is 0.6, the biased variance is 0.02 (std ≈ 0.141), and the unbiased variance is 0.025 (std ≈ 0.158). The biased one is 20% smaller than the unbiased one at N = 5, and the difference shrinks as 1/N: at N = 100, the bias is 1%; at N = 1000, 0.1%. The choice between them depends on the intent: the ML estimate is the one that maximizes the likelihood (by definition); the unbiased estimate is the one whose expectation is the true variance (a property the ML estimate does not have).

Self-explanation 1. Before reading on: in step 4, the deviations sum to 0 — (-0.1) + 0.1 + 0 + 0.2 + (-0.2) = 0. Is this a coincidence, or a property of the ML mean?

Answer. It is a property of the ML mean: ȳ is defined as (1/N)Σy_i, so Σ(y_i − ȳ) = Σy_i − Nȳ = Nȳ − Nȳ = 0 for any data. The zero sum is the first-order condition of the ML mean’s optimization (∂ℓ/∂μ = 0), and it is why the sample mean is the “center” of the data — it is the point that makes the deviations sum to zero. The property is not a coincidence; it is the definition of the mean as the minimizer of the sum of squared deviations.

Self-explanation 2. Before reading on: in step 7, the ML variance is 0.02, and in step 8, the unbiased is 0.025. Which one should the learner report?

Answer. It depends on the intent: the ML variance is the maximum of the likelihood and is the “right” answer if the question is “what variance best fits the data under the Gaussian model”; the unbiased variance is the “right” answer if the question is “what variance would we expect to see if we repeated the experiment many times” (its expectation is the true variance, which the ML variance’s is not). The reporting should state the normalization (1/N or 1/(N−1)) explicitly, because the two differ and the difference is 1/N relative. In the lab, the choice is the learner’s, and the bias is the observable.

WE-L.1.b — Completion problem

Input. Readings y = [0.52, 0.48, 0.61, 0.55, 0.44] from the same encoder, N = 5.

Trace (final 25% as prompts).

  1. Σ y_i = 0.52 + 0.48 + 0.61 + 0.55 + 0.44 = 2.60. ⟨your turn: verify⟩
  2. ȳ = 2.60 / 5 = 0.52. ⟨your turn: verify⟩
  3. Deviations: [0.00, -0.04, 0.09, 0.03, -0.08]. ⟨your turn: verify the deviations sum to zero⟩
  4. Σ (y_i − ȳ)² = 0 + 0.0016 + 0.0081 + 0.0009 + 0.0064 = 0.0170. ⟨your turn: verify⟩
  5. σ̂²_ML = 0.0170 / 5 = 0.0034; σ̂²_unbiased = 0.0170 / 4 = 0.00425. ⟨your turn⟩

Answers.

  1. Sum: 2.60. ✓
  2. ȳ = 0.52. ✓
  3. Deviations sum: 0.00 - 0.04 + 0.09 + 0.03 - 0.08 = 0.00. ✓
  4. Squared sum: 0.0170. ✓
  5. Biased: 0.0034 (std ≈ 0.058); unbiased: 0.00425 (std ≈ 0.065).

The lesson. The ML mean is a weighted center — it is the ȳ that makes the deviations’ sum zero, and the variance is the average of the squared deviations from that center. The two variance estimators differ by the N/(N−1) factor, which is 25% at N = 5, and the choice between them should be stated in any report.

WE-L.1.c — Faded: the running case with more readings

Setup. The encoder is read N = 20 times, and the readings have ȳ = 0.6 and Σ(y_i − ȳ)² = 0.38.

Question. What are the ML and unbiased variance estimates? What is the difference’s relative size?

You must produce the intermediate steps. The setup and the answer are given.

Answer. σ̂²_ML = 0.38 / 20 = 0.019; σ̂²_unbiased = 0.38 / 19 ≈ 0.020. The difference is 0.020 − 0.019 = 0.001, which is 1/20 = 5% of the unbiased value. At N = 20, the bias is smaller than at N = 5 (20% vs. 5%), and the ML estimate’s bias is negligible for large N.

The lesson. The ML variance’s bias is 1/N relative, so it shrinks with the sample size. For N ≥ 30, the bias is < 3.3%, which is often smaller than the measurement’s own noise — the bias is a small-sample problem, not a large-sample one, and the choice between the two estimators matters only when N is small.

WE-L.1.d — Contrastive near-miss: the biased variance where the unbiased is intended

Correct solution. A learner who wants the unbiased variance for a small-sample estimate uses σ̂² = (1/(N−1)) Σ (y_i − ȳ)².

Near-miss. The learner uses σ̂² = (1/N) Σ (y_i − ȳ)² (the ML estimate), believing it is unbiased because “the MLE is asymptotically unbiased.”

Why the near-miss produces a plausible but incorrect number. The ML estimate is (N−1)/N times the true variance, so it is systematically smaller than the truth. For N = 5, it is 20% smaller; for N = 10, 10% smaller; for N = 30, 3.3% smaller. The underestimation is silent (both estimators are positive and “look reasonable”), and it propagates to any downstream computation that uses the variance: a confidence interval is too narrow, a Kalman filter’s covariance is too small (overconfident), a policy’s exploration is too greedy (under-explores). The error is systematic, not random, and it shrinks with N — so a large-N test does not catch it, and a small-N deployment is where it matters.

The discrepancy, stated. The ML estimate’s expectation is (N−1)/N σ², and the unbiased estimate’s is σ². The two differ by the 1/N factor, which is the degree of freedom used by the mean’s estimation. The learner who uses the ML estimate for an unbiased purpose has confused “the MLE is asymptotically unbiased” (true, as N → ∞) with “the MLE is unbiased at every N” (false, as the (N−1)/N factor shows).

Self-explanation 3. Before reading on: why does the ML variance use N and not N−1? What does the N represent, and what does the N−1 represent?

Answer. The N in the ML variance is the number of data points, because the ML estimator is the maximizer of the likelihood, and the likelihood’s σ²-dependent term is −N/(2σ²) + (1/(2σ⁴)) Σ(y_i − μ̂)², whose derivative is zero at σ̂² = (1/N)Σ(y_i − μ̂)². The N−1 in the unbiased estimator is the number of degrees of freedom after estimating the mean — the deviations y_i − ȳ have only N−1 independent values (they sum to zero), so the “average” that corrects for the mean’s estimation divides by N−1, not N. The two are different concepts (data count vs. degrees of freedom), and the ML’s N is the data count while the unbiased’s N−1 is the degrees of freedom.

§L.1.3 — Comparative table: the two variance normalizations and the estimator they produce

EstimatorFormulaExpectationWhat it answersWhen it matters
ML varianceσ̂²_ML = (1/N) Σ (y_i − ȳ)²(N−1)/N σ²The variance that maximizes the likelihood of this datasetWhenever the question is “which parameters best fit these data under the model” — the likelihood’s own answer
Unbiased varianceσ̂² = (1/(N−1)) Σ (y_i − ȳ)²σ²The variance that would be recovered, on average, across repeated experimentsWhenever the estimate is fed forward — a covariance, an interval, an exploration rate — and the error must not be systematic
ML meanμ̂_ML = (1/N) Σ y_i = ȳμThe center that makes the deviations sum to zeroAlways; the mean is unbiased and its variance is σ²/N
Deviation sumΣ (y_i − ȳ)0 for any dataThe first-order condition of the mean’s optimizationAs the check that ȳ is the mean — the sum is zero by construction, not by luck

Design principle (P10.6). State the normalization whenever you report a variance. 0.02 and 0.025 are the same data summarized by two estimators that differ by N/(N−1); at N = 5 that is a 20% difference, which is larger than most effects a small-sample robot experiment is trying to detect. The number alone does not say which question it answers, and the reader who assumes the wrong one inherits the bias silently.

§L.1.4 — Concept check

Concept check (displaces M2). A learner says: “The MLE variance (1/N)Σ(y−ȳ)² is unbiased — the MLE is asymptotically unbiased, so it is unbiased.” Which fact has the learner failed to model, and what class of bug will that omission produce?

Answer. The learner has failed to model the difference between asymptotic and finite-N bias, and the degree of freedom that ȳ consumes. E[(1/N)Σ(y−ȳ)²] = (N−1)/N σ², which is less than σ² at every finite N — at N = 5 it is 20% less, at N = 30 it is 3.3% less, and only in the limit does the factor reach 1. The bug class is a systematic underestimate that no large-N test catches: both estimators are positive and plausible, the difference shrinks with N, so a validation run on a big dataset agrees to within noise. The error surfaces where N is small and the variance is used: a confidence interval is too narrow, a filter’s covariance is too small and it over-trusts its own state, a policy’s exploration is too greedy. It is systematic rather than random, so averaging more datasets does not remove it.

§L.1.5 — Consolidation checkpoint 1

(a) You should now be able to: (i) write the likelihood and the log-likelihood for iid data, and say why the argmax is the same for both; (ii) derive μ̂_ML and σ̂²_ML by setting the derivatives to zero, and say which normalization falls out of the likelihood; (iii) compute both variance estimates for a dataset by hand and state the relative difference 1/N.

(b) Transfer probe. The encoder is read N = 4 times with readings [0.50, 0.56, 0.44, 0.50]. Compute μ̂_ML, σ̂²_ML, and the unbiased estimate. Then a colleague reports the same dataset’s variance as 0.0018 — which estimator did they use, and by what percentage are they off? Answer sketch: Sum 2.00, so ȳ = 0.50. Deviations [0, 0.06, −0.06, 0], squared [0, 0.0036, 0.0036, 0], sum 0.0072. Biased: 0.0072/4 = 0.0018; unbiased: 0.0072/3 = 0.0024. The colleague used the ML normalization, which is (N−1)/N = 3/4 of the unbiased value — 25% low at N = 4. A learner who checks the deviations’ sum (it is zero) and the mean (correct) has verified the mean and learned nothing about the variance, which is the check that actually discriminates.

(c) Remediation rule. If you cannot do (a)(ii) or (b), return to §L.1.2 and redo WE-L.1.a and WE-L.1.d, then re-attempt (b) by computing both normalizations before comparing to the colleague’s number.


Part B — Maximum A Posteriori Estimation

§L.2 — Objects: the MAP estimator and the prior’s role

The maximum a posteriori (MAP) estimator is:

θ^MAP=arg⁡max⁡θ p(θ∣D)=arg⁡max⁡θ [log⁡p(D∣θ)+log⁡p(θ)]\hat{\theta}_{\text{MAP}} = \arg\max_\theta \, p(\theta \mid D) = \arg\max_\theta \, \left[\log p(D \mid \theta) + \log p(\theta)\right]

where p(D|θ) is the likelihood (the same as in ML), and p(θ) is the prior. The log-posterior is the sum of the log-likelihood and the log-prior, and the MAP estimator is the θ that maximizes this sum.

The relationship to ML:

  • ML is MAP with a flat prior — if p(θ) is constant (a uniform or flat prior), the log p(θ) term is constant, and the argmax is the ML argmax.
  • MAP with a Gaussian prior is a regularized ML — the log-prior is a quadratic −(1/(2σ₀²))(θ − μ₀)², which penalizes θ for moving away from μ₀, and the penalty is stronger for smaller σ₀² (a tighter prior).
  • As N → ∞, MAP → ML — the log-likelihood grows with N (each observation contributes O(1)), while the log-prior is fixed (independent of N), so for large N the likelihood dominates and the prior’s influence vanishes. The rate of the convergence is O(1/N) for the mean, and the MAP estimator’s bias relative to the ML estimator is O(1/N).

The closed form for the Gaussian case: with prior N(μ₀, σ₀²) and N iid observations y_i ~ N(μ, σ²) (σ² known), the MAP estimate of μ is:

μ^MAP=σ2μ0+Nσ02yˉσ2+Nσ02\hat{\mu}_{\text{MAP}} = \frac{\sigma^2 \mu_0 + N \sigma_0^2 \bar{y}}{\sigma^2 + N \sigma_0^2}

which is the precision-weighted average of the prior’s mean and the sample mean, with weights 1/σ₀² (for the prior) and N/σ² (for the data). This is the same expression as 0.3.1’s conjugate update’s mean, because the Gaussian posterior is symmetric and its mean equals its mode.

§L.2.1 — Result: MAP with a Gaussian prior is the precision-weighted average, and MAP equals ML in the limit of infinite data

Claim. For y_1, …, y_N ~ iid N(μ, σ²) and prior μ ~ N(μ₀, σ₀²), the MAP estimate is μ̂_MAP = (σ² μ₀ + N σ₀² ȳ)/(σ² + N σ₀²), and it satisfies:

  1. μ̂_MAP = ȳ when σ₀² → ∞ (a flat prior), and μ̂_MAP = μ₀ when N = 0.
  2. μ̂_MAP → ȳ as N → ∞ (the data dominates).
  3. μ̂_MAP is between μ₀ and ȳ for every finite N (it is a weighted average).

Proof sketch. The log-posterior is −(N/2) log(2πσ²) − (1/(2σ²)) Σ (y_i − μ)² − (1/2) log(2πσ₀²) − (1/(2σ₀²))(μ − μ₀)². Setting ∂/∂μ = 0 gives (1/σ²)(Nȳ − Nμ) − (1/σ₀²)(μ − μ₀) = 0, which rearranges to μ (N/σ² + 1/σ₀²) = Nȳ/σ² + μ₀/σ₀², giving the displayed formula. The three properties follow from the formula’s limits and the positivity of the weights.

Corollary. The MAP estimate is a shrinkage estimator: it shrinks the ML estimate toward the prior mean, with the shrinkage factor N σ₀² / (σ² + N σ₀²). For small N or tight prior (σ₀² small), the shrinkage is strong; for large N or loose prior, the shrinkage is weak. The shrinkage is a bias (the MAP estimate is biased toward μ₀), and the bias is deliberate — it trades bias for variance, which is the subject of Part D.

Corollary. The MAP estimate is not the posterior, and it is not the posterior mean in general. For a symmetric posterior (e.g., Gaussian), the mode (MAP), the mean, and the median are the same point; for a skewed posterior, they differ, and the MAP is the least representative of the three (it is the peak, not the center). The Bayesian answer is the full posterior, not the MAP; the MAP is a point summary that discards the uncertainty.

§L.2.2 — Worked-example ladder

WE-L.2.a — Fully worked: the MAP estimate for the joint-2 readings, by hand

Input. The WE-L.1.a readings y = [0.5, 0.7, 0.6, 0.8, 0.4], N = 5, ȳ = 0.6. Prior N(μ₀ = 0.5, σ₀² = 0.01), measurement variance σ² = 0.04.

Trace.

  1. σ² μ₀ = 0.04 · 0.5 = 0.02.
  2. N σ₀² ȳ = 5 · 0.01 · 0.6 = 0.03.
  3. Numerator: 0.02 + 0.03 = 0.05.
  4. Denominator: σ² + N σ₀² = 0.04 + 0.05 = 0.09.
  5. μ̂_MAP = 0.05 / 0.09 ≈ 0.5556.

Output and interpretation. The MAP estimate 0.5556 lies between the prior’s 0.5 and the ML’s 0.6, and it sits closer to the ML: 0.6 − 0.5556 = 0.0444 against 0.5556 − 0.5 = 0.0556. That is what the precisions predict. The data’s precision is N/σ² = 5/0.04 = 125, the prior’s is 1/σ₀² = 1/0.01 = 100, and the MAP is the precision-weighted average (125·0.6 + 100·0.5)/225 = 0.5556, so the heavier weight belongs to the data.

Self-explanation 4. Before reading on: in step 5, the MAP estimate 0.5556 is not the average of 0.5 and 0.6 (which would be 0.55). Why is it 0.5556 and not 0.55?

Answer. The average would be 0.55, which is the unweighted mean of the prior and the data. The MAP is the precision-weighted mean, and the weights are 1/σ₀² = 100 for the prior and N/σ² = 125 for the data. The weighted mean is (0.5·100 + 0.6·125)/(100 + 125) = (50 + 75)/225 = 125/225 ≈ 0.5556. The difference from 0.55 is the data’s excess precision (125 vs. 100), which pulls the estimate toward the data (0.6). The unweighted average would be correct only if the two precisions were equal.

WE-L.2.b — Completion problem

Input. Same readings, but the prior is tighter (σ₀² = 0.0025, μ₀ = 0.5) and the measurement variance is the same (σ² = 0.04).

Trace (final 25% as prompts).

  1. 1/σ₀² = 1/0.0025 = 400, N/σ² = 5/0.04 = 125. ⟨your turn: which weight is larger?⟩
  2. μ̂_MAP = (0.04·0.5 + 5·0.0025·0.6)/(0.04 + 5·0.0025). ⟨your turn: compute⟩
  3. Compare to WE-L.2.a’s 0.5556. ⟨your turn: is the new MAP closer to the prior or to the data?⟩

Answers.

  1. 400 (prior) vs. 125 (data) — the prior’s weight is larger, so the MAP is closer to the prior.
  2. Numerator: 0.02 + 0.0075 = 0.0275. Denominator: 0.04 + 0.0125 = 0.0525. μ̂_MAP = 0.0275 / 0.0525 ≈ 0.5238.
  3. 0.5238 is closer to the prior 0.5 (difference 0.0238) than to the data 0.6 (difference 0.0762) — the tighter prior dominates.

The lesson. The MAP estimate is a precision-weighted average, and the tighter prior (σ₀² = 0.0025 vs. 0.01) shifts the estimate toward the prior, even with the same data. The shrinkage is the deliberate effect of the prior, and it is stronger for tighter priors and weaker for looser ones.

WE-L.2.c — Faded: the large-N limit

Setup. The same prior (μ₀ = 0.5, σ₀² = 0.01), the same measurement variance (σ² = 0.04), but N = 500 readings with ȳ = 0.6.

Question. What is the MAP estimate? How close is it to the ML estimate (0.6)? What does this say about the prior’s influence as N → ∞?

You must produce the intermediate steps. The setup and the answer are given.

Answer. Numerator: 0.04·0.5 + 500·0.01·0.6 = 0.02 + 3 = 3.02. Denominator: 0.04 + 500·0.01 = 0.04 + 5 = 5.04. μ̂_MAP = 3.02 / 5.04 ≈ 0.5992. The ML is 0.6, and the MAP is 0.5992 — a difference of 0.0008, which is 0.13% of the ML. As N → ∞, the MAP → ML, and the prior’s influence vanishes as O(1/N).

The lesson. The prior’s influence is O(1/N), so it is large for small N (where it regularizes) and negligible for large N (where the data dominates). The MAP’s bias is the same O(1/N), so it shrinks with the data — the prior is a small-sample tool, and its cost is a large-sample bias that becomes negligible.

WE-L.2.d — Contrastive near-miss: the MAP where a flat prior was intended

Correct solution. A learner who wants the ML estimate uses a flat prior (p(θ) ∝ 1), and the estimate is ȳ = 0.6.

Near-miss. The learner uses a Gaussian prior (N(0.5, 0.01)) but reports the estimate as the ML estimate, believing the prior “doesn’t matter” because “the data is what counts.”

Why the near-miss produces a plausible but incorrect number. The MAP estimate 0.5556 is between the prior and the ML, and it is closer to the ML than to the prior (in WE-L.2.a), so it looks like the ML estimate with a small correction. The learner who reports it as “the ML estimate” has mislabeled a MAP estimate, and any downstream computation that uses the estimate (a policy trained on it, a controller tuned to it) is using a biased estimate that is pulled toward the prior — which is correct if the prior is real and incorrect if the prior is a guess. The error is silent (the estimate is a legal number), and the “the prior doesn’t matter” belief is the source.

The discrepancy, stated. The MAP estimate with a Gaussian prior is μ̂_MAP = (σ² μ₀ + N σ₀² ȳ)/(σ² + N σ₀²), which equals the ML estimate ȳ only when σ₀² → ∞ (a flat prior). For any finite σ₀², the MAP is a weighted average that differs from ȳ by N σ₀² (ȳ − μ₀)/(σ² + N σ₀²), which is non-zero unless ȳ = μ₀. The learner who “doesn’t believe in priors” is still using one (a flat prior), and the choice of prior changes the estimate.

Self-explanation 5. Before reading on: in the near-miss, if the prior is N(μ₀, σ₀²) and the data’s ȳ is exactly μ₀, does the MAP differ from the ML?

Answer. No — if ȳ = μ₀, then μ̂_MAP = (σ² μ₀ + N σ₀² μ₀)/(σ² + N σ₀²) = μ₀ (σ² + N σ₀²)/(σ² + N σ₀²) = μ₀ = ȳ. The MAP and ML agree at the prior’s mean, which is the fixed point of the shrinkage. The prior’s influence is zero at its own mean and grows as ȳ moves away from μ₀, with the growth’s rate set by the prior’s precision.

§L.2.3 — Comparative table: ML, MAP, and the posterior

ObjectDefinitionWhat it isWhat it keepsWhat it discards
ML estimate`argmax_θ p(Dθ)`A pointThe data’s best-fitting parameter
MAP estimate`argmax_θ [log p(Dθ) + log p(θ)]`A point — the posterior’s modeThe data’s fit and the prior’s preference
Posterior mean`E[θD]`A point — the posterior’s center of massThe full posterior’s first moment
Posterior`p(θD) ∝ p(Dθ)p(θ)`A distribution

Design principle (P10.7). Decide which of the four objects the downstream code needs before choosing an estimator. ML, MAP, and the posterior mean all return one number, and for a symmetric posterior all three coincide — which is why the choice looks free in the Gaussian case and is not free in general. The posterior is the only object that carries the uncertainty, and any decision that is non-linear in θ — a threshold, a saturation, a safety check — reads that uncertainty whether or not the code supplies it.

§L.2.4 — Concept check

Concept check (displaces M1). A learner says: “ML and MAP are the same thing — both are point estimates of θ, and with enough data MAP converges to ML anyway, so the distinction is pedantic.” Which term has the learner failed to model, and what class of bug will that omission produce?

Answer. The learner has failed to model the prior as a term in the objective, and the fact that its influence is O(1/N) rather than zero at finite N. ML is MAP with a flat prior; with any finite σ₀² the two differ by the shrinkage Nσ₀²(ȳ − μ₀)/(σ² + Nσ₀²), which is non-zero whenever ȳ ≠ μ₀. WE-L.2.a’s 0.5556 against 0.6 is a 7.4% gap at N = 5, and WE-L.2.b’s tighter prior widens it to 0.5238 — a 12.7% gap on the same data. The bug class is a mislabeled estimate: the learner reports a MAP value as “the ML estimate,” so a reader assumes no prior was used and treats the shrinkage as data. Nothing raises — the number is legal and close to the ML — and the error is invisible until the prior is wrong, at which point the estimate is wrong in the direction of the prior and the report has hidden why.

§L.2.5 — Consolidation checkpoint 2

(a) You should now be able to: (i) write the MAP objective as log-likelihood plus log-prior, and derive the Gaussian closed form by setting the derivative to zero; (ii) express the MAP estimate as a precision-weighted average and read off the two weights, 1/σ₀² and N/σ²; (iii) state the two limits — flat prior gives ML, N → 0 gives the prior mean — and the O(1/N) rate at which the prior’s influence vanishes; (iv) distinguish the MAP from the posterior, and say when the mode and the mean coincide.

(b) Transfer probe. A second sensor is added to the running case: N₂ = 20 readings with ȳ₂ = 0.75, σ₂² = 0.09, and the same prior N(0.5, 0.01). Treating the first sensor’s five readings and the second’s twenty as one dataset with a common mean, what is the MAP estimate — and what goes wrong if a learner simply averages 0.6 and 0.75? Answer sketch: The data’s precision is N₁/σ₁² + N₂/σ₂² = 5/0.04 + 20/0.09 = 125 + 222.2 = 347.2, the prior’s is 100, and the numerator is 125·0.6 + 222.2·0.75 + 100·0.5 = 75 + 166.7 + 50 = 291.7. So μ̂_MAP = 291.7/447.2 ≈ 0.652. The naive average of the two means gives 0.675 — it weights the five-reading sensor and the twenty-reading sensor equally, when the second carries 1.78× the precision and 4× the count. The precision-weighted form is the only one that respects how much each source knows, and the unweighted average is the same error WE-L.2.a’s self-explanation isolates: 0.55 against 0.5556.

(c) Remediation rule. If you cannot do (a)(ii) or (b), return to §L.2.2 and redo WE-L.2.a and WE-L.2.b, then re-attempt (b) by writing both precisions before forming any average.


Part C — Numerical Stability and the Log-Sum-Exp Trick

§L.3 — Objects: the numerical-range problem and its fix

The likelihood p(D|θ) = ∏ p(y_i|θ) is a product of N factors. When each factor is < 1, the product underflows to zero: N = 1000 with each factor 0.1 gives 0.1^1000 = 10^{-1000}, well below float64’s smallest subnormal (≈ 4.9 × 10^{-324}), so it rounds to 0.0. The log of 0.0 is -∞, and the gradient of -∞ is undefined — the optimization breaks.

But the factors are densities, and a density is not a probability. A Gaussian density at its mean is 1/√(2πσ²), which exceeds 1 whenever σ² < 1/(2π) ≈ 0.159. The lesson’s own fitted model has σ̂² = 0.02, so its per-sample density is 1.71 and the product overflows — 1.71^1000 ≈ 1.8 × 10^{233} — long before it could underflow. Both failure modes come from the same place: a product of N numbers whose log is not near zero leaves the representable range in one direction or the other, and which direction depends on the sign of the per-sample log-density, not on N alone. The sign is a property of the model, not a universal fact, so a lesson that says “products underflow” is describing half the failure.

The log-sum-exp trick is the standard fix for the log of a sum of exponentials:

log⁡∑iexp⁡(ai)=max⁡iai+log⁡∑iexp⁡(ai−max⁡iai)\log \sum_i \exp(a_i) = \max_i a_i + \log \sum_i \exp(a_i - \max_i a_i)

The subtraction of the max makes each exponent non-positive (a_i − max ≤ 0), so exp(a_i − max) ≤ 1, and the sum is between 1 and N — well-scaled and no underflow. The max is added back at the end, so the result is exact (up to floating-point rounding).

For the likelihood, the log-sum-exp is not directly needed (the log of a product is a sum of logs, which is already stable). The log-sum-exp is needed when the likelihood is a mixture or a marginalization — e.g., p(y|θ) = Σ_k π_k p(y|θ, k), where the sum is over latent components (the mixture weights). The log of that sum is the log-sum-exp of the per-component log-likelihoods.

§L.3.1 — Result: the log converts a product to a sum, and the log-sum-exp converts a sum of exponentials to a max-plus-log-sum, both for numerical stability

Claim. For N factors p_i > 0, the log-likelihood ℓ = Σ log p_i is numerically stable (each log p_i is a moderate number, and the sum is moderate), while the product ∏ p_i leaves float64’s range for large N — downward if ℓ < 0, upward if ℓ > 0. For a sum Σ exp(a_i), the log-sum-exp max a + log Σ exp(a_i − max a) is stable, while the naive log Σ exp(a_i) underflows or overflows depending on the a_i’s magnitudes.

Proof sketch. The product is exp(ℓ), so it is representable exactly when ℓ lies in [log(5 × 10^{-324}), log(1.8 × 10^{308})] ≈ [-744.4, 709.8]. Outside that window the result is 0.0 or inf, and the sign of ℓ decides which. For p_i ≈ 0.1, ℓ ≈ -2302 at N = 1000 and the product is 10^{-1000}, far below the floor. For a Gaussian density above 1, ℓ > 0 and the same N overflows instead. The log itself is always representable: ℓ is a sum of moderate terms, and the sum is a plain float64 whose range is ±1.8 × 10^{308}. For the sum: if a_i are all 1000, exp(a_i) = e^{1000} overflows to inf; if a_i are all -1000, exp(a_i) underflows to 0. The max subtraction makes the exponentials bounded (between 0 and 1), and the sum is between 1 and N.

Corollary. The threshold is not N ≳ 308/|log p| — that formula divides a base-10 exponent by a natural log. The correct condition is N |log p| > 745.13 for underflow and N log p > 709.78 for overflow, with log natural. For p = 0.1 the underflow threshold is N ≈ 745.13/2.303 ≈ 323.6, so N = 323 survives and N = 324 does not, and direct evaluation confirms 0.1^323 = 9.9 × 10^{-324} while 0.1^324 = 0.0.

Corollary. The log-likelihood is the optimization variable — not the likelihood — because the optimization’s gradient is well-defined in log-space and undefined (or NaN) in the product’s space. The gradient ∂ℓ/∂θ = Σ ∂log p_i/∂θ is a sum of moderate terms, and it is what autograd computes.

Corollary. The log-sum-exp is needed whenever the likelihood involves a sum of terms — a mixture model, a marginalization over a latent variable, a softmax — and its naive computation overflows or underflows. The trick is exact (not an approximation) and costs one extra max-and-subtract per term, which is negligible.

§L.3.2 — Worked-example ladder

WE-L.3.a — Fully worked: the underflow and its fix

Input. N = 1000 iid observations with each per-observation likelihood p(y_i|θ) = 0.1.

Trace (naive product).

  1. p(D|θ) = 0.1^{1000}.
  2. 0.1 = 10^{-1}, so 0.1^{1000} = 10^{-1000}.
  3. 10^{-1000} is below the smallest subnormal (≈ 4.9 × 10^{-324}), let alone the minimum normal (≈ 2.2 × 10^{-308}), so it underflows to 0.0.
  4. log(0.0) = -∞, and ∂(-∞)/∂θ is NaN — the optimization breaks.

Trace (log-likelihood).

  1. ℓ = Σ log(0.1) = 1000 · log(0.1) = 1000 · (-2.3026) = -2302.6.
  2. -2302.6 is a representable float64 (well within ±1.8 × 10^308), and ∂ℓ/∂θ is Σ ∂log(0.1)/∂θ, which is a finite number.
  3. The optimization works.

Output and interpretation. For this hypothetical — every per-sample likelihood 0.1 — the naive product underflows at N = 1000 while the log-likelihood is representable. The threshold is N ≈ 745.13/|log 0.1| ≈ 323.6 for p = 0.1 — verified: 0.1^323 = 9.9 × 10^{-324} is still representable as a subnormal, and 0.1^324 = 0.0. The log is required for any N above the threshold, which is the typical case in robot learning (datasets with thousands of examples).

The direction of the failure is not fixed. A density above 1 makes ℓ > 0, and the product overflows instead of underflowing. The lesson’s own fitted Gaussian has σ̂² = 0.02 and a mean per-sample log-density of +0.537 (density 1.711), so at N = 1000 the product is 1.77 × 10^{233} — still representable, but climbing, and it overflows to inf at N ≈ 1322. A deliberately bad model (μ = 0, σ² = 4) has a mean log-density of −1.660 and underflows to 0.0 at N ≈ 449. Both are the same code path; the sign of ℓ picks the direction. That is why the lab’s step 4 measures both.

Self-explanation 6. Before reading on: in step 4, log(0.0) = -∞ and the gradient is NaN. Why is NaN the worse outcome than -∞?

Answer. -∞ is a value that signals the problem (the likelihood is zero), and a human can see the issue. NaN is a result of arithmetic (0/0, ∞ − ∞, ∞ · 0) and it propagates through every subsequent operation — once a gradient is NaN, the optimizer’s step is NaN, the parameter is NaN, and the entire training is broken with no way to recover except a restart. -∞ is a signal, NaN is a contagion; the log-likelihood’s finite value is the only way to avoid both.

WE-L.3.b — Completion problem

Input. A mixture model p(y|θ) = 0.5 · N(y; μ₁, 1) + 0.5 · N(y; μ₂, 1), with μ₁ = 0, μ₂ = 10, and y = 0.

Trace (final 25% as prompts).

  1. N(0; 0, 1) = 0.3989, N(0; 10, 1) = (1/√(2π)) exp(-50) ≈ 0.3989 · 1.9 × 10^{-22} ≈ 7.7 × 10^{-23}. ⟨your turn: verify⟩
  2. p(y|θ) = 0.5 · 0.3989 + 0.5 · 7.7 × 10^{-23} ≈ 0.1995 + 3.8 × 10^{-23} ≈ 0.1995. ⟨your turn: is the second term’s contribution visible in the sum?⟩
  3. log p(y|θ) = log(0.1995) ≈ -1.612. ⟨your turn: verify⟩
  4. The log-sum-exp version: log(0.5·exp(-0.9189) + 0.5·exp(-50 - 0.9189)). ⟨your turn: compute the two log-terms and the log-sum-exp, and compare to step 3⟩

Answers.

  1. N(0; 10, 1) = 0.3989 · exp(-50) ≈ 7.7 × 10^{-23}. ✓
  2. The second term is 3.8 × 10^{-23}, which is negligible relative to 0.1995 (the ratio is 2 × 10^{-22}), so the sum is 0.1995 to 22 decimal places — the second term’s contribution is invisible.
  3. log(0.1995) ≈ -1.612. ✓
  4. The two log-terms are log(0.5) + log(N(0; 0, 1)) = -0.6931 - 0.9189 = -1.612, and log(0.5) + log(N(0; 10, 1)) = -0.6931 - 50 - 0.9189 = -51.612. The log-sum-exp: max(-1.612, -51.612) = -1.612, and -1.612 + log(exp(0) + exp(-50)) = -1.612 + log(1 + 1.9 × 10^{-22}) ≈ -1.612 + 1.9 × 10^{-22} ≈ -1.612. ✓ Matches step 3.

The lesson. The mixture’s second component is negligible for y = 0 (it is 22 orders of magnitude smaller), and the naive sum handles it correctly — the underflow is not an issue here because the larger term is well-scaled. The log-sum-exp is needed when the largest term is itself tiny (e.g., all a_i are -1000), which is the case when the observation is far from all the components — a rare but important case (it is the “outlier” case, and the log-sum-exp handles it gracefully where the naive sum would underflow to 0).

WE-L.3.c — Faded: the log-sum-exp with all-negative log-terms

Setup. A mixture with three components, log-terms [-750, -752, -751] (all tiny).

Question. What is the naive log-sum-exp, and what is the stable version? What is the difference?

You must produce the intermediate steps. The setup and the answer are given.

Answer. Naive: log(exp(-750) + exp(-752) + exp(-751)). Each exponential underflows to 0.0 — exp(-750) = 0.0 in float64, because rounding to zero happens below log(2^{-1075}) ≈ −745.133 — so the sum is 0.0 and the log raises a domain error (or returns -inf in a vectorized library). Stable: max = -750, and log(exp(0) + exp(-2) + exp(-1)) = log(1 + 0.1353 + 0.3679) = log(1.5032) ≈ 0.4076, plus the max -750 gives -749.5924. The stable result is finite; the naive one is not.

The threshold matters here, and the first version of this example missed it. With [-500, -502, -501] the naive computation is fine: exp(-500) = 7.12 × 10^{-218} is comfortably representable, the naive sum is 1.071 × 10^{-217}, and its log is -499.5924 — identical to the stable value. The naive method fails only when the largest term is itself below the underflow line, max a_i < log(2^{-1075}) ≈ −745.13. Choosing a “tiny” example without checking that number produces a lesson that demonstrates nothing.

The lesson. The log-sum-exp’s max subtraction makes the sum finite when every term is below the representable floor — not merely when the terms are small. The condition is max a_i < −745.13, and a mixture whose largest term is above that line needs no trick at all. The stable version is exact (it adds the max back) and costs one subtract and one add per term.

WE-L.3.d — Contrastive near-miss: the log-likelihood computed as log(sum(pdf)) instead of sum(log(pdf))

Correct solution. ℓ = Σ log p(y_i|θ) — the sum of logs, term by term.

Near-miss. ℓ = log(Σ p(y_i|θ)) — the log of the sum, which is a different quantity (the log of the sum of the likelihoods, not the log-likelihood of the data).

Why the near-miss produces a plausible but incorrect number. The sum Σ p(y_i|θ) is a sum of probabilities, not the likelihood’s product; for N iid observations, the likelihood is ∏ p(y_i|θ), whose log is Σ log p(y_i|θ). The near-miss computes log(Σ p(y_i|θ)), which is log(N · average p) — a different function of θ — and its argmax is not the ML estimate. The near-miss’s value is positive (the sum of probabilities is > 0, and its log is finite), and it looks like a log-likelihood, but it is not. The error is silent (the function is differentiable, and its optimization runs), and the resulting θ is the wrong estimate.

The discrepancy, stated. Σ log p(y_i|θ) vs. log(Σ p(y_i|θ)) differ by the order of the log and the sum. The log of a sum is the log-sum-exp (Part C’s trick), and the log of a product is a sum of logs. The likelihood is a product (the iid assumption), so its log is a sum of logs — the learner who writes log(Σ p) has computed the log of a different likelihood (the sum-likelihood, which is the likelihood of a mixture with equal weights, not the iid likelihood). The two are not the same, and the argmax is different.

Self-explanation 7. Before reading on: in the near-miss, if N = 1, does the near-miss equal the correct log-likelihood?

Answer. Yes — for N = 1, Σ p(y_i|θ) = p(y_1|θ) and Σ log p(y_i|θ) = log p(y_1|θ), so log(Σ p) = log(p) = Σ log p. The two coincide at N = 1 and diverge for N ≥ 2. The error is invisible in a single-observation test and visible in a multi-observation one — the same “test at N = 1 is insufficient” pattern as the q = 0 test in the FK lessons.

§L.3.3 — Comparative table: four ways to compute a likelihood, and what each one costs

FormExpressionRepresentable whenFailure signatureWhen it is the right form
Naive product`∏ p(y_iθ)`ℓ = Σ log p_i lies inside ≈ [−744.4, 709.8]0.0 if ℓ < 0, inf if ℓ > 0 — no exception, no warning
Sum of logs`Σ log p(y_iθ)`Always, while each term is moderateNone; ℓ is an ordinary float64
Log of a sum (near-miss)`log(Σ p(y_iθ))`AlwaysNone — but the argmax is wrong
Log-sum-expmax a + log Σ exp(a_i − max a)Always, including when every a_i < −745.13Naive version gives -inf; this one is finite and exactA mixture, a marginalization over a latent variable, a softmax

Design principle (P10.8). Measure the sign of the per-sample log-density before predicting which way the product fails. The product is exp(ℓ), so the direction of its failure is the sign of ℓ, and ℓ’s sign is a property of the fitted model: a Gaussian density above 1 gives ℓ > 0 and an overflow, a density below 1 gives ℓ < 0 and an underflow. A lesson or a test that assumes underflow has described half the failure mode, and the lab’s own fitted model goes the other way.

§L.3.4 — Concept check

Concept check (displaces M3). A learner says: “The log is monotone, so maximizing the log-likelihood gives the same answer as maximizing the likelihood — the log is a convenience, not a necessity, and I can skip it for small N.” Which claim has the learner failed to model, and what class of bug will that omission produce?

Answer. The learner has correctly modeled the argmax and failed to model the numerical claim. The argmax is indeed unchanged; the computation is not. The product is exp(ℓ), so it is representable only while ℓ stays inside ≈ [−744.4, 709.8], and ℓ grows linearly in N — so the window closes at a finite N regardless of how moderate each factor is. For p = 0.1 the threshold is N ≈ 323.6: 0.1^323 = 9.9 × 10^{-324} survives and 0.1^324 = 0.0. The bug class is a silent range loss that becomes a NaN gradient: the product underflows to zero, log(0.0) = -∞, and the derivative of -∞ is NaN, which then contaminates every parameter in the next optimizer step. The failure arrives at a N the learner chose for statistical reasons, not numerical ones, so “small N” is not a defense — it is a bound the learner has not computed.

§L.3.5 — Consolidation checkpoint 3

(a) You should now be able to: (i) say why the log converts the product to a sum, and why the sum is the quantity that stays in range; (ii) state the underflow and overflow thresholds in natural log terms — N|log p| > 745.13 and N log p > 709.78 — and compute the threshold for a given per-sample likelihood; (iii) write the log-sum-exp identity, state why the max subtraction bounds every exponential, and say why the result is exact rather than approximate; (iv) say which likelihoods need the log-sum-exp at all (mixtures, marginalizations, softmaxes) and which do not (plain iid products).

(b) Transfer probe. A mixture has three components with log-terms [-400, -402, -401], and the learner’s naive code returns -399.5924 without complaint. Is the code wrong? Now run it on [-800, -802, -801]: what does the naive version return, and what does the stable version return? Answer sketch: At [-400, -402, -401] the code is correct — exp(-400) = 1.9 × 10^{-174} is comfortably representable, the naive sum is 1.07 × 10^{-174}, and its log is -399.5924, identical to the stable value. The trick is unnecessary there, and a test at that input proves nothing. At [-800, -802, -801] every exponential underflows to 0.0, so the naive sum is 0.0 and the log is -inf; the stable version takes max = -800 and returns -800 + log(1 + exp(-2) + exp(-1)) = -800 + 0.4076 = -799.5924. The discriminating condition is max a_i < −745.13, and it is a fact about the largest term, not about how small the terms look.

(c) Remediation rule. If you cannot do (a)(iii) or (b), return to §L.3.2 and redo WE-L.3.c, then re-attempt (b) by computing exp of the largest term before deciding whether the naive path survives.


Part D — Bias, Regularization, and the MAP as a Statistical Tool

§L.4 — Objects: bias, variance, and the MAP’s role

The bias of an estimator θ̂ is E[θ̂] − θ* (the difference between the estimator’s expectation and the true value). The variance of an estimator is E[(θ̂ − E[θ̂])²] (the estimator’s spread). The mean squared error (MSE) decomposes as:

MSE(θ^)=Bias(θ^)2+Var(θ^)\text{MSE}(\hat{\theta}) = \text{Bias}(\hat{\theta})^2 + \text{Var}(\hat{\theta})

The ML estimator of the variance σ² is biased (Part A), and the ML estimator of the mean is unbiased. The MAP estimator is biased (toward the prior), and its variance is smaller than the ML’s (the shrinkage reduces the spread). The bias-variance trade-off: the MAP’s deliberate bias reduces the variance, and the MSE may be lower for the MAP than for the ML, especially for small N.

The MAP as regularization: the log-prior’s role is a penalty on θ, and the Gaussian prior’s quadratic penalty (1/(2σ₀²))(θ − μ₀)² is the L2 regularization (weight decay) with strength 1/σ₀². The Laplace prior’s |θ − μ₀|/σ₀ penalty is the L1 regularization (lasso), which promotes sparsity. The regularization’s strength is the prior’s precision, and the regularization’s form is the prior’s shape.

§L.4.1 — Result: the MAP’s bias reduces its variance, and the MSE trade-off favors the MAP for small N

Claim. For the Gaussian mean estimation, the ML estimator’s bias is 0 and its variance is σ²/N; the MAP estimator’s bias is −(N σ₀²/(σ² + N σ₀²)) (ȳ − μ₀) (relative to the true μ), and its variance is (N σ₀²/(σ² + N σ₀²))² (σ²/N) — smaller than the ML’s by the square of the shrinkage factor. The MAP’s MSE is bias² + variance, which is lower than the ML’s for μ close to μ₀ and higher for μ far from μ₀ — the trade-off is favorable when the prior is accurate and unfavorable when it is wrong.

Proof sketch. The ML estimator’s bias is zero by construction (the sample mean is unbiased), and its variance is σ²/N (from the sum of independent Gaussians). The MAP estimator is a weighted average of the unbiased ML and the fixed prior mean, and its bias is the weight times the difference μ₀ − μ (if the true μ differs from μ₀, the MAP is pulled toward μ₀, which is a bias). Its variance is the square of the ML weight times the ML’s variance, which is smaller than the ML’s by the shrinkage factor squared. The MSE is the sum, and the trade-off’s sign depends on (μ − μ₀)² vs. the variance reduction.

Corollary. The regularization interpretation: the Gaussian prior’s log is −(1/(2σ₀²))(θ − μ₀)², which is exactly the L2 penalty’s form (a quadratic penalty around μ₀). The penalty’s strength 1/σ₀² is the prior’s precision — a tighter prior is a stronger regularizer, and an infinitely loose prior (σ₀² → ∞) is no regularization (ML). The regularization’s interpretation is the MAP’s statistical content, and the statistical interpretation is the MAP’s regularization form.

§L.4.2 — Worked-example ladder

WE-L.4.a — Fully worked: the bias and variance of the ML and MAP estimates, by hand

Input. True μ = 0.6, σ² = 0.04, prior μ₀ = 0.5, σ₀² = 0.01, N = 5. Compute the ML’s bias and variance, and the MAP’s bias and variance.

Trace (ML).

  1. E[μ̂_ML] = μ = 0.6, so bias = 0.6 − 0.6 = 0.
  2. Var[μ̂_ML] = σ²/N = 0.04/5 = 0.008.
  3. MSE_ML = 0 + 0.008 = 0.008.

Trace (MAP).

  1. The shrinkage factor α = N σ₀²/(σ² + N σ₀²) = 5·0.01/(0.04 + 0.05) = 0.05/0.09 ≈ 0.5556. Note μ̂_MAP = α·μ₀ + (1−α)·μ̂_ML.
  2. E[μ̂_MAP] = α·μ₀ + (1−α)·μ = 0.5556·0.5 + 0.4444·0.6 = 0.2778 + 0.2667 = 0.5444. Bias = 0.5444 − 0.6 = −0.0556.
  3. Var[μ̂_MAP] = (1−α)² · Var[μ̂_ML] = 0.4444² · 0.008 ≈ 0.1975 · 0.008 ≈ 0.00158.
  4. MSE_MAP = 0.0556² + 0.00158 ≈ 0.00309 + 0.00158 ≈ 0.00467.

Output and interpretation. The ML’s MSE is 0.008 (all variance), and the MAP’s is 0.00467 (bias squared 0.00309 plus variance 0.00158) — smaller than the ML’s by 42%. The MAP is better by MSE because the prior μ₀ = 0.5 is close to the truth μ = 0.6 (the bias is small), and the variance reduction is substantial. The bias scales with the prior’s distance from the truth: at μ = 5 it is α(μ₀ − μ) = 0.5556·(0.5 − 5) = −2.5, so the bias squared is 6.25 against the ML’s variance of 0.008 — the MAP would be 800× worse, and the prior that helped at μ = 0.6 would be actively harmful.

Self-explanation 8. Before reading on: in step 5, the MAP’s bias is −0.0556 (negative), meaning the MAP underestimates the true 0.6. Why negative, and what does it tell you about the prior’s direction?

Answer. The bias is negative because the prior 0.5 is below the true 0.6, and the shrinkage pulls the estimate toward the prior — the MAP is 0.5444, which is below the truth 0.6, so the bias E[μ̂_MAP] − μ = 0.5444 − 0.6 = −0.0556 is negative. The sign of the bias is the sign of (μ₀ − μ), and its magnitude is the shrinkage factor times |μ₀ − μ|. The MAP’s bias is toward the prior, which is good if the prior is accurate and bad if it is wrong.

WE-L.4.b — Completion problem

Input. Same as above, but the true μ is 0.5 (equal to the prior’s μ₀). Compute the MAP’s bias and variance.

Trace (final 25% as prompts).

  1. E[μ̂_MAP] = α·μ₀ + (1−α)·μ = 0.5556·0.5 + 0.4444·0.5. ⟨your turn⟩
  2. Bias: ⟨your turn⟩
  3. Variance: (1−α)²·σ²/N = 0.4444²·0.008. ⟨your turn⟩
  4. Compare to the ML’s MSE: ⟨your turn⟩

Answers.

  1. E[μ̂_MAP] = 0.5556·0.5 + 0.4444·0.5 = 0.5 — the MAP’s mean equals the true μ when μ = μ₀.
  2. Bias = 0.5 − 0.5 = 0 — the MAP is unbiased when the prior is exactly correct.
  3. Variance: 0.1975·0.008 ≈ 0.00158 — smaller than the ML’s 0.008.
  4. MSE: 0 + 0.00158 = 0.00158, which is 5× smaller than the ML’s 0.008. The MAP is much better when the prior is exactly correct — the shrinkage is pure gain (bias 0, variance reduced).

The lesson. The MAP’s bias vanishes when the prior’s mean equals the truth, and its variance is always smaller than the ML’s (by the shrinkage squared). The MAP is never worse than the ML in variance, and its MSE is better when the prior is accurate and worse when it is wrong — the trade-off is the prior’s accuracy.

WE-L.4.c — Faded: the regularization interpretation

Setup. A linear regression y = Xθ + ε with ε ~ N(0, σ²I), prior θ ~ N(0, σ₀²I). The ML loss is ‖y − Xθ‖², and the MAP loss is ‖y − Xθ‖² + λ‖θ‖² with λ = σ²/σ₀².

Question. Derive the MAP loss from the log-posterior, state the relationship between λ and the prior’s precision, and interpret λ = 0 and λ = ∞.

You must produce the intermediate steps. The setup and the answer are given.

Answer. Log-likelihood: −(N/2) log(2πσ²) − (1/(2σ²)) ‖y − Xθ‖². Log-prior: −(d/2) log(2πσ₀²) − (1/(2σ₀²)) ‖θ‖². The log-posterior is the sum, and the MAP is argmin_θ [‖y − Xθ‖²/σ² + ‖θ‖²/σ₀²] = argmin_θ [‖y − Xθ‖² + (σ²/σ₀²)‖θ‖²]. So λ = σ²/σ₀² — the ratio of the measurement’s variance to the prior’s variance. λ = 0 means σ₀² → ∞ (a flat prior), which is ML (no regularization). λ = ∞ means σ₀² → 0 (an infinitely tight prior at 0), which forces θ = 0 (the prior dominates completely). The intermediate λ is a trade-off, and its value is the prior’s precision relative to the data’s.

The lesson. The L2 regularization’s strength λ is the prior’s precision 1/σ₀² times the measurement’s variance σ² — a ratio, not an absolute value. The MAP’s regularization is statistically interpretable (a Gaussian prior), and the choice of λ is the choice of the prior’s tightness.

WE-L.4.d — Contrastive near-miss: the MAP reported as the posterior mean

Correct solution. The full posterior p(μ|D) = N(μ; μ̂_MAP, σ̂²_post) with σ̂²_post = 1/(1/σ₀² + N/σ²) is the Bayesian answer, and the MAP is its mode (which equals its mean for the symmetric posterior).

Near-miss. The learner reports the MAP estimate as “the posterior,” with no uncertainty, and any downstream computation that uses the MAP treats it as certain.

Why the near-miss produces a plausible but incorrect number. The MAP estimate is a single number, and the posterior is a distribution; the uncertainty (σ̂²_post) is discarded by the MAP, and any downstream decision that uses the estimate is overconfident — the robot trusts the estimate more than it should, and the cost is a grasp that is too aggressive, a navigation that is too direct, an exploration that is too greedy. The error is silent (the estimate is correct, the uncertainty is missing), and the decision’s behavior is wrong in a way that is only visible against a Bayesian decision that uses the uncertainty.

The discrepancy, stated. The MAP’s value is the posterior’s mode, and the posterior’s spread σ̂²_post is not zero — it is 1/(1/σ₀² + N/σ²), which is smaller than either input’s variance but not zero. The learner who reports the MAP as “the posterior” has replaced a distribution with a point, and the point’s value is correct but the uncertainty is lost. The fix is to report the full posterior (mean and covariance) and to use the distribution in any decision that is sensitive to uncertainty.

Self-explanation 9. Before reading on: in the near-miss, if the downstream decision is linear in the estimate (e.g., a control law u = Kμ̂), does the missing uncertainty matter?

Answer. For a linear decision, the expected behavior depends only on the mean of the posterior, and the MAP’s value is the mean (for the symmetric Gaussian), so the expected control is correct. But the variance of the control is non-zero (the estimate’s uncertainty propagates to the control’s uncertainty), and any non-linear decision (a saturating control, a threshold, a safety check) is sensitive to the variance. The missing uncertainty matters for non-linear decisions and for risk-sensitive ones, and the linear case is the exception, not the rule.

§L.4.3 — Comparative table: the prior, the penalty, and the estimator it produces

PriorLog-prior termEquivalent penaltyPenalty strengthEffect on the estimate
Flat (p(θ) ∝ 1)ConstantNone0ML; no shrinkage, full variance σ²/N
Gaussian N(μ₀, σ₀²)−(1/(2σ₀²))(θ − μ₀)²L2 (weight decay), λ‖θ − μ₀‖² with λ = σ²/σ₀²1/σ₀², the prior’s precisionShrinks toward μ₀; bias α(μ₀ − μ), variance (1−α)²σ²/N
Laplace Laplace(μ₀, b)`−θ − μ₀/b`L1 (lasso)
Point prior (σ₀² → 0)0 at θ = μ₀, −∞ elsewhereInfinite penalty∞θ̂ = μ₀; the prior dominates completely, no data consulted

Design principle (P10.9). The prior’s precision is the regularizer’s strength, and the prior’s shape is the regularizer’s form. 1/σ₀² sets how much the estimate is pulled, and the exponent on |θ − μ₀| sets where it is pulled from — a quadratic penalty tapers smoothly, an absolute-value penalty has a kink at zero and therefore produces exact zeros. Reading a regularization coefficient without its prior is reading half a specification: the same λ means different things under different data scales, because λ = σ²/σ₀² is a ratio.

§L.4.4 — Concept check

Concept check (displaces M4). A learner says: “MAP is Bayesian — it uses p(θ|D), it comes from Bayes’ rule, so a MAP estimate is the Bayesian answer.” Which distinction has the learner failed to model, and what class of bug will that omission produce?

Answer. The learner has failed to model the difference between a point and a distribution. MAP is the mode of the posterior — one number — and the Bayesian answer is the whole posterior p(θ|D), whose spread is σ̂²_post = 1/(1/σ₀² + N/σ²), which is smaller than either input’s variance but not zero. For the symmetric Gaussian posterior the mode, the mean, and the median coincide, so the value is right and the omission looks harmless; for a skewed posterior they differ, and the mode is the least representative of the three. The bug class is overconfident action: the downstream decision receives a point and treats it as certain, so a grasp is too aggressive, a navigation is too direct, an exploration is too greedy. It is silent because the estimate is correct — only the uncertainty is missing — and it becomes visible only against a decision that reads the variance: a saturating control, a threshold, a safety check. For a decision that is linear in θ, the expected behavior is right and the missing spread does not matter; that is the exception, not the rule.

§L.4.5 — Consolidation checkpoint 4

(a) You should now be able to: (i) decompose the MSE into bias squared plus variance, and state which term the MAP trades against which; (ii) compute the shrinkage factor α, the MAP’s bias, and the MAP’s variance for a concrete problem, and compare its MSE to the ML’s; (iii) state the two regularization correspondences — Gaussian prior to L2, Laplace prior to L1 — and the λ = σ²/σ₀² relationship; (iv) state why the MAP is not the posterior, and name the decision class for which the missing spread matters.

(b) Transfer probe. The running case’s prior is N(0.5, 0.01) and the data have ȳ = 0.6, σ² = 0.04, N = 5, giving α ≈ 0.5556, MSE_MAP ≈ 0.00467 against the ML’s 0.008. A second engineer proposes tightening the prior to σ₀² = 0.001 on the strength of that result. Compute α and the MAP’s MSE under the assumption that the prior mean 0.5 is correct, then under the assumption that the truth is μ = 0.8. What does the comparison say about when tightening is safe? Answer sketch: With σ₀² = 0.001, α = 5·0.001/(0.04 + 0.005) = 0.005/0.045 ≈ 0.111. If μ = μ₀ = 0.5, the bias is zero and MSE = (1−α)²σ²/N = 0.79·0.008 ≈ 0.0063 — worse than the loose prior’s 0.00467, because tightening also cuts the weight on the data that would otherwise reduce the variance. If μ = 0.8, the bias is α(μ₀ − μ) = 0.111·(−0.3) = −0.0333, squared 0.00111, plus variance 0.0063, giving 0.0074 — still below the ML’s 0.008, but the loose prior at the same wrong truth gives 0.5556·(−0.3) = −0.1667, squared 0.0278, which is 3.5× worse than the ML. The lesson the probe isolates: tightening helps only when the prior mean is right, and a prior that is right at μ = 0.5 becomes actively harmful at μ = 0.8 — the trade-off’s sign is the prior’s accuracy, not its tightness.

(c) Remediation rule. If you cannot do (a)(i) or (b), return to §L.4.2 and redo WE-L.4.a and WE-L.4.d, then re-attempt (b) by writing the bias and variance terms separately before summing them.


Part E — Failure Taxonomy, Lab, Cheat Sheet, Exercises, Exit Test

§L.5 — Failure taxonomy

SymptomThree most likely causes, rankedDiscriminating testFix
The optimization produces NaN after a few steps(1) The likelihood is computed as a product and underflows; (2) the log is taken of a zero; (3) the gradient is ∞ - ∞Print the likelihood’s value; if it is 0 or inf, it is (1)Compute the sum of logs, not the log of the product
The ML estimate equals the prior’s mean exactly(1) The prior is mistakenly used as the data’s mean; (2) N = 0 (no data); (3) the likelihood’s gradient is zero (flat likelihood)Check the log-likelihood’s value at the estimate; if it is constant, it is (3)Verify the data, the prior, and the likelihood
The MAP equals the ML(1) The prior is flat; (2) the prior’s precision is zero (σ₀² = ∞); (3) the data’s N is so large that the prior is negligibleCheck 1/σ₀²; if it is 0, it is (1) or (2)Use a proper prior, or verify the flat prior is intended
The variance estimate is biased downward(1) ML variance (1/N); (2) small N; (3) the data is drawn from a distribution with less variance than expectedCheck the normalization; if 1/N, it is (1)Use 1/(N-1) for the unbiased estimator
The posterior’s variance is larger than the prior’s(1) Variances added (should be precisions); (2) the “sum of independent” formula used for the posterior; (3) the posterior is the wrong directionCheck the formula; if it is σ² + σ_y², it is (1)Use 1/σ² + 1/σ_y² for precisions
The MAP is reported with no uncertainty(1) MAP-as-posterior confusion; (2) the posterior’s covariance is not computed; (3) the decision is assumed linearCheck whether the downstream computation uses a variance; if not, it is (1)Report the full posterior (mean and covariance), not just the MAP
The estimate is “correct” but the robot over-acts(1) The estimate’s uncertainty is ignored; (2) the decision is non-linear and risk-sensitive; (3) the prior is overconfidentCheck the decision’s dependence on σ̂²_post; if it is ignored, it is (1)Use the full posterior in the decision, not just the MAP

§L.6 — The lab

§L.6.1 — Task

Implement mle_gaussian, map_gaussian, log_likelihood, and log_posterior; verify the ML = MAP with a flat prior, the MAP = precision-weighted average, and the bias of the ML variance; then break the pipeline three ways.

§L.6.2 — Conditions

Python 3.10+, PyTorch 2.2+, CPU. Seed 42, dtype = float64. Runtime under 2 minutes. Data: y = [0.5, 0.7, 0.6, 0.8, 0.4] (the running case’s joint-2 readings). Prior: μ₀ = 0.5, σ₀² = 0.01. Measurement variance σ² = 0.04.

§L.6.3 — Controls held fixed

The seed, dtype, the data, the prior, and the measurement variance are fixed across all steps.

§L.6.4 — Metrics table

MetricDefinitionHow measured
mu_mlML estimate of μ(1/N) Σ y_i
sigma2_mlML estimate of σ² (biased)(1/N) Σ (y_i − ȳ)²
sigma2_unbiasedUnbiased estimate of σ²(1/(N-1)) Σ (y_i − ȳ)²
mu_mapMAP estimate of μ(σ² μ₀ + N σ₀² ȳ)/(σ² + N σ₀²)
log_likelihoodΣ log N(y_i; μ̂_ML, σ̂²_ML)The sum of logs
product_likelihood∏ N(y_i; μ̂_ML, σ̂²_ML)The product (leaves float64’s range for large N; direction set by the sign of ℓ)
mu_map_flat_priorMAP with σ₀² → ∞Should equal mu_ml
bias_ml_varianceE[σ̂²_ML] − σ² (empirical, over many datasets)Average over M datasets

§L.6.5 — Pre-registered hypotheses

  • H1. mu_ml = 0.6, sigma2_ml = 0.02, sigma2_unbiased = 0.025 (from WE-L.1.a). The prediction is the triple of exact values.
  • H2. mu_map ≈ 0.5556 (from WE-L.2.a). The prediction is 0.556 ± 0.001.
  • H3. log_likelihood is a finite number while product_likelihood is not representable, and the direction depends on the model: on the fitted model (mean per-sample log-density +0.537, density 1.711 > 1) the product overflows to inf at N ≈ 1322; on a deliberately bad model (μ = 0, σ² = 4, mean log-density −1.660) it underflows to 0.0 at N ≈ 449. At N = 1000 the fitted product is still representable (1.77 × 10^{233}), so the lab measures the direction rather than asserting one. The prediction is the pair: log finite, product out of range.
  • H4. mu_map_flat_prior ≈ mu_ml (within 1e-10 when σ₀² = 1e10). The prediction is the approximate equality.
  • H5. bias_ml_variance is negative (the ML variance underestimates the truth) and its magnitude is ≈ σ²/N = 0.008 for N = 5. The prediction is -0.008 ± 0.005 (the average over M = 1000 datasets).

§L.6.6 — The break-it sequence

  1. The product’s range (silent). Compute the likelihood as a product of per-sample densities for N = 1000 (extend the data by repeating it). Predicted symptom: the product leaves float64’s range with no exception and no warning — 0.0 if the mean log-density is negative, inf if it is positive — while the log-likelihood stays finite. Actual: as predicted, and the direction had to be measured: the fitted model’s mean log-density is +0.537, so its product climbs to 1.77 × 10^{233} at N = 1000 and overflows at N ≈ 1322; a bad model (μ = 0, σ² = 4) has mean log-density −1.660 and underflows to 0.0 at N ≈ 449. A lab that only ever runs the fitted model would report “the product is fine”, which is the wrong lesson.
  2. The biased variance (silent). Use 1/N instead of 1/(N-1) for the variance, and compare to the unbiased. Predicted symptom: the estimate is 20% smaller at N = 5, and the downstream σ² is systematically too small. Actual: as predicted — silent (both are positive).
  3. The flat-prior MAP (silent). Set σ₀² = 1e10 (effectively a flat prior) and verify the MAP equals the ML. Predicted symptom: the two are equal to 1e-10, and the prior’s influence is zero — which is correct for a flat prior but is the near-miss if a Gaussian prior was intended. Actual: as predicted.

All three are silent.

§L.6.7 — Acceptance criteria

CriterionThresholdArtefact it is tested on
mu_ml, sigma2_ml, sigma2_unbiased match the predictionsWithin 1e-10The lab’s ML prints
mu_map matches the predictionWithin 1e-3The lab’s MAP print
log_likelihood finite, product_likelihood out of float64’s range at N = 1000The pair, with the measured directionThe lab’s comparison
mu_map_flat_prior ≈ mu_mlWithin 1e-10The lab’s flat-prior check
bias_ml_variance negative and ≈ -0.008Within 0.005The lab’s bias print
Three break-its predicted3/3 timestampedThe prediction file

§L.6.8 — What good looks like

A strong submission has five properties:

  1. The closed forms are derived, not just coded. The submission shows the derivative ∂ℓ/∂μ = 0 and its solution, and the MAP’s derivation, so the formulas’ provenance is visible.
  2. The ML and MAP are compared, not just computed. The submission prints both, states the difference (0.6 vs. 0.5556), and explains the shrinkage.
  3. The log-likelihood’s stability is demonstrated, not just asserted. The submission shows the product leaving float64’s range at N = 1000 — in whichever direction its own model’s mean log-density dictates — and the log-likelihood finite, with the numeric values.
  4. The bias is measured empirically, not just computed. The submission runs M = 1000 datasets and reports the average ML variance’s deviation from the true σ².
  5. The three break-its are predicted before running, with the symptom (underflow, 20% bias, flat-prior equality) stated.

§L.6.9 — Runnable lab code

# lab_0_3_2.py — MLE, MAP, the log-likelihood, and the bias of the ML variance
# Requires: Python 3.10+, PyTorch 2.2+

import torch, math
torch.manual_seed(42); dtype, device = torch.float64, torch.device("cpu")

# --- Step 1: the data and the closed-form estimators -----------------
y = torch.tensor([0.5, 0.7, 0.6, 0.8, 0.4], dtype=dtype, device=device)
N = y.shape[0]
ybar = y.mean()
deviations = y - ybar
sigma2_ml = (deviations ** 2).mean()           # 1/N
sigma2_unb = (deviations ** 2).sum() / (N - 1) # 1/(N-1)
print(f"mu_ml = {ybar.item():.6f}, sigma2_ml = {sigma2_ml.item():.6f}, sigma2_unb = {sigma2_unb.item():.6f}")

# --- Step 2: the MAP estimate ---------------------------------------
sigma2 = 0.04; mu0 = 0.5; sigma0_sq = 0.01
mu_map = (sigma2 * mu0 + N * sigma0_sq * ybar) / (sigma2 + N * sigma0_sq)
print(f"mu_map = {mu_map.item():.6f}")

# Precision-weighted check
prec_data = N / sigma2; prec_prior = 1 / sigma0_sq
mu_map_check = (prec_data * ybar + prec_prior * mu0) / (prec_data + prec_prior)
print(f"mu_map (precision form) = {mu_map_check.item():.6f}")

# --- Step 3: the log-likelihood and the product ---------------------
def log_pdf(y, mu, var):
    return -0.5 * math.log(2 * math.pi * var) - (y - mu) ** 2 / (2 * var)

log_lik = sum(log_pdf(yi.item(), ybar.item(), sigma2_ml.item()) for yi in y)
prod_lik = torch.tensor(1.0, dtype=dtype)
for yi in y:
    prod_lik = prod_lik * torch.exp(torch.tensor(log_pdf(yi.item(), ybar.item(), sigma2_ml.item()), dtype=dtype))
print(f"log_likelihood = {log_lik:.6f}, product_likelihood = {prod_lik.item():.6e}")

# --- Step 4: the product leaves float64's range ------------------
# A density is not a probability. The fitted sigma2_ml = 0.02 gives a density
# of 1/sqrt(2*pi*0.02) = 3.99 at the mean, and the data's MEAN log-density is
# +0.537 — POSITIVE. So the product of N such factors OVERFLOWS, not underflows.
# Measure the sign first, then let it choose the direction.
def product_of_likelihoods(data, mu, var):
    """Naive product of per-sample densities, with the range outcome reported."""
    total = 0.0
    prod = 1.0
    for yi in data:
        lp = log_pdf(yi, mu, var)
        total += lp
        try:
            prod *= math.exp(lp)
        except OverflowError:
            return total, float('inf'), 'overflow'
    if prod == 0.0:
        return total, 0.0, 'underflow'
    return total, prod, 'representable'

y_many = y.repeat(200)  # N = 1000
# The product's fate is set by the MEAN log-density, not by one sample's.
mean_log = sum(log_pdf(yi.item(), ybar.item(), sigma2_ml.item()) for yi in y_many) / len(y_many)
print(f"mean per-sample log-density (fitted model) = {mean_log:.4f} "
      f"-> density {math.exp(mean_log):.4f}")
ll_fit, prod_fit, outcome_fit = product_of_likelihoods(y_many.tolist(), ybar.item(), sigma2_ml.item())
print(f"fitted model, N=1000: log_lik = {ll_fit:.4f}, product = {outcome_fit}")

# The same code path on a model whose density is below 1 goes the other way.
ll_bad, prod_bad, outcome_bad = product_of_likelihoods(y_many.tolist(), 0.0, 4.0)
print(f"bad model (mu=0, var=4), N=1000: log_lik = {ll_bad:.4f}, product = {outcome_bad}")

# The threshold is a NATURAL-log bound on N*|mean log p|, not the base-10 308.
print(f"overflow threshold N = {709.78 / mean_log:.1f} for the fitted model; "
      f"underflow threshold N = {745.13 / abs(ll_bad / len(y_many)):.1f} for the bad one")

# --- Step 5: the flat-prior limit --------------------------------
sigma0_sq_flat = 1e10
mu_map_flat = (sigma2 * mu0 + N * sigma0_sq_flat * ybar) / (sigma2 + N * sigma0_sq_flat)
print(f"mu_map (flat prior) = {mu_map_flat.item():.10f}, mu_ml = {ybar.item():.10f}")

# --- Step 6: the ML variance's bias, empirically ------------------
M = 1000
true_sigma2 = 0.04
biased_vars = torch.zeros(M, dtype=dtype)
for m in range(M):
    sample = torch.randn(N, dtype=dtype) * math.sqrt(true_sigma2) + 0.5
    biased_vars[m] = ((sample - sample.mean()) ** 2).mean()
print(f"E[sigma2_ml] over {M} datasets: {biased_vars.mean().item():.6f}")
print(f"true sigma2 = {true_sigma2}, bias = {biased_vars.mean().item() - true_sigma2:.6f}")
print(f"predicted bias = {-true_sigma2/N:.6f}")

# --- Step 7: break-it #1 — the product's range -------------------
print(f"[break-it #1] product: fitted model -> {outcome_fit}, bad model -> {outcome_bad} "
      f"(log finite in both: {ll_fit:.1f}, {ll_bad:.1f})")

# --- Step 8: break-it #2 — the biased variance -------------------
print(f"[break-it #2] sigma2_ml = {sigma2_ml.item():.6f}, sigma2_unb = {sigma2_unb.item():.6f}, ratio = {sigma2_ml.item()/sigma2_unb.item():.6f}")
# Expected ratio: (N-1)/N = 4/5 = 0.8

# --- Step 9: break-it #3 — the flat-prior MAP --------------------
print(f"[break-it #3] mu_map (flat) = {mu_map_flat.item():.10f}, mu_ml = {ybar.item():.10f}, equal: {abs(mu_map_flat.item() - ybar.item()) < 1e-8}")

# --- Step 10: summary --------------------------------------------
print("\n--- lab 0.3.2 summary ---")
print(f"H1: mu_ml={ybar.item():.4f}, sigma2_ml={sigma2_ml.item():.4f}, sigma2_unb={sigma2_unb.item():.4f}")
print(f"H2: mu_map={mu_map.item():.4f}")
print(f"H3: log_lik finite at N=1000 ({ll_fit:.2f}); product: fitted={outcome_fit}, bad={outcome_bad}")
print(f"H4: mu_map(flat) = {mu_map_flat.item():.10f} ≈ mu_ml = {ybar.item():.10f}")
print(f"H5: empirical bias = {biased_vars.mean().item() - true_sigma2:.6f} (predicted {-true_sigma2/N:.6f})")

§L.7 — The cheat sheet

Term or resultStatement or definitionRole
Likelihood`p(Dθ)`, the probability of the data given the parameter
Log-likelihood`ℓ(θ) = log p(Dθ) = Σ log p(y_i
MLE`θ̂_ML = argmax_θ p(Dθ)`
MAP`θ̂_MAP = argmax_θ [log p(Dθ) + log p(θ)]`
ML = MAP with flat priorThe two coincide when p(θ) is constantIntroduced here; consumed by 0.3.4
ML = MAP in the limitAs N → ∞, the prior’s influence is O(1/N)Introduced here; consumed by 0.3.4
Gaussian MLE meanμ̂ = ȳ = (1/N)Σy_iIntroduced here; consumed by 0.3.4
Gaussian MLE varianceσ̂² = (1/N)Σ(y_i − ȳ)² (biased)Introduced here; consumed by 0.3.4
Unbiased varianceσ̂² = (1/(N-1))Σ(y_i − ȳ)²Introduced here; consumed by 0.3.4
Gaussian MAP meanμ̂_MAP = (σ²μ₀ + Nσ₀²ȳ)/(σ² + Nσ₀²)Introduced here; consumed by 0.3.4
Shrinkage factorα = Nσ₀²/(σ² + Nσ₀²); μ̂_MAP = α μ₀ + (1−α) μ̂_MLIntroduced here; consumed by 0.3.4
Bias–variance trade-offMSE = Bias² + Var; MAP trades bias for varianceIntroduced here; consumed by 0.3.4
L2 regularization as MAPGaussian prior’s log = −λ‖θ‖²; λ = σ²/σ₀²Introduced here; consumed by every regularized learning lesson
L1 regularization as MAPLaplace prior’s log = −λ‖θ‖₁Introduced here; consumed by 3.1.2 (sparse VLA)
Log-sum-exp tricklog Σ exp(a_i) = max a + log Σ exp(a_i − max a)Introduced here; consumed by every mixture and marginalization lesson
Underflow threshold`Nlog p
MAP ≠ posteriorMAP is a point (the mode); the posterior is a distributionIntroduced here; consumed by 0.3.4 and every Bayesian decision lesson
Design principle P10.1ML is MAP with a flat prior; the prior is never “absent,” only flatIntroduced here; consumed by 0.3.4
Design principle P10.2Use the log-likelihood, not the likelihood, for any N > 100Introduced here; consumed by every training loop
Design principle P10.3The ML variance is biased; use 1/(N-1) for the unbiased estimatorIntroduced here; consumed by 0.3.4
Design principle P10.4MAP’s regularization strength is the prior’s precision relative to the data’sIntroduced here; consumed by every regularized learning lesson
Design principle P10.5Report the full posterior, not just the MAP, when uncertainty mattersIntroduced here; consumed by 3.3.4 (risk-sensitive decisions)
Design principle P10.6State the normalization whenever you report a varianceIntroduced here; consumed by 0.3.4
Design principle P10.7Decide which object the downstream code needs — point or distribution — before choosing an estimatorIntroduced here; consumed by 3.3.4
Design principle P10.8Measure the sign of the per-sample log-density before predicting the product’s failure directionIntroduced here; consumed by every training loop with large N
Design principle P10.9The prior’s precision is the regularizer’s strength; the prior’s shape is its formIntroduced here; consumed by every regularized learning lesson

§L.8 — Exercises with solutions

E 0.3.2.1 — The MLE by hand (numerical). Readings y = [1.0, 1.2, 0.8, 1.1, 0.9]. Compute μ̂_ML and σ̂²_ML (biased) and σ̂²_unbiased.

S 0.3.2.1. Sum = 5.0, μ̂_ML = 1.0. Deviations [0, 0.2, -0.2, 0.1, -0.1], squared [0, 0.04, 0.04, 0.01, 0.01], sum 0.10. Biased: 0.10/5 = 0.02. Unbiased: 0.10/4 = 0.025. The lesson: the ML mean is the sample average, and the variance is the average squared deviation (or its N/(N-1) correction for the unbiased version).

E 0.3.2.2 — The MAP by hand (numerical). Same readings, prior N(0.9, 0.04), measurement variance σ² = 0.05.

S 0.3.2.2. Numerator: σ²·μ₀ + N·σ₀²·ȳ = 0.05·0.9 + 5·0.04·1.0 = 0.045 + 0.2 = 0.245. Denominator: 0.05 + 5·0.04 = 0.05 + 0.2 = 0.25. μ̂_MAP = 0.245/0.25 = 0.98. The lesson: the MAP 0.98 is between the prior 0.9 and the ML 1.0, and it is closer to the ML because the data’s precision (N/σ² = 100) is much larger than the prior’s (1/σ₀² = 25). The shrinkage factor α = 0.04·5/0.25 = 0.2 gives μ̂_MAP = 0.2·0.9 + 0.8·1.0 = 0.18 + 0.8 = 0.98. ✓

E 0.3.2.3 — ML = MAP with a flat prior (numerical). Same readings, prior N(0.9, 1e10).

S 0.3.2.3. σ₀² = 1e10, so Nσ₀² = 5e10, which dominates the denominator (σ² + 5e10 ≈ 5e10). μ̂_MAP = (0.05·0.9 + 5e10·1.0)/(0.05 + 5e10) ≈ 5e10/5e10 = 1.0 = ȳ. The flat prior gives the ML. The lesson: the flat prior’s MAP is the ML, and the prior’s precision must be 0 (equivalently, σ₀² = ∞) for the two to coincide exactly; a finite σ₀² gives a weighted average, which is the MAP.

E 0.3.2.4 — The bias (numerical). For N = 10, what is the ML variance’s relative bias? For N = 100?

S 0.3.2.4. The relative bias is 1/N. For N = 10, it is 10% (the ML variance is 0.9 times the truth); for N = 100, it is 1%. The lesson: the ML variance’s bias is 1/N relative, so it is negligible for N > 100 and substantial for N < 10. The choice between the biased and unbiased estimators matters for small samples, and the unbiased is the default in most statistical software because it is correct for every N.

E 0.3.2.5 — The log-likelihood (numerical). For N = 500 iid observations each with per-sample likelihood 0.05, what is the product, and what is the log-likelihood? Which one is representable in float64?

S 0.3.2.5. Product: 0.05^{500} = (5 × 10^{-2})^{500} = 5^{500} × 10^{-1000}. log₁₀(5^{500}) = 500·log₁₀(5) ≈ 500·0.699 = 349.5, so the product is 10^{349.5 - 1000} = 10^{-650.5}. That is below even the smallest subnormal (10^{-323.3}), let alone the minimum normal (2.2 × 10^{-308}), so it underflows to 0.0. Log-likelihood: 500·ln(0.05) ≈ 500·(-2.996) ≈ -1498, which is representable — the log is an ordinary float64, and -1498 sits well inside the ±1.8 × 10^{308} range. The lesson: the log-likelihood is always representable when the per-sample log is moderate, even when the product is astronomically small — the log is the correct variable for the optimization, and the product is the wrong one for N above the underflow threshold.

E 0.3.2.6 — The log-sum-exp (numerical, transfer probe). Compute log(exp(-1000) + exp(-1001)) naively and with the log-sum-exp trick.

S 0.3.2.6. Naive: exp(-1000) ≈ 0 (underflows to 0 in float64), so the sum is 0 and the log is -inf. Stable: max = -1000, log(exp(0) + exp(-1)) = log(1 + 0.368) = log(1.368) ≈ 0.313, plus max = -1000 gives -999.687. The lesson: the log-sum-exp trick is exact (it adds the max back), and it is necessary when all the terms are tiny; the naive sum underflows and the log is -inf.

E 0.3.2.7 — The MAP vs. posterior (diagnostic). A learner reports “the posterior is N(0.5, 0.01)” — but they computed the MAP and the prior’s variance, not the posterior’s. What is the correct posterior, and what is the error?

S 0.3.2.7. The correct posterior is N(μ̂_MAP, σ̂²_post) with σ̂²_post = 1/(1/σ₀² + N/σ²) = 1/(100 + 125) = 1/225 ≈ 0.00444. The learner’s “posterior” has the prior’s variance 0.01, not the posterior’s 0.00444 — the posterior is more certain than the prior, and the learner’s variance is too large. The error is the “MAP-as-posterior” confusion combined with the wrong variance: the MAP’s value is right, but the uncertainty is the prior’s, not the posterior’s. The lesson: the posterior’s variance is the harmonic combination of the prior’s and the data’s, and it is smaller than either; reporting the prior’s variance as the posterior’s is the overconfidence’s opposite (underconfidence), and it is the same mislabeling of a summary as a distribution.

E 0.3.2.8 — The regularization interpretation (constructive, transfer probe). A learner’s linear regression uses loss = ‖y − Xθ‖² + λ‖θ‖². What is the statistical interpretation of λ, and what is the prior that produces it?

S 0.3.2.8. The λ‖θ‖² penalty is the L2 regularization, and its statistical interpretation is a Gaussian prior θ ~ N(0, σ₀²I) with λ = σ²/σ₀². The prior’s precision is 1/σ₀² = λ/σ², and the prior’s variance is σ₀² = σ²/λ. A larger λ means a tighter prior (smaller σ₀²), which is stronger regularization; a smaller λ means a looser prior, which is weaker regularization. λ = 0 is a flat prior (ML), and λ = ∞ is a point prior at θ = 0 (the prior dominates). The lesson: the regularization’s strength is the prior’s precision relative to the data’s, and the λ’s units are (data variance)/(prior variance), which is why the same λ can mean different things for different data scales.


§L.9 — Exit test

You are ready for Lesson 0.3.3 when you can, without notes:

  1. State the ML and MAP estimators, and the one term that distinguishes them.
  2. Derive the ML estimate of a Gaussian’s mean and variance by setting the derivative to zero.
  3. Compute by hand μ̂_ML and σ̂²_ML for a given dataset, and the unbiased σ̂².
  4. Derive the MAP estimate of the mean under a Gaussian prior, and state the precision-weighted average’s form.
  5. State the relationship ML = MAP with a flat prior, and the O(1/N) convergence of MAP to ML.
  6. State the shrinkage factor α, and interpret it as the prior’s precision relative to the data’s.
  7. State the ML variance’s bias (1/N) relative, and the unbiased estimator’s normalization 1/(N−1).
  8. State the log-sum-exp trick and why it is exact (not an approximation).
  9. State why the log-likelihood is the optimization variable, and the underflow threshold for the product.
  10. State the MAP’s regularization interpretation (Gaussian prior = L2, Laplace prior = L1) and the λ = σ²/σ₀² relationship.
  11. State the “MAP is not the posterior” distinction and its consequence for a risk-sensitive decision.
  12. Diagnose the “product underflow” (log missing), the “biased variance” (wrong normalization), and the “flat prior where Gaussian intended” (under-regularized estimate) from their symptoms.

If any item is not answerable without notes, return to the section that introduced it: 1 to §L.1 and §L.2, 2 to §L.1.1, 3 to §L.1.2, 4 to §L.2.1, 5 to §L.2.1, 6 to §L.2.1, 7 to §L.4.1, 8 to §L.3.1, 9 to §L.3.1, 10 to §L.4.1, 11 to §L.4.4, 12 to §L.5.

§L.10 — Provenance and evidence tiers

What this lesson reproduces. Every hypothesis in §L.6.5 is executed by the lab in §L.6.9. These claims are tier A, status R.

ClaimTierStatusWhere it is reproduced
ML on the five-point dataset gives μ = 0.6, σ²_ML = 0.02, σ²_unbiased = 0.025ARH1, all three printed
The MAP mean with a Gaussian prior (μ₀ = 0.5, σ₀² = 0.01) is 0.5556ARH2, matching WE-L.2.a
The log-likelihood stays finite at N = 1000 (537.07) while the product is not representableARH3, both branches printed
A flat prior reproduces ML exactly (0.6000000000 both ways)ARH4
The ML variance estimator is biased low: −0.008472 measured against −0.008000 predicted over 1000 datasetsARH5, ensemble mean and prediction both printed
The ML variance underestimates by the factor (N−1)/N = 0.8ARBreak-it #2, printed as 0.020000 / 0.025000 = 0.800000

What this lesson does not reproduce. H3’s overflow direction is stated in the text with a predicted crossing at N ≈ 1322 and an underflow at N ≈ 449 on the bad model; the lab prints the per-sample log-densities (+0.537 and −1.660) and the representability verdict, but not the crossing counts themselves. Those two integers are tier B, status S — derived from the printed per-sample value, not printed by the run. The arithmetic is one line (log10(N) = 308/0.537), which is why they are stated with a number rather than a bound, but a reader checking the lab for them will not find them.

The bias’s magnitude is the lesson’s sharpest measured claim, and it is an ensemble. −0.008472 is a mean over 1000 datasets, and the prediction −σ²/N = −0.008 is the asymptotic value; the 6% gap between them is the O(1/N²) correction, not a defect. A single dataset’s deviation from 0.04 is dominated by that dataset’s own noise, exactly as the lesson says for the smoother in 0.3.3. Testing this claim on one sample would return a number in the tens of percent and prove nothing.

No external sources are cited in this lesson. ML and MAP are defined in the lesson and exercised by the lab; the bias identity E[σ²_ML] = σ²(N−1)/N is derived in §L.3 rather than cited. No paper is named, and none is implied.