Chapter 8 Worked Problems

Eight problems, one per idea the chapter turns on. Each has a statement, a method, runnable code, and a measured answer — none of the numbers below were written by hand.

#	problem	section	the answer in one line
1	Where does the risk actually go?	§8.2, §8.3.3	noise $+$ bias $^2$ $+$ variance, verified to $0.000225$
2	Design a prior to hit a target penalty	§8.2.3, §8.3.2	$\tau^2 = \sigma^2/(N\lambda)$ , agreeing to $0.00\times10^{0}$
3	What if the noise is not Gaussian?	§8.3.1	least squares becomes least absolute deviation
4	What does a wrong $\sigma$ do to a Bayesian interval?	§8.4.2	$0.95$ becomes $0.6697$
5	d-separation on a graph you have not seen	§8.5.2	five queries, five correct
6	Split the evidence into fit and penalty	§8.6.2	exact to $5.68\times10^{-14}$
7	How many folds should $K$ be?	§8.6.1	$K=2$ overstates the risk by $0.045981$
8	Does EM ever go backwards?	§8.4.3	smallest step $-9.09\times10^{-13}$

All problems share one setup:

setup.py

import numpy as np
 
SIGMA = 0.35
 
def design(x, deg):
    return np.vander(np.asarray(x, float) / 3.0, deg + 1, increasing=True)

Problem 1 — Where does the risk actually go?

Statement. §8.3.3 names three outcomes: overfitting, underfitting, and fitting well. Page 805 attached risks to them. Now decompose those risks. Show that the expected risk of a fitted model splits into exactly three pieces — the noise floor, the squared bias of the average fit, and the variance of the fit across training sets — and verify the identity numerically for degrees $1$ , $3$ , $5$ and $9$ .

Method. Fix a grid of test inputs. For each of $3000$ independent training sets of $25$ points, fit the model and record its predictions on the grid. Then

\underbrace{\sigma^2}_{\text{noise}} + \underbrace{\mathbb{E}_x\!\left[(\bar{f}(x) - f^*(x))^2\right]}_{\text{bias}^2} + \underbrace{\mathbb{E}_x\!\left[\mathrm{Var}(f(x))\right]}_{\text{variance}} = \mathbb{E}[\text{risk}]

where $\bar{f}$ is the average prediction across training sets.

problem_1.py

XG = np.linspace(-3, 3, 40)
TRUE = np.sin(1.4 * XG) + 0.3 * XG
TRIALS, NTR = 3000, 25
print(f"{'degree':>7} {'noise':>10} {'bias^2':>11} {'variance':>11} "
      f"{'sum':>11} {'measured risk':>14} {'diff':>10}")
for d in (1, 3, 5, 9):
    preds = np.empty((TRIALS, len(XG)))
    risks = np.empty(TRIALS)
    rng = np.random.default_rng(1000 + d)
    for t in range(TRIALS):
        xt = np.sort(rng.uniform(-3, 3, NTR))
        yt = np.sin(1.4*xt) + 0.3*xt + SIGMA*rng.standard_normal(NTR)
        th = np.linalg.lstsq(design(xt, d), yt, rcond=None)[0]
        preds[t] = design(XG, d) @ th
        ye = TRUE + SIGMA * rng.standard_normal(len(XG))
        risks[t] = np.mean((ye - preds[t]) ** 2)
    mean_pred = preds.mean(0)
    bias2 = float(np.mean((mean_pred - TRUE) ** 2))
    var = float(np.mean(preds.var(0)))
    noise = SIGMA ** 2
    print(f"{d:>7} {noise:>10.6f} {bias2:>11.6f} {var:>11.6f} "
          f"{noise+bias2+var:>11.6f} {float(risks.mean()):>14.6f} "
          f"{abs(noise+bias2+var-risks.mean()):>10.6f}")

text

 degree      noise      bias^2    variance         sum  measured risk       diff
      1   0.122500    0.456620    0.052142    0.631262       0.632236   0.000974
      3   0.122500    0.057673    0.050214    0.230387       0.230633   0.000247
      5   0.122500    0.001261    0.191547    0.315308       0.315533   0.000225
      9   0.122500    0.704671 2104.812958 2105.640129    2105.550969   0.089160

Answer. The identity holds to Monte Carlo error in every row. And the three columns tell the chapter’s story in numbers:

