""" ============================================================================= 第 3 讲 · 代码 1：Bias–variance —— 重新采样训练集，直接"看见"偏差和方差 ============================================================================= 对应讲义：lectures/lecture03_generalization_regularization.md 第 1 节 ----------------------------------------------------------------------------- 【这个文件要回答的问题 What this file answers】 ----------------------------------------------------------------------------- 第 1 讲 §2.5 只用了**一个**训练集：多项式次数升高，测试误差先降后升。我们说"低次欠拟合（high bias）、 高次过拟合（high variance）"，但 bias 和 variance 是"对训练集的随机性取期望"才有定义的量，一个训练集看不见它们。 合成数据可以无限地重新采样训练集，于是： 演示 1：多项式次数 k = 0..9，每个 k 重新采样 2000 个训练集（n = 30），把期望测试误差拆成 bias² + variance + noise，并检查三项之和等于直接测出来的测试误差。 顺便看到：高次多项式的 variance 是"重尾"的 —— 平均值被极少数倒霉的训练集主导。 演示 2：同一件事换一个旋钮：k = 14、n = 15（第 1 讲测试 MSE 爆到 2×10⁹ 的那个设置）上扫 ridge 的 λ。 λ 从 0 往上加，variance 一路降、bias² 一路升 —— 这回答了第 1 讲 README 的"动手改一改 3"。 演示 3：Approximation error vs estimation error：它们和 bias² / variance 是什么关系？ n 变大时哪一项变、哪一项不变？最优的 k 怎样随 n 变化？ ----------------------------------------------------------------------------- 【运行方式】 ----------------------------------------------------------------------------- python3 01_bias_variance.py # 约 3 秒 依赖本目录的 gen_tools.py（它会自动导入第 1 讲的 datasets.py 和 linear_models.py）。 """ import numpy as np from gen_tools import (bias_variance, fit_least_squares, fit_ridge, make_sine, poly_features, print_header, sine_target, to_unit_interval) NOISE_STD = 0.3 # make_sine 的默认噪声；σ² = 0.09 是任何模型都去不掉的 irreducible error NOISE_VAR = NOISE_STD ** 2 def poly_design(x, k): """多项式特征 [1, u, …, u^k]，u = 2x − 1 ∈ [−1, 1]（和第 1 讲完全相同的特征）。""" return poly_features(to_unit_interval(x), k) def resample_predictions(fit, k, n, R, x_test, rng): """ 重复 R 次："采一个新训练集 → 训练 → 在固定的测试点上预测"。 【参数】 fit ：函数 fit(Phi, y) → θ（比如 fit_least_squares，或者一个包装了 λ 的 ridge） k ：多项式次数；n：每个训练集的大小；R：重复次数 x_test：(m,) 固定的测试点（所有训练集共用，这样才能在每个点上比较不同训练集的预测） rng ：随机数生成器 【返回】(preds, train_mse, theta_norm) preds ：(R, m)，第 r 行是第 r 个训练集训练出的模型在 x_test 上的预测 train_mse ：(R,) 每个模型的训练误差 theta_norm ：(R,) 每个模型的 ‖θ‖ """ Phi_test = poly_design(x_test, k) preds = np.empty((R, len(x_test))) train_mse = np.empty(R) theta_norm = np.empty(R) for r in range(R): x, y = make_sine(n, rng, noise_std=NOISE_STD) Phi = poly_design(x, k) th = fit(Phi, y) preds[r] = Phi_test @ th train_mse[r] = np.mean((Phi @ th - y) ** 2) theta_norm[r] = np.linalg.norm(th) return preds, train_mse, theta_norm def direct_test_mse(preds, f_test, rng): """ 直接测量每个模型的测试误差：给测试点配上**新的**带噪声标签 y = f(x) + ε，算 MSE。 【为什么要单独算一遍？】bias² + variance + σ² 是从公式来的；这里不用公式，只是老老实实地"在新数据上测"， 两者应该对得上（差别只来自有限个噪声样本的随机性）。这是对分解公式的一个数值检验。 【返回】(R,) 每个训练集的模型各自的测试 MSE """ y_test = f_test + rng.normal(0.0, NOISE_STD, size=preds.shape) return np.mean((preds - y_test) ** 2, axis=1) # ============================================================================= # 演示 1：多项式次数 —— bias² 降、variance 升 # ============================================================================= def demo_degree(): print_header("演示 1：多项式次数 k 的 bias–variance 分解（n = 30，每个 k 重采样 R = 2000 个训练集）") rng = np.random.default_rng(31) n, R = 30, 2000 x_test = rng.uniform(0.0, 1.0, size=1000) f_test = sine_target(x_test) print("数据：y = sin(2πx) + N(0, 0.3²)。测试点固定为 1000 个 x；每个 k 都重新采 2000 个训练集。") print("bias² = 平均预测离真相多远；variance = 换一个训练集，预测抖多少；noise = σ² = 0.09。\n") print(f"{'k':>3}{'train':>9}{'bias^2':>11}{'variance':>11}{'sum+noise':>11}{'test(mean)':>12}" f"{'test(median)':>14}{'test(p90)':>11}") rows = {} for k in range(10): preds, tr, _ = resample_predictions(fit_least_squares, k, n, R, x_test, rng) b2, var = bias_variance(preds, f_test) te = direct_test_mse(preds, f_test, rng) rows[k] = (b2, var, float(np.mean(te)), float(np.median(te)), float(np.mean(tr))) print(f"{k:>3}{np.mean(tr):>9.4f}{b2:>11.4g}{var:>11.4g}{b2 + var + NOISE_VAR:>11.4g}{np.mean(te):>12.4g}" f"{np.median(te):>14.4f}{np.quantile(te, 0.9):>11.4f}") best_mean = min(rows, key=lambda k: rows[k][2]) best_med = min(rows, key=lambda k: rows[k][3]) print(f""" 观察 Observations： · "sum+noise" 一列（bias² + variance + 0.09）和直接测出来的 "test(mean)" 一列对得上（差别只来自测试噪声的随机性）： 期望测试误差 = bias² + variance + irreducible noise，这是恒等式。 · bias² 随 k 下降：k = 0 时 {rows[0][0]:.3f}（常数拟合正弦），k = 3 时 {rows[3][0]:.4f}，k = 5 时 {rows[5][0]:.1e}。 variance 随 k 上升：k = 1 时 {rows[1][1]:.3f}，k = 5 时 {rows[5][1]:.3f}，k = 9 时 {rows[9][1]:.3g}。 按平均测试误差，最好的是 k = {best_mean}（{rows[best_mean][2]:.4f}）；按中位数，最好的是 k = {best_med}（{rows[best_med][3]:.4f}）。 · k = 2 比 k = 1 差、k = 4 比 k = 3 差：在中心化的坐标 u = 2x − 1 下，sin(2πx) = −sin(πu) 是奇函数， 偶次项对逼近毫无帮助（bias² 几乎不变：{rows[1][0]:.4f} → {rows[2][0]:.4f}，{rows[3][0]:.4f} → {rows[4][0]:.4f}），却要多估一个参数（variance 变大）。 "多一个参数"不是免费的。 · k ≥ 6 时平均值和中位数分道扬镳（k = 9：平均 {rows[9][2]:.3g}，中位数 {rows[9][3]:.3f}）：variance 是**重尾**的。 大多数训练集上高次多项式还过得去；少数训练集恰好在 [0, 1] 的两端留下空白，多项式在空白处外推、飞出去， 一次就贡献巨大的误差。"期望测试误差"被这些倒霉的训练集主导 —— 这也是第 1 讲里单个训练集看不到的。 · k ≥ 6 的 bias² 不再单调下降，是蒙特卡洛误差：平均预测 f̄ 只是 2000 个样本的平均， bias² 的估计里混进了约 variance / R 的噪声（gen_tools.bias_variance 的【边界】）。真实的 bias² 只会更小。 · 训练误差一路下降（k = 9 时 {rows[9][4]:.4f}，已经低于噪声 0.09）：训练误差衡量的是"拟合了多少"，不是"学到了多少"。""") # ============================================================================= # 演示 2：ridge 的 λ —— 同一个权衡，换一个旋钮 # ============================================================================= def demo_ridge(): print_header("演示 2：k = 14、n = 15 上扫 ridge 的 λ（第 1 讲 §2.5 测试 MSE = 2×10⁹ 的设置）") rng = np.random.default_rng(32) n, k, R = 15, 14, 2000 x_test = rng.uniform(0.0, 1.0, size=1000) f_test = sine_target(x_test) print("15 个参数、15 个训练点：λ = 0 时恰好插值（穿过每个点）。J(θ) = (1/2n)‖Φθ − y‖² + (λ/2)‖θ_{1:}‖²（不惩罚截距）。\n") print(f"{'lambda':>8}{'train':>9}{'bias^2':>11}{'variance':>11}{'sum+noise':>11}{'test(mean)':>12}" f"{'test(median)':>14}{'median|theta|':>15}") rows = {} for lam in [0.0, 1e-6, 1e-4, 1e-3, 3e-3, 1e-2, 3e-2, 0.1, 1.0, 10.0]: fit = (lambda Phi, y, lam=lam: fit_ridge(Phi, y, lam)) # 【Python 小知识】lambda 的默认参数 lam=lam 在定义时就把当前的 λ "冻结"进去； # 不这样写的话，所有 lambda 都会引用循环变量 lam 的最终值（Python 闭包的 late binding）。 preds, tr, nrm = resample_predictions(fit, k, n, R, x_test, rng) b2, var = bias_variance(preds, f_test) te = direct_test_mse(preds, f_test, rng) rows[lam] = (b2, var, float(np.mean(te)), float(np.median(te)), float(np.mean(tr))) print(f"{lam:>8g}{np.mean(tr):>9.4f}{b2:>11.4g}{var:>11.4g}{b2 + var + NOISE_VAR:>11.4g}{np.mean(te):>12.4g}" f"{np.median(te):>14.4g}{np.median(nrm):>15.4g}") best = min((lam for lam in rows if lam > 0), key=lambda lam: rows[lam][2]) print(f""" 观察 Observations： · λ = 0：插值，训练误差 = 0，测试误差的平均值是天文数字（{rows[0.0][2]:.1e}），中位数也有 {rows[0.0][3]:.1e} —— 第 1 讲的 2×10⁹ 是一个"典型"的训练集。 · 加一点点 λ 就天翻地覆：λ = 1e−6 时平均测试误差降到 {rows[1e-6][2]:.3g}，λ = {best:g} 时最低（{rows[best][2]:.4f}）。 variance 从 {rows[1e-6][1]:.3g}（λ = 1e−6）降到 {rows[best][1]:.4f}（λ = {best:g}），代价是 bias² 从 {rows[1e-6][0]:.4f} 升到 {rows[best][0]:.4f}。 · λ 继续增大，bias² 主导：λ = 10 时 bias² = {rows[10.0][0]:.3f}，几乎就是"只剩截距"的常数模型（演示 1 的 k = 0：bias² ≈ 0.49）。 · 所以 λ 和多项式次数 k 是同一个"复杂度旋钮"的两种形式：ridge 不删特征，而是让 θ 不能变大， 曲线就没法在数据点之间剧烈摆动（‖θ‖ 的中位数从 10⁷ 量级降到个位数）。 · 最好的 ridge（{rows[best][2]:.4f}）仍然比演示 1 的最好结果（n = 30、k = 3：约 0.11）差得多：这里只有 15 个点 —— 正则化能控制 variance，但不能凭空变出数据。""") # ============================================================================= # 演示 3：approximation error vs estimation error # ============================================================================= def demo_approx_estimation(): print_header("演示 3：Approximation error vs estimation error —— 哪一项随 n 变？") rng = np.random.default_rng(33) # 用一个很密的均匀网格近似 x ~ U(0, 1) 上的积分 E_x[·]。 grid = (np.arange(20000) + 0.5) / 20000 f_grid = sine_target(grid) x_eval = grid[::10] # 2000 个点用来评估 estimation error ks, ns, R = [1, 3, 5, 7], [30, 100, 300, 1000], 400 print("h*_k = 次数 ≤ k 的多项式里离 sin(2πx) 最近的那个（在 U(0,1) 上的 L² 投影，用 20000 个无噪声网格点算）。") print(" approximation error = E_x[(h*_k − f)²] —— 假设类本身的局限，和 n 无关") print(" estimation error = E_S E_x[(ĥ_S − h*_k)²] —— 有限样本让我们找不到 h*_k，随 n 下降") print(" excess risk = approximation + estimation（勾股定理：ĥ − h*_k 在假设类里，h*_k − f 和假设类正交）") print(f" 每格是 R = {R} 个训练集的平均；括号里是经验公式 (k+1)σ²/n。\n") head = "".join(f"{'n=' + str(n):>20}" for n in ns) print(f"{'k':>3}{'approx':>11}{head}") excess = {} for k in ks: Phi_grid = poly_design(grid, k) h_star = fit_least_squares(Phi_grid, f_grid) # 无噪声、2 万个点：就是 L² 投影 approx = float(np.mean((Phi_grid @ h_star - f_grid) ** 2)) h_star_eval = poly_design(x_eval, k) @ h_star cells = [] for n in ns: preds, _, _ = resample_predictions(fit_least_squares, k, n, R, x_eval, rng) est = float(np.mean((preds - h_star_eval) ** 2)) excess[(k, n)] = (approx, est) cells.append(f"{est:.3g} ({(k + 1) * NOISE_VAR / n:.2g})") print(f"{k:>3}{approx:>11.3g}" + "".join(f"{c:>20}" for c in cells)) print("\n 每个 n 下 excess risk（approx + estimation）最小的 k：") best_k = {} for n in ns: best_k[n] = min(ks, key=lambda k: sum(excess[(k, n)])) detail = ", ".join(f"k={k}: {sum(excess[(k, n)]):.3g}" for k in ks) print(f" n = {n:>4}：best k = {best_k[n]} [{detail}]") a1, e1 = excess[(1, 1000)] print(f""" 观察 Observations： · approximation error 只取决于 k：k = 1 时 {excess[(1, 30)][0]:.3f}（直线逼近正弦），k = 7 时 {excess[(7, 30)][0]:.1e}。 数据再多它也不会变 —— 想降低它只能换更大的假设类。 · estimation error 随 n 大约按 1/n 下降；n ≥ 100、k ≥ 3 时和经验公式 (k+1)σ²/n 很接近： "每多一个参数，就多付 σ²/n 的估计误差"。k = 1 时比公式大（n = 1000：{e1:.2g} vs {2 * NOISE_VAR / 1000:.2g}）， 因为模型错设（misspecified）时，逼近不了的那部分 f − h*（大小 ≈ {a1:.2f}）对估计来说也像噪声： 有效噪声 ≈ σ² + approx = {NOISE_VAR + a1:.2f}，2 × {NOISE_VAR + a1:.2f} / 1000 ≈ {2 * (NOISE_VAR + a1) / 1000:.2g}。 · 最优的 k 随 n 变大（n = {ns[0]} 时 k = {best_k[ns[0]]}，n = {ns[-1]} 时 k = {best_k[ns[-1]]}）：数据多了，付得起更多参数的估计误差， 就该用更复杂的模型。这是"模型复杂度应该随数据量增长"的定量版本。 · n = 30 时 k = 5、7 的 estimation error 又一次被重尾的训练集抬高（k = 7：{excess[(7, 30)][1]:.3g}，远大于公式的 {8 * NOISE_VAR / 30:.2g}）： 经验公式是 n ≫ k 时的渐近结果，n 只有参数个数的 4 倍时并不可靠。 · 和 bias–variance 的关系：bias² ≈ approximation error（n 大时平均预测 → h*）， variance ≈ estimation error。两种分解从不同角度切同一块蛋糕：bias–variance 对训练集取期望， approximation–estimation 以"假设类里最好的那个 h*"为参照点 —— 后者是学习理论（§3）的语言。""") def main(): demo_degree() demo_ridge() demo_approx_estimation() if __name__ == "__main__": main()