""" ============================================================================= 第 4 讲 · 代码 2：Bagging 与随机森林（bagging & random forests） ============================================================================= 对应讲义：lectures/lecture04_trees_ensembles_nn.md 第 3 节 ----------------------------------------------------------------------------- 【这个文件要回答的问题】 ----------------------------------------------------------------------------- 代码 1 说明单棵深树的方差很大。把很多棵树平均起来，方差能降多少？降不下去的那部分是什么？ 演示 1：Bagging 降方差、不降偏差。在 Friedman #1 数据上重复 20 次“抽训练集 → 训练 30 棵 bootstrap 树”， 直接量出单棵树的方差 σ²、树与树之间的相关系数 ρ，检验公式 Var(平均) = ρσ² + (1 − ρ)σ²/B。 对照：bagging 一个稳定的模型（线性回归）几乎什么都不改变。 演示 2：随机森林的特征子采样（max_features = m < d）。m 越小，ρ 越小，方差的地板 ρσ² 越低 —— 但单棵树变差（σ² 和偏差上升）。净效果取决于数据：有一半特征是噪声时，小 m 反而更差。 演示 3：Out-of-bag（OOB）误差。每棵树大约有 36.8% 的样本没见过；用它们估计的误差和真正的测试误差很接近， 而训练误差严重低估。 ----------------------------------------------------------------------------- 【运行方式】 ----------------------------------------------------------------------------- python3 02_bagging_forests.py # 约 25 秒（一共训练约 4000 棵树） 依赖本目录的 tree_models.py、nonlinear_data.py，以及第 1 讲的 linear_models.py。 """ import numpy as np from nonlinear_data import friedman1_target, make_friedman1 from tree_models import RandomForest from linear_models import add_intercept, fit_least_squares, mse, print_header # 模拟实验的规模：R 个独立的训练集，每个训练集上训练 B 棵树。改大会更准，但更慢（时间 ∝ R·B）。 N_TRAIN, R_SETS, B_TREES, N_TEST = 200, 20, 30, 400 _SIM_CACHE = {} # ============================================================================= # 工具：重复实验，测量 σ²、ρ、偏差和方差 # ============================================================================= def simulate_forests(d, m, seed): """ 在 R_SETS 个独立的 Friedman #1 训练集上各训练一个 B_TREES 棵树的森林（每个节点抽 m 个特征）。 【思路】要谈“方差”，就必须有很多个训练集：方差指的是“换一份训练数据，预测会变多少”。 这里每份训练集都从同一个分布里重新抽，所以 P[r, b, i] 的随机性来自两处：训练集 r 和 bootstrap / 特征子采样。 【参数】d：特征数（≥ 5；前 5 个有用）；m：max_features；seed：训练集和测试点的随机种子 【返回】(P, f)：P 形状 (R, B, N_TEST)，第 r 个训练集上第 b 棵树在第 i 个测试点的预测；f 是测试点的真实函数值 【边界】同一组 (d, m, seed) 只算一次（存在 _SIM_CACHE 里），演示 1 和演示 2 可以共用结果。 同一个 d、seed 下不同 m 用的是同一批训练集，所以比较不同 m 时没有“训练集不同”带来的额外噪声。 """ key = (d, m, seed) if key not in _SIM_CACHE: rng = np.random.default_rng(seed) X_test = rng.uniform(0.0, 1.0, size=(N_TEST, d)) train_sets = [make_friedman1(N_TRAIN, rng, d=d) for _ in range(R_SETS)] P = np.zeros((R_SETS, B_TREES, N_TEST)) for r, (X, y) in enumerate(train_sets): forest = RandomForest(B_TREES, "regression", max_features=m, seed=1000 + r).fit(X, y) P[r] = forest.predict_each(X_test) _SIM_CACHE[key] = (P, friedman1_target(X_test)) return _SIM_CACHE[key] def variance_stats(P): """ 从 P (R, B, n) 估计单棵树的方差 σ² 和两棵树之间的协方差 ρσ²（都对测试点取平均）。 【思路】对每个测试点 x： σ²(x) = 单棵树预测的方差（对训练集 r 取方差，再对 B 棵树平均） ρσ²(x) = 同一训练集上两棵**不同**的树的协方差，用恒等式 Var(Σ_b T_b) = Σ_b Var(T_b) + Σ_{b≠b'} Cov(T_b, T_b') 反解出所有 B(B − 1) 对的平均协方差。 ρ 就是 (平均协方差) / (平均方差)。 【返回】(sigma2, rho_sigma2, rho)：三个 float 【ML 联系】两棵树为什么相关？因为它们看的是**同一份**训练集（只是 bootstrap 不同）：训练集里的噪声 会同时影响所有的树，这部分误差平均多少棵树都消不掉。 """ R, B, _ = P.shape var_each = P.var(axis=0) # (B, n)：每棵树各自的方差（对 r） sigma2 = var_each.mean(axis=0) # (n,) var_sum = P.sum(axis=1).var(axis=0) # (n,)：B 棵树之和的方差 cov = (var_sum - var_each.sum(axis=0)) / (B * (B - 1)) s2, c = float(sigma2.mean()), float(cov.mean()) return s2, c, c / s2 def bias2_var(preds, f): """ 预测 preds (R, n)（R 个训练集上的预测）相对真实函数 f (n,) 的偏差² 和方差，对测试点取平均。 【思路】偏差² = (E[预测] − f)²。只有 R 个训练集时，样本均值本身带着方差 / R 的噪声， 直接平方会把偏差高估 var/R；所以减去 s²/R（s² 是 ddof = 1 的样本方差），得到无偏的估计。 R = 20、单棵树方差约 10 时，这个修正约 0.5，不能忽略。 【返回】(bias2, var) """ R = preds.shape[0] mean_pred = preds.mean(axis=0) s2 = preds.var(axis=0, ddof=1) bias2 = (mean_pred - f) ** 2 - s2 / R return float(bias2.mean()), float(preds.var(axis=0).mean()) # ============================================================================= # 演示 1：Bagging 降方差、不降偏差 # ============================================================================= def demo_bagging_variance(): print_header("演示 1：Bagging = 平均 B 棵 bootstrap 树 —— 方差公式 ρσ² + (1 − ρ)σ²/B") d, seed = 10, 21 print(f"Friedman #1（d = {d}，前 5 个特征有用，噪声方差 1），每个训练集 {N_TRAIN} 个点；") print(f"重复 {R_SETS} 个独立训练集，每个训练集上 {B_TREES} 棵不限深度的 bootstrap 回归树（max_features = d，即纯 bagging）。") P, f = simulate_forests(d, d, seed) s2, c, rho = variance_stats(P) print(f"\n单棵树的方差 σ² = {s2:.2f}，两棵树的协方差 ρσ² = {c:.2f}，相关系数 ρ = {rho:.3f}\n") print(f"{'B':>4}{'Var (measured)':>16}{'rho*s2 + (1-rho)*s2/B':>24}{'bias^2':>9}{'MSE vs f':>10}") rows = [] for B in [1, 2, 5, 10, B_TREES]: ens = P[:, :B].mean(axis=1) # (R, n)：前 B 棵树的平均 b2, v = bias2_var(ens, f) pred = c + (s2 - c) / B err = float(np.mean((ens - f) ** 2)) rows.append((B, v, pred, b2, err)) print(f"{B:>4}{v:>16.2f}{pred:>24.2f}{b2:>9.2f}{err:>10.2f}") # 对照：bagging 线性回归（稳定、低方差、高偏差的模型） rng = np.random.default_rng(seed) # 同一个种子 → 和上面完全相同的测试点和训练集 X_test = rng.uniform(0.0, 1.0, size=(N_TEST, d)) train_sets = [make_friedman1(N_TRAIN, rng, d=d) for _ in range(R_SETS)] A_test = add_intercept(X_test) boot_rng = np.random.default_rng(22) ols, bagged = np.zeros((R_SETS, N_TEST)), np.zeros((R_SETS, N_TEST)) for r, (X, y) in enumerate(train_sets): A = add_intercept(X) ols[r] = A_test @ fit_least_squares(A, y) boots = [A_test @ fit_least_squares(A[idx], y[idx]) for idx in (boot_rng.integers(0, N_TRAIN, N_TRAIN) for _ in range(B_TREES))] bagged[r] = np.mean(boots, axis=0) b_ols, v_ols = bias2_var(ols, f) b_bag, v_bag = bias2_var(bagged, f) print(f"\n对照：线性回归（低方差、高偏差）也做同样的 bagging（{B_TREES} 个 bootstrap 样本）：\n") print(f"{'model':<22}{'bias^2':>9}{'variance':>10}") print(f"{'linear regression':<22}{b_ols:>9.2f}{v_ols:>10.3f}") print(f"{'bagged linear reg.':<22}{b_bag:>9.2f}{v_bag:>10.3f}") one, full = rows[0], rows[-1] print(f""" 观察 Observations： · 实测方差和公式几乎重合：B = {rows[1][0]}、{rows[2][0]}、{rows[3][0]} 时 {rows[1][1]:.2f} vs {rows[1][2]:.2f}、{rows[2][1]:.2f} vs {rows[2][2]:.2f}、{rows[3][1]:.2f} vs {rows[3][2]:.2f}。 （B = {B_TREES} 那一行按构造必然相等：ρσ² 就是用这 {B_TREES} 棵树估出来的；真正的检验是 B < {B_TREES} 的几行。）平均 {B_TREES} 棵树把方差从 {one[1]:.2f} 降到 {full[1]:.2f}，但再多的树也降不到 ρσ² = {c:.2f} 以下：两棵树都是从**同一份**训练集里抽的 bootstrap， 训练集里的噪声是它们共有的，平均消不掉。 · 偏差几乎不变（单棵 {one[3]:.2f} → 平均 {full[3]:.2f}）：平均很多个“同样有偏”的模型，结果还是同样有偏。 所以 bagging 的零件应该是**低偏差、高方差**的模型 —— 不剪枝的深树正好如此（ρ = {rho:.3f}，树和树之间相关性很低）。 · 线性回归恰好相反：偏差² {b_ols:.2f} 主导，方差只有 {v_ols:.3f}。bagging 之后偏差 {b_bag:.2f}、方差 {v_bag:.3f}，几乎没变 —— 对一个稳定的、（近似）线性的估计量，bootstrap 平均基本还是它自己。Breiman（1996）的结论：bagging 帮助的是“不稳定”的学习器。 · 最终误差 MSE vs f ≈ 偏差² + 方差：bagging 后 {full[4]:.2f}，主要剩下的是偏差（{full[3]:.2f}）。要降偏差，得换思路（§4 boosting）。""") # ============================================================================= # 演示 2：随机森林 —— 特征子采样降低 ρ # ============================================================================= def demo_random_forest(): print_header("演示 2：Random forest —— 每个节点只看 m 个随机特征，让树“被迫”不同") print(f"同样的重复实验（{R_SETS} 个训练集 × {B_TREES} 棵树）。两份数据：") print(" d = 5 ：Friedman #1 的 5 个特征全部有用") print(" d = 10：再加 5 个纯噪声特征（和演示 1 是同一批训练集）\n") print(f"{'data':<7}{'m':>4}{'sigma^2':>9}{'rho':>8}{'rho*s2':>9}{'bias^2':>9}{'var(B=30)':>11}{'MSE vs f':>10}") results = {} for d, seed, ms in [(5, 23, [5, 3, 1]), (10, 21, [10, 3, 1])]: for m in ms: P, f = simulate_forests(d, m, seed) s2, c, rho = variance_stats(P) ens = P.mean(axis=1) b2, v = bias2_var(ens, f) err = float(np.mean((ens - f) ** 2)) results[d, m] = (s2, rho, c, b2, v, err) print(f"{'d=' + str(d):<7}{m:>4}{s2:>9.2f}{rho:>8.3f}{c:>9.2f}{b2:>9.2f}{v:>11.2f}{err:>10.2f}") r5 = r10 = results # 两个别名，只是为了让下面的 f-string 好读 best5 = min((m for m in [5, 3, 1]), key=lambda m: results[5, m][5]) best10 = min((m for m in [10, 3, 1]), key=lambda m: results[10, m][5]) print(f""" 观察 Observations： · 机制如预期：m 越小，树与树越不相关（d = 10：ρ 从 {r10[10, 10][1]:.3f} 降到 {r10[10, 1][1]:.3f}），方差的地板 ρσ² 从 {r10[10, 10][2]:.2f} 降到 {r10[10, 1][2]:.2f}。 · 代价同样明显：每棵树只能在随机的少数特征里挑切点，单棵树变差 —— σ² 从 {r10[10, 10][0]:.1f} 涨到 {r10[10, 1][0]:.1f}，偏差² 从 {r10[10, 10][3]:.2f} 涨到 {r10[10, 1][3]:.2f}。 · 净效果取决于数据： - d = 5（全部特征有用）：最好的是 m = {best5}（MSE vs f {r5[5, best5][5]:.2f}，纯 bagging {r5[5, 5][5]:.2f}），子采样的好处很小； - d = 10（一半是噪声特征）：最好的是 m = {best10}（{r10[10, best10][5]:.2f}）；m = 1 时一半的节点只能在噪声特征上切，偏差² 翻倍多，MSE {r10[10, 1][5]:.2f}。 这和 ESL 第 15.3 节的结论一致：相关特征占比小时，小 m 的随机森林表现差。 · 为什么好处不大？bagging 之后剩下的误差主要是偏差（d = 10：偏差² {r10[10, 10][3]:.2f} vs 方差 {r10[10, 10][4]:.2f}）。 特征子采样最多只能省下方差那一部分（就算把地板 ρσ² = {r10[10, 10][2]:.2f} 压到 0 也省不了更多），却会推高偏差。 而且深树 + 有噪声的回归，bagging 的 ρ 本来就只有 {r10[10, 10][1]:.2f} 左右。ρ 很高的问题（比如有一两个特别强的特征、 每棵树的第一刀都一样）上，特征子采样的作用更大。Breiman 建议回归用 m ≈ d/3、分类用 m ≈ √d， 但 m 终究是一个要调的超参数（用 OOB 误差调很方便，见演示 3）；sklearn 的 RandomForestRegressor 默认就用全部特征。""") # ============================================================================= # 演示 3：Out-of-bag 误差 # ============================================================================= def demo_oob(): print_header("演示 3：Out-of-bag（OOB）误差 —— 不需要验证集的“免费”测试误差") rng = np.random.default_rng(24) X, y = make_friedman1(300, rng) Xt, yt = make_friedman1(5000, rng) print("Friedman #1（d = 10），训练 300 个点，测试 5000 个点。对每个 B，比较三种误差估计：") print(" train MSE：在训练集上的误差（每棵树都见过这些点）；OOB MSE：每个点只用没见过它的树预测；test MSE：真正的测试误差。\n") print(f"{'max_features':<14}{'B':>5}{'train MSE':>11}{'OOB MSE':>10}{'test MSE':>10}{'OOB - test':>12}") summary = {} for m in [10, 3]: forest = RandomForest(200, "regression", max_features=m, seed=25).fit(X, y) oob_frac = float(np.mean(forest.inbag == 0)) for B in [5, 10, 25, 50, 100, 200]: sub = RandomForest(B, "regression", max_features=m, seed=25) sub.trees, sub.inbag, sub._X_train = forest.trees[:B], forest.inbag[:B], X # 复用前 B 棵树 oob, n_oob = sub.oob_value() ok = n_oob > 0 tr, te = mse(y, sub.predict(X)), mse(yt, sub.predict(Xt)) ob = mse(y[ok], oob[ok]) summary[m, B] = (tr, ob, te, int(np.sum(~ok))) note = "" if ok.all() else f" ({np.sum(~ok)} points never OOB)" print(f"{m:<14}{B:>5}{tr:>11.2f}{ob:>10.2f}{te:>10.2f}{ob - te:>12.2f}{note}") summary[m, "frac"] = oob_frac print(f""" 观察 Observations： · 每棵树没见过的样本比例 = {summary[10, 'frac']:.3f}，和理论值 (1 − 1/n)ⁿ ≈ e⁻¹ ≈ 0.368 一致。 · 训练误差严重低估：200 棵树时 train MSE 只有 {summary[10, 200][0]:.2f}，而测试误差是 {summary[10, 200][2]:.2f}。 每棵树都是长到纯的，它见过的点几乎被“背”了下来（只有 ~1/3 的树没见过某个点，它们的预测才把训练误差拉高一点）。 · OOB 误差和测试误差很接近（B = 200：{summary[10, 200][1]:.2f} vs {summary[10, 200][2]:.2f}；m = 3：{summary[3, 200][1]:.2f} vs {summary[3, 200][2]:.2f}）。 B 小的时候 OOB 偏高（B = 5：{summary[10, 5][1]:.2f} vs {summary[10, 5][2]:.2f}）：每个点只有约 0.37·B 棵 OOB 树，OOB 预测 是“更小的森林”的预测，方差更大 —— 所以 OOB 估计的是一个比实际森林小的模型，偏悲观；B 大了以后这个差别就消失了。 · In practice OOB is free cross-validation (Lecture 3): use it to choose m or min_samples_leaf without holding out a validation set. It only works for bagging-type methods — in boosting every tree depends on all earlier trees, so there are no clean "unseen" examples.""") def main(): demo_bagging_variance() demo_random_forest() demo_oob() if __name__ == "__main__": main()