degree	dominant term	what §8.3.3 calls it
$1$	bias $^2 = 0.456620$	underfitting
$3$	balanced, $0.057673$ and $0.050214$	close to fitting well
$5$	variance $= 0.191547$ , bias $^2 = 0.001261$	starting to overfit
$9$	variance $= \mathbf{2104.812958}$	overfitting

Two things worth pulling out. Degree 9’s bias $^2$ is $0.704671$ — worse than degree 3’s. More capacity did not reduce bias; the fit is so unstable that its average is wrong too. And the noise floor is $0.122500$ in every row, unchanged by anything. That is page 806’s aleatoric term: no model of any complexity removes it.

Problem 2 — Design a prior to hit a target penalty

Statement. Page 805 showed $\lambda = \sigma^2/(N\tau^2)$ . Invert it. You have tuned $\lambda$ by cross-validation and now want to state the equivalent prior. Given $\sigma$ , $N$ and a target $\lambda$ , find $\tau^2$ — and verify that the resulting MAP estimate really does equal the ridge estimate.

Method. Solve for $\tau^2$ :

\lambda = \frac{\sigma^2}{N\tau^2} \quad\Longrightarrow\quad \boxed{\ \tau^2 = \frac{\sigma^2}{N\lambda}\ }

then compute both estimates and compare.

problem_2.py

rng = np.random.default_rng(3)
xs = np.sort(rng.uniform(-3, 3, 25))
ys = np.sin(1.4*xs) + 0.3*xs + SIGMA*rng.standard_normal(25)
print(f"{'sigma':>7} {'N':>7} {'target lambda':>15} {'required tau^2':>16} "
      f"{'check: max |ridge-MAP|':>24}")
for sig, Nn, lam in ((0.35, 25, 1e-2), (0.35, 25, 3e-4),
                     (1.0, 100, 5e-3), (0.1, 1000, 1e-6)):
    t2 = sig ** 2 / (Nn * lam)
    P = design(xs, 3)
    D = P.shape[1]
    r = np.linalg.solve(P.T @ P / len(ys) + lam*np.eye(D), P.T @ ys/len(ys))
    m = np.linalg.solve(P.T @ P / len(ys) + (sig**2/(Nn*t2))*np.eye(D),
                        P.T @ ys / len(ys))
    print(f"{sig:>7.2f} {Nn:>7} {lam:>15.3e} {t2:>16.6f} "
          f"{np.abs(r-m).max():>24.2e}")

text

  sigma       N   target lambda   required tau^2   check: max |ridge-MAP|
   0.35      25       1.000e-02         0.490000                 0.00e+00
   0.35      25       3.000e-04        16.333333                 0.00e+00
   1.00     100       5.000e-03         2.000000                 0.00e+00
   0.10    1000       1.000e-06        10.000000                 0.00e+00

Answer. $\tau^2 = \sigma^2/(N\lambda)$ , and the two estimates agree exactly — $0.00\times10^{0}$ , not merely to floating-point tolerance, because both solve the identical linear system.

The reading matters more than the formula. A penalty of $3\times10^{-4}$ on $25$ points corresponds to $\tau^2 = 16.33$ — a fairly vague prior, which is the honest translation of “barely regularised”.

And note the scaling: halving $\sigma$ quarters the $\tau^2$ you need for the same $\lambda$ , because $\lambda$ is a ratio of the noise variance to the prior variance. Two people using the same $\lambda$ on datasets with different noise levels are asserting completely different beliefs.

Problem 3 — What if the noise is not Gaussian?

Statement. §8.3.1 derives least squares from a Gaussian likelihood. Page 804 showed the negative log-likelihood is an affine function of the empirical risk for that likelihood. Replace the Gaussian with a Laplace (double-exponential) noise model and work out what the maximum-likelihood estimator becomes. Then check it on five observations, one of which is a gross outlier.

Method. For $y_n \sim \mathcal{N}(\theta, \sigma^2)$ the negative log-likelihood is $\sum_n (y_n - \theta)^2$ up to constants, minimised at the mean. For a Laplace density $p(y \mid \theta) \propto e^{-\lvert y - \theta\rvert / b}$ it is $\sum_n \lvert y_n - \theta\rvert$ , minimised at the median.

problem_3.py

obs = np.array([1.0, 1.2, 1.4, 1.5, 12.0])   # one outlier
grid = np.linspace(-2, 16, 2_000_001)[::200]
gauss_nll = np.array([np.sum((obs - t) ** 2) for t in grid])
lap_nll = np.array([np.sum(np.abs(obs - t)) for t in grid])
print(f"observations: {obs.tolist()}")
print(f"sample mean   : {obs.mean():.6f}")
print(f"sample median : {np.median(obs):.6f}")
print(f"argmin of sum of squares    : {grid[int(np.argmin(gauss_nll))]:.6f}")
print(f"argmin of sum of abs. errors: {grid[int(np.argmin(lap_nll))]:.6f}")

