""" ============================================================================= CS329T 第 2 讲 · 工具模块 1：扩展的信贷世界 + 一个用 numpy 写的梯度提升树（gradient-boosted trees） ============================================================================= 对应讲义：lectures/lecture02_explanations_shapley.md（§5、§7、§8 的实验都用到这里的模型） ----------------------------------------------------------------------------- 【为什么需要这个模块 Why this module】 ----------------------------------------------------------------------------- 第 1 讲的 Model A / B 都是**线性**逻辑回归：每个特征的贡献就是 θ_j·x̃_j，加起来正好等于 logit，解释"免费"。 第 2 讲要回答的问题是：模型**不是线性**的时候怎么办？银行真正上线的往往是树集成（gradient boosting、random forest）， 它们的输出是很多棵树的叶子值之和，特征之间有大量交互（"收入高 **并且** 没有逾期 → 加分"）。 所以本模块提供： 1. make_credit_ext：在第 1 讲的世界上，再加 4 个"申请表上有、但不影响真实还款"的字段： savings 存款（k$），和 income 强相关（存款多的人通常收入高） utilization 信用卡额度使用率，和 debt_ratio 强相关 accounts 信用账户数，和 history 相关（档案越久，账户越多） app_minute 提交申请时是几点几分的"分钟"（0~59），纯噪声 真实还款概率**完全不变**（仍然是第 1 讲 true_logit 的公式）—— 这 4 个字段只通过"和真正原因相关"携带信息。 这样我们就有了**相关特征**（correlated features），第 7 节"把一个特征去掉是什么意思"需要它。 2. BoostedTrees：一个从零写的、很小的梯度提升决策树（XGBoost 风格的二阶 gain），逻辑损失，深度 3。 不追求速度和花样，只追求：读得懂、结果可复现、输出是一个真正"非线性 + 有交互"的黑盒。 ----------------------------------------------------------------------------- 【约定 Conventions】 ----------------------------------------------------------------------------- · 扩展世界的 X 是 n × 10：前 6 列和第 1 讲完全相同（FEATURES），后 4 列见 EXT_FEATURES[6:]。 同一个 seed 下，前 6 列、y、group、p_true 和第 1 讲的 make_credit **逐位相同**（先调用第 1 讲的生成器，再抽新字段）。 · FORM = 申请表上的 9 个字段（扩展世界里除 officer_flag 以外的全部列）；LEGIT（第 1 讲）= 前 5 个"正当"字段。 · 模型都有 logit(X) / proba(X) / predict(X) 三个方法，接收**完整的** X（n × 6 或 n × 10）， 自己用 cols 选出要读的列 —— 这样"模型读了哪些列"是模型的属性，解释代码不需要知道。 【注释标签】【思路】【参数】【返回】【边界】【可信联系】【numpy 小知识】【Python 小知识】同第 1 讲。 """ import sys from pathlib import Path import numpy as np # 把第 1 讲的目录加进模块搜索路径，导入共享的数据生成器和工具。 # 【Python 小知识】Path(__file__).resolve() 是本文件的绝对路径；.parent 是所在目录（lecture02）， # 再 .parent 是 code/；拼上 "lecture01" 就是第 1 讲的代码目录。sys.path.insert(0, ...) 把它放在搜索路径最前面。 sys.path.insert(0, str(Path(__file__).resolve().parent.parent / "lecture01")) from trust_data import FEATURES, FLAG, LEGIT, make_credit # noqa: E402 （第 1 讲的模块） from trust_tools import sigmoid # noqa: E402 EXT_FEATURES = FEATURES + ["savings", "utilization", "accounts", "app_minute"] FORM = [0, 1, 2, 3, 4, 6, 7, 8, 9] # 申请表上的 9 个字段（不含 officer_flag，第 5 列） def make_credit_ext(n, rng, deployment=False, minority_frac=0.3): """ 生成扩展的信贷数据：第 1 讲的 6 列 + 4 个"相关但不是原因"的字段。 【思路】 1. 先调用第 1 讲的 make_credit(n, rng, ...)：它消耗 rng 的前一段随机数，所以前 6 列和第 1 讲同 seed 时完全一样； 2. 再用同一个 rng 抽 4 个新字段。每个新字段 = "某个真实特征的函数" + 独立噪声： savings = income × 0.25 × lognormal(0, 0.3) → log(savings) 和 log(income) 的相关系数约 0.8 utilization = clip(0.40 + 0.9·(debt_ratio − 0.35) + N(0, 0.07), 0, 1) → 和 debt_ratio 相关约 0.84 accounts = Poisson(1 + 0.5·history) → 和 history 相关 app_minute = 0~59 的均匀整数 → 和一切都无关 3. 真实还款概率 p_true 和标签 y 都来自第 1 讲的公式，**不依赖**新字段。 所以从"世界"的角度，新字段对还款没有任何因果作用；但它们和真正的原因相关，模型可能会用它们。 【参数】n：人数；rng：np.random.Generator；deployment：同第 1 讲（True 时 officer_flag 全为 0）；minority_frac：同第 1 讲 【返回】dict：X（n × 10）、y、group、p_true（同第 1 讲） 【边界】同一个 rng 调用两次会得到不同的数据；想要同一份数据，就用同一个 seed 新建 rng。 【可信联系】相关特征是 Shapley 解释最难的地方（讲义 §7）："去掉 income"时，savings 该不该跟着变？ 不同的回答（interventional vs conditional）会给出不同的拒绝理由。 """ d = make_credit(n, rng, minority_frac=minority_frac, deployment=deployment) X = d["X"] income, debt, hist = X[:, 0], X[:, 1], X[:, 2] savings = income * 0.25 * rng.lognormal(0.0, 0.3, n) utilization = np.clip(0.40 + 0.9 * (debt - 0.35) + rng.normal(0.0, 0.07, n), 0.0, 1.0) accounts = rng.poisson(1.0 + 0.5 * hist).astype(float) app_minute = rng.integers(0, 60, n).astype(float) d["X"] = np.column_stack([X, savings, utilization, accounts, app_minute]) return d class BoostedTrees: """ 一个很小的梯度提升决策树分类器（gradient-boosted decision trees, GBDT），逻辑损失，numpy 实现。 【用法】 model = BoostedTrees(cols=LEGIT).fit(X_train, y_train) model.proba(X) # 还款概率 model.logit(X) # 所有树的叶子值之和 + 初始值（"原始分数"） model.predict(X) # 0/1 决策（阈值 0.5） 【思路】Boosting = 一棵接一棵地种小树，每棵新树去拟合"当前模型还错在哪里"： 1. 初始分数 F₀ = log(ȳ / (1 − ȳ))（所有人都预测平均还款率）； 2. 第 t 轮：用当前分数算每个样本的梯度 g = p − y 和二阶导 h = p(1 − p)（逻辑损失对 F 的导数）； 3. 种一棵深度为 depth 的树：在每个结点，枚举"特征 j、阈值 τ"，选使 XGBoost gain 最大的切分： gain = G_L²/(H_L + λ) + G_R²/(H_R + λ) − G²/(H + λ) （G、H 是结点里 g、h 的和）。叶子的值 = −G/(H + λ)，就是这一片样本上"牛顿步"的最优增量； 4. F ← F + lr × 新树的输出。重复 n_trees 轮。 为什么候选阈值只用分位数（n_bins 个）？这是 LightGBM / XGBoost 的 histogram 技巧：把每个特征预先分箱， 之后每个结点只需要 np.bincount 一次就能得到所有切分的 G_L、H_L，比逐个阈值扫描快得多。 【参数】 n_trees 树的棵数；depth：每棵树的深度（深度 3 = 最多 8 个叶子，能表示 3 个特征之间的交互） lr 学习率（shrinkage），越小越需要更多的树 lam 叶子值的 L2 正则 λ；min_leaf：每个叶子至少多少个样本 n_bins 每个特征的候选阈值个数（分位数） cols 读 X 的哪些列；None = 全部 【边界】 · 树是分段常数函数：在两个阈值之间，改变一个特征**完全不会**改变输出 —— 梯度处处为 0 或不存在。 这正是"梯度类解释"（第 3 讲）对树模型无用、而 Shapley 这类"扰动类解释"有用的原因。 · 训练集以外的组合（比如 income 很低但 savings 很高）落在哪个叶子，取决于训练时恰好怎么切的 —— 树在"没见过的区域"的行为是任意的外推。讲义 §7.3 会看到这件事对解释的影响。 · 某个结点找不到合法切分时（样本太少或 gain ≤ 0），它把所有样本送到左孩子（阈值 = +inf）， 这样所有树都是满二叉树，预测时可以统一处理。 【可信联系】Merrick & Taly（2020）的 Explanation Game 案例研究、Lundberg et al. 的 TreeSHAP 都以 GBDT 为对象： 它是表格数据上最常用的"黑盒"，也是信贷模型里最常见的非线性模型之一。 """ def __init__(self, n_trees=150, depth=3, lr=0.05, lam=1.0, min_leaf=80, n_bins=32, cols=None): self.n_trees = n_trees self.depth = depth self.lr = lr self.lam = lam self.min_leaf = min_leaf self.n_bins = n_bins self.cols = cols def _select(self, X): """选出模型要读的列。""" return X if self.cols is None else X[:, self.cols] def _grow_tree(self, Z, bins, g, h): """ 种一棵满二叉树（深度 self.depth），返回 (feat, thr, val) 三个数组。 【思路】结点按"层序"（BFS）编号：根是 0，结点 k 的左右孩子是 2k+1、2k+2；内部结点 0 … 2^depth − 2， 叶子接在后面。按编号从小到大处理，保证先处理父结点、再处理孩子。 node[i] 记录第 i 个样本当前在哪个结点。 【参数】Z：n × d 特征；bins：长度 d 的列表，bins[j][i] = 样本 i 的第 j 个特征落在哪个箱；g、h：梯度和二阶导 【返回】feat（内部结点用哪个特征）、thr（阈值）、val（每个叶子的值） 【numpy 小知识】np.bincount(b, weights=w, minlength=K)：对每个箱 k，求 b == k 的样本的 w 之和。 一次调用就得到所有箱的 G；再 np.cumsum 累加，就是"阈值取第 k 个时，左边的 G_L"。 """ n, d = Z.shape n_int = 2 ** self.depth - 1 feat = np.zeros(n_int, dtype=int) thr = np.full(n_int, np.inf) # +inf：所有样本都去左边（"不切分"） node = np.zeros(n, dtype=int) for k in range(n_int): idx = np.where(node == k)[0] best_gain, best_j, best_b = 1e-12, -1, -1 if len(idx) >= 2 * self.min_leaf: gk, hk = g[idx], h[idx] G, H = gk.sum(), hk.sum() parent = G * G / (H + self.lam) for j in range(d): nb = len(self.thresholds[j]) + 1 if nb < 2: continue bj = bins[j][idx] GL = np.cumsum(np.bincount(bj, weights=gk, minlength=nb))[:-1] HL = np.cumsum(np.bincount(bj, weights=hk, minlength=nb))[:-1] NL = np.cumsum(np.bincount(bj, minlength=nb))[:-1] GR, HR, NR = G - GL, H - HL, len(idx) - NL gain = GL ** 2 / (HL + self.lam) + GR ** 2 / (HR + self.lam) - parent gain[(NL < self.min_leaf) | (NR < self.min_leaf)] = -np.inf b = int(np.argmax(gain)) if gain[b] > best_gain: best_gain, best_j, best_b = float(gain[b]), j, b if best_j >= 0: feat[k] = best_j thr[k] = self.thresholds[best_j][best_b] go_right = bins[best_j][idx] > best_b # 箱号 > b ⇔ x > 阈值 else: go_right = np.zeros(len(idx), dtype=bool) node[idx] = 2 * k + 1 + go_right leaf = node - n_int val = np.zeros(n_int + 1) for l in range(n_int + 1): m = leaf == l if m.any(): val[l] = -g[m].sum() / (h[m].sum() + self.lam) return feat, thr, val def _tree_output(self, tree, Z): """ 一棵树对每一行的输出。 【思路】从根出发走 depth 步：每一步看当前结点的特征是否 > 阈值，> 就去右孩子（2k+2），否则左孩子（2k+1）。 【numpy 小知识】Z[np.arange(n), f] 是"花式索引"（fancy indexing）：第 i 行取第 f[i] 列， 一次取出每个样本在"它当前结点的切分特征"上的值。 """ feat, thr, val = tree n = len(Z) node = np.zeros(n, dtype=int) rows = np.arange(n) for _ in range(self.depth): node = 2 * node + 1 + (Z[rows, feat[node]] > thr[node]) return val[node - (2 ** self.depth - 1)] def fit(self, X, y): """ 训练。 【参数】X：n × D 原始特征（模型只读 cols 指定的列）；y：0/1 标签 【返回】self（可以链式调用） 【边界】某一列只有一个取值时没有候选阈值（nb < 2），这一列永远不会被用来切分。 """ Z = self._select(X).astype(float) n, d = Z.shape qs = np.linspace(0.0, 1.0, self.n_bins + 1)[1:-1] self.thresholds = [np.unique(np.quantile(Z[:, j], qs)) for j in range(d)] # side="left"：箱号 b = 第一个 ≥ x 的阈值的下标，所以 x ≤ thr[b'] ⇔ b ≤ b'。 bins = [np.searchsorted(self.thresholds[j], Z[:, j], side="left") for j in range(d)] p0 = float(np.clip(np.mean(y), 1e-6, 1 - 1e-6)) self.base = np.log(p0 / (1.0 - p0)) F = np.full(n, self.base) self.trees = [] for _ in range(self.n_trees): p = sigmoid(F) g, h = p - y, p * (1.0 - p) tree = self._grow_tree(Z, bins, g, h) F = F + self.lr * self._tree_output(tree, Z) self.trees.append(tree) return self def logit(self, X, chunk=100000): """ 原始分数 F(x) = F₀ + lr × Σ_t tree_t(x)。 【边界】行数很多时（解释算法会一次送进几十万行"混合输入"），按 chunk 行分块计算，避免占用太多内存。 """ Z = self._select(X) out = np.empty(len(Z)) for s in range(0, len(Z), chunk): Zc = Z[s:s + chunk] F = np.full(len(Zc), self.base) for tree in self.trees: F += self.lr * self._tree_output(tree, Zc) out[s:s + chunk] = F return out def proba(self, X): """还款概率 P(y = 1 | x) = sigmoid(logit)。""" return sigmoid(self.logit(X)) def predict(self, X, thr=0.5): """0/1 决策 1{p ≥ thr}。有了这个方法，第 1 讲的 permutation_importance 可以直接用在树模型上。""" return (self.proba(X) >= thr).astype(int) def used_features(self): """ 返回模型**实际用来切分**的特征（在完整 X 里的列号），以及每个特征被用了多少次。 【思路】遍历所有树的内部结点，thr < inf 的结点才真正切分了。 【可信联系】"模型读了这一列" ≠ "模型用了这一列"：一个从未被用来切分的特征，对输出没有任何影响（dummy player）。 """ counts = {} for feat, thr, _ in self.trees: for f, t in zip(feat, thr): if np.isfinite(t): col = int(f) if self.cols is None else int(self.cols[int(f)]) counts[col] = counts.get(col, 0) + 1 return dict(sorted(counts.items())) if __name__ == "__main__": # 自检：扩展世界的前 6 列和第 1 讲同 seed 时相同；新字段的相关性；树模型的准确率。 from trust_tools import accuracy, auc a = make_credit(2000, np.random.default_rng(5)) b = make_credit_ext(2000, np.random.default_rng(5)) print("自检 self-check") print(f" 前 6 列与第 1 讲完全相同：{bool(np.array_equal(a['X'], b['X'][:, :6]))}；y 相同：{bool(np.array_equal(a['y'], b['y']))}") X = b["X"] print(f" corr(log income, log savings) = {np.corrcoef(np.log(X[:, 0]), np.log(X[:, 6]))[0, 1]:.3f}，" f"corr(debt_ratio, utilization) = {np.corrcoef(X[:, 1], X[:, 7])[0, 1]:.3f}，" f"corr(history, accounts) = {np.corrcoef(X[:, 2], X[:, 8])[0, 1]:.3f}") tr = make_credit_ext(4000, np.random.default_rng(6)) te = make_credit_ext(4000, np.random.default_rng(7)) m = BoostedTrees(cols=FORM).fit(tr["X"], tr["y"]) print(f" BoostedTrees（9 个表格字段）测试准确率 = {accuracy(te['y'], m.predict(te['X'])):.3f}，" f"AUC = {auc(te['y'], m.proba(te['X'])):.3f}；oracle 准确率 = {accuracy(te['y'], (te['p_true'] >= 0.5).astype(int)):.3f}") print(f" 用来切分的列及次数：{ {EXT_FEATURES[k]: v for k, v in m.used_features().items()} }") print(f" FLAG 列号 = {FLAG}，LEGIT = {LEGIT}，FORM = {FORM}")