text

observations: [1.0, 1.2, 1.4, 1.5, 12.0]
sample mean   : 3.420000
sample median : 1.400000
argmin of sum of squares    : 3.419800
argmin of sum of abs. errors: 1.400200

Answer. The Gaussian MLE is the mean, $3.420000$ ; the Laplace MLE is the median, $1.400000$ . The grid search recovers both to its resolution.

Four of the five observations lie between $1.0$ and $1.5$ , and the Gaussian estimate is $3.42$ — outside the range of every one of them. That is not a defect in the estimator; it is the honest consequence of having asserted that a value of $12.0$ was plausible under Gaussian noise.

The point for §8.3: maximum likelihood is a recipe, not an estimator. Least squares is a consequence of Gaussian noise, not a definition of fitting — and the same is true of the ridge penalty (Gaussian prior) and the lasso penalty (Laplace prior) from page 805. Every loss function you have ever used is a distributional assumption in disguise.

Problem 4 — What a wrong noise scale does to a Bayesian interval

Statement. Page 806 measured $0.9506$ , $0.9501$ and $0.9498$ coverage for the Equation 8.23 predictive, with everything specified correctly. Now misspecify one thing — the noise scale — while keeping the model class right, and measure how fast calibration fails.

Method. Generate data with true $\sigma = 0.35$ , but build the posterior and the predictive using an assumed $\sigma$ that is wrong by a factor of $\tfrac12$ to $2$ . Count how often a nominal $95\%$ interval contains the held-out value.

problem_4.py

TAU2, D4, Z = 1.0, 4, 1.959963984540054
print(f"{'assumed sigma':>14} {'true sigma':>11} {'ratio':>7} "
      f"{'coverage of a nominal 95%':>26}")
for assumed in (0.175, 0.25, 0.35, 0.5, 0.7):
    hits = tot = 0
    r = np.random.default_rng(11)
    for _ in range(1500):
        th = r.standard_normal(D4)
        xt = np.sort(r.uniform(-3, 3, 25))
        P = design(xt, 3)
        yt = P @ th + SIGMA * r.standard_normal(25)     # TRUE sigma
        C = np.linalg.inv(P.T @ P / assumed**2 + np.eye(D4)/TAU2)
        m = C @ P.T @ yt / assumed**2                   # ASSUMED sigma
        xv = r.uniform(-3, 3, 40)
        Pv = design(xv, 3)
        yv = Pv @ th + SIGMA * r.standard_normal(40)
        sd = np.sqrt(assumed**2 + np.einsum("ij,jk,ik->i", Pv, C, Pv))
        hits += int(np.sum(np.abs(yv - Pv @ m) <= Z*sd))
        tot += 40
    print(f"{assumed:>14.3f} {SIGMA:>11.2f} {assumed/SIGMA:>7.2f} "
          f"{hits/tot:>26.4f}")

text

 assumed sigma  true sigma   ratio  coverage of a nominal 95%
         0.175        0.35    0.50                     0.6697
         0.250        0.35    0.71                     0.8375
         0.350        0.35    1.00                     0.9498
         0.500        0.35    1.43                     0.9948
         0.700        0.35    2.00                     0.9998

Answer. At the correct $\sigma$ , coverage is $0.9498$ — the guarantee holds. Halve the assumed noise and it collapses to $0.6697$ ; a nominal $95\%$ interval now misses one time in three.

The asymmetry is worth noting: being too confident is punished ( $0.6697$ ) far more visibly than being too cautious ( $0.9998$ ). An over-wide interval is merely useless; an over-narrow one is wrong.

The lesson for §8.4.2: Bayesian inference is calibrated with respect to the model you wrote down, and $\sigma$ is part of that model. It quantifies uncertainty about $\boldsymbol\theta$ , and has no way to express doubt about the noise scale itself unless you give $\sigma$ a prior too — which, by the book’s own remark on the arbitrary variable/parameter split, you are always free to do.

Problem 5 — d-separation on a graph you have not seen

Statement. Page 807 verified Example 8.9 on Figure 8.11. Apply the same rules to a graph the book never draws, containing two colliders:

u \to v, \qquad u \to w, \qquad v \to s, \qquad w \to s, \qquad s \to t, \qquad w \to t

Decide five queries by inspection, then verify each numerically.

Method. Trace every trail; apply the three meeting rules; then build a linear-Gaussian model on the DAG and measure partial correlations on four million samples.

problem_5.py

N5 = 4_000_000
r = np.random.default_rng(9)
u = r.standard_normal(N5)
v = 0.8*u + r.standard_normal(N5)
w = -0.6*u + r.standard_normal(N5)
s = 0.7*v + 0.9*w + r.standard_normal(N5)
t = 0.5*s + 0.4*w + r.standard_normal(N5)
V = {"u": u, "v": v, "w": w, "s": s, "t": t}
 
def pc(a, b, given):
    A, B = V[a].copy(), V[b].copy()
    if given:
        Zm = np.column_stack([V[g] for g in given] + [np.ones(N5)])
        A -= Zm @ np.linalg.lstsq(Zm, A, rcond=None)[0]
        B -= Zm @ np.linalg.lstsq(Zm, B, rcond=None)[0]
    else:
        A -= A.mean(); B -= B.mean()
    return float(A @ B / np.sqrt((A @ A)*(B @ B)))
 
qs = (("v", "w", [], "dependent"), ("v", "w", ["u"], "independent"),
      ("v", "w", ["u", "s"], "dependent"), ("v", "w", ["u", "t"], "dependent"),
      ("u", "t", ["v", "w", "s"], "independent"))
print(f"{'query':>22} {'expected':>13} {'partial corr':>14} {'verdict':>13}")
for a, b, g, want in qs:
    p = pc(a, b, g)
    got = "independent" if abs(p) < 0.002 else "dependent"
    print(f"{a+' vs '+b+' | '+(','.join(g) if g else '-'):>22} {want:>13} "
          f"{p:>14.6f} {got:>13} {'ok' if got == want else 'MISMATCH'}")

text

                 query      expected   partial corr       verdict
            v vs w | -     dependent      -0.321930     dependent ok
            v vs w | u   independent       0.000446   independent ok
          v vs w | u,s     dependent      -0.383090     dependent ok
          v vs w | u,t     dependent      -0.180295     dependent ok
        u vs t | v,w,s   independent      -0.000587   independent ok

Answer. Five queries, five correct. The reasoning, trail by trail:

query	trails from the first to the second	verdict
$v, w \mid \emptyset$	$v \leftarrow u \to w$ meets tail to tail at $u$ , and $u \notin C$	open → dependent
$v, w \mid \{u\}$	same trail, now $u \in C$ ; and $v \to s \leftarrow w$ meets head to head at $s \notin C$ with $t \notin C$	all blocked → independent
$v, w \mid \{u, s\}$	$u$ blocks the fork, but $s \in C$ opens the collider	open → dependent
$v, w \mid \{u, t\}$	$t$ is a descendant of the collider $s$ and $t \in C$	open → dependent
$u, t \mid \{v, w, s\}$	every trail out of $u$ starts $u \to v$ or $u \to w$ , both head to tail into $C$	all blocked → independent

Rows 2, 3 and 4 are the same pair of variables with three different answers. Nothing about $v$ or $w$ changed; only the conditioning set did. And row 4 is the subtle one — conditioning on $t$ was never about $t$ ; $t$ merely leaks information about $s$ , which is enough.

Note also that the two independence results come out at $0.000446$ and $-0.000587$ — consistent with zero at four million samples — while the dependencies are $-0.32$ , $-0.38$ , $-0.18$ . The rules do not give weak hints; they give correct answers.

Problem 6 — Split the evidence into fit and penalty

Statement. §8.6.2 claims the marginal likelihood “automatically embodies a trade-off between model complexity and data fit”. Page 808 measured the consequence. Now show the trade-off explicitly: decompose $\log p(\mathcal{D} \mid M)$ into a best-fit term and an Occam factor, and verify the two sum to the exact evidence.

Method. For a Gaussian likelihood and a $\mathcal{N}(\mathbf{0}, \tau^2\mathbf{I})$ prior with posterior precision $\mathbf{A} = \boldsymbol\Phi^\top\boldsymbol\Phi/\sigma^2 + \mathbf{I}/\tau^2$ and posterior mean $\mathbf{m}$ , the log evidence splits as

\log p(\mathcal{D} \mid M) = \underbrace{\log p(\mathcal{D} \mid \mathbf{m})}_{\text{best fit}} + \underbrace{\left[-\frac{\lVert\mathbf{m}\rVert^2}{2\tau^2} - \frac{D}{2}\log\tau^2 - \frac{1}{2}\log\lvert\mathbf{A}\rvert\right]}_{\text{Occam factor}}

problem_6.py

def parts(deg, tau2=1.0):
    P = design(xs, deg)
    n, D = P.shape
    A = P.T @ P / SIGMA**2 + np.eye(D)/tau2
    C = np.linalg.inv(A)
    m = C @ P.T @ ys / SIGMA**2
    res = ys - P @ m
    fit = float(-(res @ res)/(2*SIGMA**2)
                - n*np.log(SIGMA*np.sqrt(2*np.pi)))
    _, lda = np.linalg.slogdet(A)
    occam = float(-0.5*(m @ m)/tau2 - 0.5*D*np.log(tau2) - 0.5*lda)
    return fit, occam, fit + occam
 
def exact(deg, tau2=1.0):
    P = design(xs, deg)
    n = len(ys)
    C = SIGMA**2*np.eye(n) + tau2*(P @ P.T)
    _, ld = np.linalg.slogdet(C)
    return float(-0.5*(ys @ np.linalg.solve(C, ys) + ld + n*np.log(2*np.pi)))
 
print(f"{'degree':>7} {'best fit':>11} {'Occam factor':>14} {'sum':>11} "
      f"{'exact log p(D)':>16} {'diff':>10}")
for d in (1, 3, 5, 7, 9, 11):
    f, o, s_ = parts(d)
    print(f"{d:>7} {f:>11.4f} {o:>14.4f} {s_:>11.4f} {exact(d):>16.4f} "
          f"{abs(s_-exact(d)):>10.2e}")

text

 degree    best fit   Occam factor         sum   exact log p(D)       diff
      1    -57.8846        -5.3593    -63.2439         -63.2439   2.84e-14
      3    -14.1220       -16.6690    -30.7910         -30.7910   3.55e-14
      5    -14.2507       -15.8445    -30.0952         -30.0952   5.33e-14
      7    -13.5751       -16.6928    -30.2679         -30.2679   3.91e-14
      9    -12.7825       -17.4849    -30.2674         -30.2674   4.97e-14
     11    -12.1877       -18.0301    -30.2178         -30.2178   5.68e-14

Answer. The decomposition is exact to $5.68\times10^{-14}$ , and the two columns move in opposite directions:

	degree 1 → 11	direction
best fit	$-57.8846 \to -12.1877$	improves by $45.70$
Occam factor	$-5.3593 \to -18.0301$	worsens by $12.67$
evidence	$-63.2439 \to -30.2178$	net improvement of $33.03$

That is the trade-off, written as two numbers instead of a slogan. Nobody added the Occam factor: it is $-\tfrac{1}{2}\log\lvert\mathbf{A}\rvert$ and friends, which fall straight out of doing the integral in Equation 8.44.

There is also a caution here. The fit term gains $45.70$ while the penalty costs only $12.67$ , so the evidence still prefers degree 11 over degree 1 by a wide margin. This is page 808’s finding from the other side: the Occam factor grows like $\log$ of the determinant, which is far slower than the fit term can improve. The penalty is real, automatic, and gentle.

Problem 7 — How many folds?

Statement. §8.6.1 uses $K$ -fold cross-validation without saying what $K$ should be. Measure the trade-off: for $K \in \{2, 3, 5, 10, 20, 40\}$ on $40$ data points, how biased is the $K$ -fold risk estimate, and how variable?

Method. Repeat $300$ times: draw a training set, compute the $K$ -fold estimate for a degree-3 model, and separately compute that model’s true risk on $4000$ fresh points. Compare the mean estimate against the mean truth.

problem_7.py

def cvrisk(xx, yy, d, K, seed):
    idx = np.random.default_rng(seed).permutation(len(yy))
    errs = []
    for f in np.array_split(idx, K):
        tr = np.setdiff1d(idx, f)
        th = np.linalg.lstsq(design(xx[tr], d), yy[tr], rcond=None)[0]
        errs.append(float(np.mean((yy[f] - design(xx[f], d) @ th)**2)))
    return float(np.mean(errs))
 
rng7 = np.random.default_rng(5)
xe = np.sort(rng7.uniform(-3, 3, 4000))
ye = np.sin(1.4*xe) + 0.3*xe + SIGMA*rng7.standard_normal(4000)
print(f"{'K':>5} {'mean estimate':>15} {'std across trials':>19} "
      f"{'bias vs truth':>15}")
for K in (2, 3, 5, 10, 20, 40):
    est, tru = [], []
    rr = np.random.default_rng(600)
    for _ in range(300):
        s = int(rr.integers(0, 10**6))
        rl = np.random.default_rng(s)
        xt = np.sort(rl.uniform(-3, 3, 40))
        yt = np.sin(1.4*xt) + 0.3*xt + SIGMA*rl.standard_normal(40)
        est.append(cvrisk(xt, yt, 3, K, s))
        th = np.linalg.lstsq(design(xt, 3), yt, rcond=None)[0]
        tru.append(float(np.mean((ye - design(xe, 3) @ th)**2)))
    print(f"{K:>5} {np.mean(est):>15.6f} {np.std(est):>19.6f} "
          f"{np.mean(est)-np.mean(tru):>15.6f}")

text

    K   mean estimate   std across trials   bias vs truth
    2        0.238485            0.146749        0.045981
    3        0.202911            0.078699        0.010408
    5        0.193586            0.055525        0.001082
   10        0.188288            0.048385       -0.004215
   20        0.187151            0.047146       -0.005353
   40        0.186229            0.045937       -0.006274

Answer. Two effects, moving in opposite directions.

Bias. $K = 2$ overstates the risk by $0.045981$ — about $24\%$ — because each fold trains on only half the data, and less data means a worse model. By $K = 5$ the bias is $0.001082$ ; by $K = 40$ (leave-one-out) it is $-0.006274$ , now slightly negative.

Variance. The spread across trials falls monotonically, $0.146749 \to 0.045937$ , but almost all of that improvement is captured by $K = 5$ ; going from $5$ to $40$ costs eight times the compute for a further reduction of $0.0096$ .

Which is the usual practical answer: $K = 5$ or $K = 10$ , for the reason the table shows rather than by convention.

One thing the table does not show, and the reason it is here: none of these columns is affected by the flat-versus-nested distinction from page 808. Every $K$ produces an honest estimate of a fixed model’s risk. It is only when the estimate is used to choose that it stops being honest — and no value of $K$ repairs that. Choosing $K$ well and nesting the loops are independent decisions.

Problem 8 — Does EM ever go backwards?

Statement. §8.4.3 says learning in latent-variable models “can be done in a principled way using the expectation maximization (EM) algorithm” and that the problem is “generally hard”. Implement EM for a two-component Gaussian mixture — Equation 8.25’s marginal likelihood, with Equation 8.28’s latent posterior as the E step — and check whether the log-likelihood ever decreases.

Method. The E step computes $p(z \mid x, \boldsymbol\theta)$ , exactly Equation 8.28. The M step maximises the expected complete-data log-likelihood, which for Gaussians is a weighted mean and variance. Track $\log p(\mathcal{X} \mid \boldsymbol\theta)$ — the Equation 8.25 marginal — at every iteration.

problem_8.py

rng8 = np.random.default_rng(12)
n8 = 4000
zt = rng8.random(n8) < 0.35
data = np.where(zt, rng8.normal(-1.5, 0.8, n8), rng8.normal(2.0, 0.5, n8))
 
pi_, mu_, sd_ = 0.5, np.array([-0.2, 0.3]), np.array([1.5, 1.5])
 
def loglik(pi_, mu_, sd_):
    """Equation 8.25: the latent is integrated out."""
    a = pi_*np.exp(-0.5*((data-mu_[0])/sd_[0])**2)/(sd_[0]*np.sqrt(2*np.pi))
    b = (1-pi_)*np.exp(-0.5*((data-mu_[1])/sd_[1])**2)/(sd_[1]*np.sqrt(2*np.pi))
    return float(np.sum(np.log(a + b)))
 
prev = loglik(pi_, mu_, sd_)
print(f"{'iteration':>10} {'log likelihood':>17} {'increase':>12}")
print(f"{0:>10} {prev:>17.6f} {'-':>12}")
worst = 0.0
for it in range(1, 41):
    a = pi_*np.exp(-0.5*((data-mu_[0])/sd_[0])**2)/(sd_[0]*np.sqrt(2*np.pi))
    b = (1-pi_)*np.exp(-0.5*((data-mu_[1])/sd_[1])**2)/(sd_[1]*np.sqrt(2*np.pi))
    g = a/(a+b)                                   # E step = Equation 8.28
    pi_ = float(g.mean())                         # M step
    mu_ = np.array([float(g @ data / g.sum()),
                    float((1-g) @ data / (1-g).sum())])
    sd_ = np.array([float(np.sqrt(g @ (data-mu_[0])**2 / g.sum())),
                    float(np.sqrt((1-g) @ (data-mu_[1])**2 / (1-g).sum()))])
    cur = loglik(pi_, mu_, sd_)
    worst = min(worst, cur - prev)
    if it in (1, 2, 3, 5, 10, 20, 40):
        print(f"{it:>10} {cur:>17.6f} {cur-prev:>12.6f}")
    prev = cur
print(f"\nsmallest step over 40 iterations: {worst:.2e}")
print(f"recovered: pi = {pi_:.4f}, mu = {np.round(mu_,4).tolist()}, "
      f"sd = {np.round(sd_,4).tolist()}")

text

 iteration    log likelihood     increase
         0      -8598.208515            -
         1      -7985.290303   612.918212
         2      -7964.671489    20.618814
         3      -7926.075534    38.595956
         5      -7637.071196   205.985505
        10      -6169.214889    85.216912
        20      -6085.162145     0.000000
        40      -6085.162145     0.000000
 
smallest step over 40 iterations: -9.09e-13
recovered: pi = 0.3507, mu = [-1.5083, 2.0205], sd = [0.7957, 0.495]
true:      pi = 0.3500, mu = [-1.5, 2.0], sd = [0.8, 0.5]

Answer. No. The smallest step over $40$ iterations is $-9.09\times10^{-13}$ — floating-point noise at a converged optimum, not a decrease. This is EM’s guarantee, and it holds without a step size, a line search, or any tuning at all.

Two details worth noticing.

Convergence is not monotone in speed. Iteration 1 gains $612.92$ , iteration 2 gains $20.62$ , iteration 3 gains $38.60$ , iteration 5 gains $205.99$ . The algorithm slows, then accelerates again as the two components separate, then stops dead by iteration 20. The likelihood never falls; the rate does whatever it likes.

Recovery is good but not exact: $\pi = 0.3507$ against $0.3500$ , $\mu = [-1.5083, 2.0205]$ against $[-1.5, 2.0]$ , $\sigma = [0.7957, 0.4950]$ against $[0.8, 0.5]$ . That residual is sampling error at $n = 4000$ , not an EM failure — page 804’s $1/\sqrt{N}$ again.

And the caveat the book attaches: this run converged to the right answer from a deliberately poor start. EM guarantees you will not go downhill; it guarantees nothing about which local maximum you reach. That is the content of “learning in latent-variable models is generally hard, as we will see in Chapter 11.”

🧪 Try It Yourself

Exercise 1 – The three pieces of the risk

Exercise 2 – Design the prior

Exercise 3 – The noise model chooses the estimator

Exercise 4 – The Occam factor, isolated

Exercise 5 – EM never goes downhill

pch.quizTag Did the problems land?

In Problem 1, the degree-9 model had a squared bias of 0.704671 — worse than degree 3's 0.057673. Why does more capacity give more bias?

pch.quizShowAnswer

B — The fit is so unstable across training sets that even its AVERAGE is wrong — Bias is measured on the average prediction across 3000 training sets. With a variance of 2104.81, individual fits swing so wildly that their mean does not settle near the truth on a finite number of trials. The noise floor, meanwhile, is 0.122500 in every row and never moves.
Problem 3 found the Gaussian MLE at 3.420000 and the Laplace MLE at 1.400000 on the same five numbers. What does that show?

pch.quizShowAnswer

B — That least squares is a CONSEQUENCE of assuming Gaussian noise, not a definition of fitting — Four of the five observations lie between 1.0 and 1.5 and the Gaussian estimate sits outside all of them — the honest consequence of having asserted that 12.0 was plausible. Every loss function is a distributional assumption in disguise: squared error is Gaussian, absolute error is Laplace, the ridge penalty is a Gaussian prior, the lasso penalty a Laplace one.
Problem 4 halved the assumed noise scale while keeping the model class correct. What happened to a nominal 95 percent interval?

pch.quizShowAnswer

B — Coverage fell to 0.6697 — it missed one time in three — And the failure is asymmetric: assuming sigma twice too large gives 0.9998, merely useless, while assuming it half too small gives 0.6697, which is wrong. Bayesian inference is calibrated with respect to the model you wrote down, and sigma is part of that model.
In Problem 5, v and w came out dependent given {u, t} at -0.180295, even though t is neither on the path nor the collider. Why?

pch.quizShowAnswer

B — Because t is a DESCENDANT of the collider s, and the head-to-head rule mentions descendants — Conditioning on t was never about t. It leaks information about s, which is enough to open the head-to-head meeting at s. Rows 2, 3 and 4 of that table are the same pair of variables with three different answers, and nothing about v or w changed between them.
Problem 6 split the log evidence into a fit term and an Occam factor. From degree 1 to 11, what happened?

pch.quizShowAnswer

B — The fit improved by 45.70 and the Occam factor worsened by 12.67, so the evidence still favours degree 11 over degree 1 — The decomposition is exact to 5.68e-14. Nobody added the Occam factor; it falls out of doing Equation 8.44's integral. But it grows like a log determinant, far slower than the fit term can improve, which is why page 808 found the evidence penalises extra capacity automatically and gently.
Problem 7 found K = 2 overstated the risk by 0.045981 and K = 40 understated it by 0.006274. Does picking K well fix the flat-versus-nested problem from page 808?

pch.quizShowAnswer

B — No — every K gives an honest estimate of a FIXED model's risk, and none of them repairs the bias from using that estimate to CHOOSE — Choosing K well and nesting the loops are independent decisions. K controls how much data each fold trains on, which is a bias-variance question. Nesting controls whether the number you report was also the number you optimised, which is a different failure entirely — and page 808 measured it at 13.5 percent better than chance on data with no signal at all.

Recall card

The risk decomposes into three pieces: the noise floor, the squared bias of the average fit, and the variance of the fit across training sets. Verified to within Monte Carlo error at degrees 1, 3, 5 and 9.
Underfitting is bias; overfitting is variance. Degree 1 has bias squared 0.456620 and variance 0.052142; degree 9 has variance 2104.812958. The noise floor is 0.122500 in every row and no model removes it.
More capacity can raise bias too. Degree 9’s squared bias, 0.704671, is worse than degree 3’s, because a fit that unstable has a wrong average.
To translate a tuned penalty into a prior, invert the correspondence: tau squared equals sigma squared over N lambda. A penalty of 3e-4 on 25 points is a prior variance of 16.33.
Lambda scales with the NOISE variance, so halving sigma quarters the tau squared you need for the same lambda. Two people using the same lambda on differently noisy data are asserting different beliefs.
Least squares is a consequence of Gaussian noise, not a definition of fitting. Laplace noise gives the median instead: measured, 3.420000 against 1.400000 on the same five numbers. Every loss function is a distributional assumption in disguise.
A wrong noise scale breaks calibration even with the right model class. Halving the assumed sigma takes a nominal 95 percent interval down to 0.6697 coverage. Being too confident is punished far more visibly than being too cautious.
d-separation on a two-collider graph: five queries, five correct. The same pair of variables came out dependent, independent, and dependent again under three conditioning sets, with nothing about the variables changing.
Conditioning on a descendant of a collider opens the path, measured at -0.180295 where conditioning on the collider itself gave -0.383090. Weaker, and nowhere near zero.
The log evidence splits exactly into a best-fit term and an Occam factor, verified to 5.68e-14. From degree 1 to 11 the fit gains 45.70 and the Occam factor costs 12.67.
Nobody adds the Occam factor. It is minus half a log determinant and the prior’s own terms, falling straight out of Equation 8.44’s integral — which is what “automatically embodies a trade-off” means.
It is also gentle. A log determinant grows far more slowly than a fit term improves, which is why the evidence rarely condemns an over-large model outright.
K = 2 overstates the risk by 0.045981 because each fold trains on half the data. By K = 5 the bias is 0.001082 and most of the variance reduction is already captured; K = 40 costs eight times the compute for another 0.0096.
Choosing K and nesting the loops are independent decisions. Every K honestly estimates a fixed model’s risk; none of them repairs the bias created by using that estimate to choose.
EM never decreases the likelihood. Smallest step over 40 iterations: -9.09e-13, which is floating-point noise at a converged optimum. No step size, no line search, no tuning.
The rate is not monotone even though the likelihood is. Iteration 1 gained 612.92, iteration 2 gained 20.62, iteration 5 gained 205.99 as the components separated.
EM guarantees you will not go downhill, and nothing about which local maximum you reach. That is the content of the book’s warning that learning in latent-variable models is generally hard.

Next: the whole chapter on one page. Chapter 8 Formula Sheet

pch.coffeeTagline

pch.coffeeCta

pch.feedbackHeading

pch.feedbackSubheading