""" ============================================================================= 第 5 讲 · 数据模块：在第 1 讲的信贷世界上加"代理变量、年龄、有偏的标签"，以及一个可控 base rate 的极简世界 ============================================================================= 对应讲义：lectures/lecture05_measuring_fairness.md（§4、§6、§7、§9 的实验都从这里取数据） ----------------------------------------------------------------------------- 【为什么要一个新的数据模块？ Why a new data module?】 ----------------------------------------------------------------------------- 第 1 讲的 trust_data.py 已经有一个"知道真相"的信贷世界：两个群体（group 1 占 30%）、真实还款概率 p_true、 thin-file 效应（信用历史长度对 group 1 几乎不代表风险）。第 5 讲要"度量公平"，还需要四样东西，第 1 讲没有： 1. 一个**代理变量**（proxy）：邮编 zip。美国历史上的 redlining（按街区划红线拒贷）让居住地和种族高度相关， 所以"模型不看种族、只看邮编"并不能让群体信息消失 —— 这就是 Barocas & Selbst (2016) 说的 redundant encoding。 2. **年龄** age：第二个受保护属性（美国 ECOA 也禁止按年龄歧视信贷申请人，有例外），用来做交叉群体（intersectionality）。 年龄和信用历史长度相关（年轻人的信用档案自然更"薄"）。 3. **有偏的历史标签**： · officer_decisions：过去的信贷员**决定**（批 / 不批）—— 他们对 group 1 额外打了折扣。 如果银行"训练模型去模仿过去的决定"，模型就会学到这个折扣。 · observed_label：记录下来的还款结果是真实结果的一个**有偏代理**（group 1 的按时还款有一部分被记成违约）， 两种误记方式："shift"（按群体平移的代理）和 "flip"（与特征无关的随机误记）。 对应 Obermeyer et al. (2019)：用"医疗花费"代替"医疗需求"。 4. 一个**极简的潜变量世界** make_latent_world：只有一个特征 x，两个群体的 base rate 可以随意设定。 用它扫描 base rate 的差距，干净地看不可能定理（第 1 讲的信贷世界里 base rate 差距是固定的）。 【重要约定】 · 本模块**不修改**第 1 讲的任何文件，只是 import 它们。 · make_credit_plus 先调用第 1 讲的 make_credit(n, rng)，**再**用同一个 rng 抽 zip 和 age。 所以 make_credit_plus(4000, default_rng(103)) 的 X、y、group、p_true 和第 1 讲的部署数据**逐位相同**， 第 5 讲的数字可以直接和第 1 讲 §5.2 对照。 · 每个演示用自己的 np.random.default_rng(seed)。 【numpy 小知识】 rng.choice(k, size=n, p=probs) 从 0..k−1 里按概率 probs 抽 n 次 rng.gamma(shape, scale, size) Gamma 分布（正的、右偏），这里用来生成"毕业 / 工作后多久才开始有信用记录" np.eye(k)[idx] one-hot 编码：第 i 行是单位矩阵的第 idx[i] 行 """ import sys from pathlib import Path import numpy as np # 把第 1 讲的目录加到模块搜索路径里，这样就能 import 第 1 讲的 trust_data / trust_tools。 # 【Python 小知识】Path(__file__).resolve() 是"这个文件的绝对路径"；.parent 是它所在的目录； # .parent.parent / "lecture01" 就是 code/lecture01。sys.path.insert(0, ...) 让 Python 先在那里找模块。 sys.path.insert(0, str(Path(__file__).resolve().parent.parent / "lecture01")) from trust_data import LEGIT, FEATURES, make_credit, true_logit # noqa: E402 （import 必须在 sys.path 修改之后） N_ZIPS = 10 # 10 个"街区"（邮编） GROUP1_ZIPS = [7, 8, 9] # group 1 主要住在这 3 个街区（历史上的居住隔离） AGE_BUCKETS = [(18, 30, "18-29"), (30, 40, "30-39"), (40, 200, "40+")] # [lo, hi) 年龄段 # ============================================================================= # 1. 第 1 讲的世界 + 邮编 + 年龄 # ============================================================================= def make_credit_plus(n, rng, minority_frac=0.3, deployment=True, segregation=0.8): """ 生成 n 个申请人：第 1 讲的全部字段，再加上邮编 zip 和年龄 age。 【思路】 1. 调用第 1 讲的 make_credit(n, rng, ...)：得到 X（6 列）、y、group、p_true。 2. 邮编：以概率 segregation，一个人住在"自己群体常住的街区"里（group 1 → 街区 7–9，group 0 → 街区 0–6，均匀）； 否则住在 10 个街区中的任意一个（均匀）。segregation 越大，邮编越能"猜出"群体。 **邮编本身对还款没有任何因果作用**：真实还款概率只由第 1 讲的公式决定，不含 zip。 3. 年龄：age = 19 + history + Gamma(2, 4)。信用历史不可能比"成年后的时间"长， 所以历史越长的人年龄越大；Gamma 项是"成年后多久才开始有信用记录"（平均 8 年）。 【参数】 n 申请人数 rng np.random.Generator minority_frac group 1 的比例（传给第 1 讲） deployment True = 决策时刻的数据（officer_flag 全为 0；默认，因为第 5 讲不讨论泄露） segregation 0~1，居住隔离的程度；0 = 邮编和群体无关 【返回】第 1 讲的 dict，再加三个键： "zip" 长度 n 的 int 数组，0..9 "age" 长度 n 的 float 数组（岁） "age_bucket" 长度 n 的 int 数组，0/1/2 对应 AGE_BUCKETS 【边界】 · 必须先调 make_credit 再抽 zip/age（顺序决定了 rng 的状态）：这样前四个字段和第 1 讲同 seed 的数据完全一样。 · segregation = 1 时 group 1 只住在 7–9 号街区、group 0 只住在 0–6 号：邮编 = 群体标签的完美代理。 【可信联系】§6：fairness through unawareness 为什么失败。 """ d = make_credit(n, rng, minority_frac=minority_frac, deployment=deployment) g = d["group"] own = np.where(g == 1, rng.choice(GROUP1_ZIPS, size=n), rng.choice([z for z in range(N_ZIPS) if z not in GROUP1_ZIPS], size=n)) anywhere = rng.integers(0, N_ZIPS, n) d["zip"] = np.where(rng.random(n) < segregation, own, anywhere) d["age"] = 19.0 + d["X"][:, 2] + rng.gamma(2.0, 4.0, n) d["age_bucket"] = age_bucket(d["age"]) return d def age_bucket(age): """ 把年龄映射到 AGE_BUCKETS 的下标 0/1/2。 【numpy 小知识】np.digitize(age, [30, 40]) 返回每个年龄落在哪个区间：<30 → 0，30–39 → 1，≥40 → 2。 """ return np.digitize(age, [AGE_BUCKETS[1][0], AGE_BUCKETS[2][0]]) def zip_onehot(zip_codes): """ 邮编的 one-hot 编码（去掉第 0 列作为参照类别，避免和截距完全共线）。 【参数】zip_codes：长度 n 的 int 数组（0..N_ZIPS−1） 【返回】n × (N_ZIPS − 1) 的 0/1 矩阵 【numpy 小知识】np.eye(k)[idx] 取单位矩阵的若干行，正好是 one-hot。 """ return np.eye(N_ZIPS)[zip_codes][:, 1:] def design(d, use_group=False, use_zip=False, cols=None, group_history=False, log_income=False): """ 从一份数据 d 里拼出模型的输入矩阵（"给模型看哪些列"由这个函数一处决定）。 【思路】第 5 讲反复比较"模型看得到什么"：看群体（aware）、不看群体（unaware）、看邮编（proxy）…… 把拼矩阵的逻辑集中在一个函数里，脚本里只需要写 design(d, use_zip=True) 之类。 【参数】 d make_credit_plus 返回的 dict use_group 是否加入群体标签 a 这一列（"aware"） use_zip 是否加入邮编的 one-hot（9 列） cols 用第 1 讲 X 的哪几列；None = 5 个正当特征 LEGIT group_history 是否加入交互项 a × history（让线性模型能学"history 对两个群体的意义不同"） log_income 是否把 income 换成 log(income)（和真实公式一致，模型就几乎"设定正确"） 【返回】(M, names)：n × m 的矩阵和列名列表 【边界】group_history=True 时会自动加入 a 本身（只有交互项没有主效应的模型很难解释）。 """ cols = LEGIT if cols is None else cols X = d["X"][:, cols].copy() names = [FEATURES[j] for j in cols] if log_income and 0 in cols: X[:, cols.index(0)] = np.log(X[:, cols.index(0)]) names[cols.index(0)] = "log_income" blocks, a = [X], d["group"].astype(float) if use_group or group_history: blocks.append(a[:, None]) names.append("group") if group_history: blocks.append((a * d["X"][:, 2])[:, None]) names.append("group×history") if use_zip: blocks.append(zip_onehot(d["zip"])) names += [f"zip{z}" for z in range(1, N_ZIPS)] return np.column_stack(blocks), names # ============================================================================= # 2. 有偏的历史标签 # ============================================================================= def officer_decisions(d, rng, penalty=1.0): """ 过去的信贷员做出的批 / 不批决定（historical human decisions）。 【思路】信贷员对风险的判断大体是对的（看得到真实 logit 里的信息），但对 group 1 额外减去 penalty 个 logit： P(批准) = sigmoid(true_logit − penalty · a) 这模拟 Barocas & Selbst 说的"训练数据反映了过去的偏见"：如果银行把"过去的决定"当成标签， 模型学到的就是"过去的人会怎么决定"，包括他们的偏见。 【参数】d：make_credit_plus 的输出；rng：随机数生成器；penalty：对 group 1 的折扣（logit 单位） 【返回】长度 n 的 0/1 数组（1 = 当年批准了） 【边界】penalty = 0 时信贷员是"公正的"，标签只是真实还款概率的一次随机实现。 """ z = true_logit(d["X"], d["group"]) - penalty * d["group"] return (rng.random(len(z)) < 1.0 / (1.0 + np.exp(-z))).astype(int) def observed_label(d, rng, mode="shift", strength=1.0): """ "记录下来的"还款结果：真实结果的一个有偏代理（biased proxy label / measurement bias）。 【思路】真实结果 y（会不会按时还款）是我们关心的；数据库里记下来的是 y_obs。对 group 1，一部分**按时还款**被记成了"违约" （比如这些社区的贷款更多由只上报负面记录的机构发放、或者小额逾期被更激进地上报）。group 0 的记录是准确的。 y_obs 永远 ≤ y：只会把 1 记成 0，不会反过来。两种误记方式： · mode = "shift"（默认，Obermeyer 型）：记录的结果像是来自一个"对 group 1 打了折扣"的世界： P(y_obs = 1 | x, a) = sigmoid(true_logit − strength · a) 实现方法：对 y = 1 且 a = 1 的人，以概率 1 − sigmoid(z − strength)/sigmoid(z) 把 1 改成 0 （z = true_logit）。可以验证 P(y_obs = 1 | x) = sigmoid(z) · sigmoid(z − strength)/sigmoid(z) = sigmoid(z − strength)。 越接近决策边界的人越容易被误记。代理标签是真实目标的一个"按群体平移过的版本" —— 和 Obermeyer et al. (2019) 的"医疗花费 = 健康需求打了折扣"是同一个结构。 · mode = "flip"（随机误记）：对 y = 1 且 a = 1 的人，以固定概率 strength 把 1 改成 0，与 x 无关。 【参数】 d make_credit_plus 的输出 rng 随机数生成器 mode "shift" 或 "flip" strength shift：logit 折扣（0 = 没有偏差）；flip：误记概率（0~1） 【返回】长度 n 的 0/1 数组 y_obs 【边界】strength = 0 时 y_obs = y。flip 模式下 strength 必须在 [0, 1] 内。 【可信联系】§7：依赖标签的公平指标（separation、sufficiency）会继承标签的偏差；两种误记方式让审计员看到的假象不同（代码 04）。 """ y, a = d["y"], d["group"] if mode == "shift": z = true_logit(d["X"], a) # keep = sigmoid(z − s) / sigmoid(z)，在对数空间里算：log sigmoid(u) = −logaddexp(0, −u)，避免 e^{−z} 上溢。 keep = np.exp(np.logaddexp(0.0, -z) - np.logaddexp(0.0, -(z - strength))) flip_prob = np.where(a == 1, 1.0 - keep, 0.0) elif mode == "flip": flip_prob = np.where(a == 1, strength, 0.0) else: raise ValueError(mode) flip = (y == 1) & (rng.random(len(y)) < flip_prob) return np.where(flip, 0, y) # ============================================================================= # 3. 极简的潜变量世界：base rate 可控 # ============================================================================= def make_latent_world(n_per_group, rng, p0, p1, sep=1.5, sep1=None): """ 一个只有一个特征的世界，用来干净地研究不可能定理（§4）。 【思路】"标签决定特征"的生成方式（class-conditional）： y ~ Bernoulli(p_a) base rate 由群体决定，可以随意设 x | y ~ N(+sep/2, 1) 若 y = 1；N(−sep/2, 1) 若 y = 0 —— **两个群体相同** 这样 x 对两个群体的"信息量"完全一样，唯一的差别是 base rate。sep 越大，x 越能区分 y（AUC 越高）。 如果给了 sep1，group 1 用 sep1（"对 group 1 的测量更粗糙"），用来研究"base rate 相同、信息量不同"的情形。 【参数】 n_per_group 每个群体的人数 rng 随机数生成器 p0, p1 两个群体的 base rate，0 < p < 1 sep 两类均值之差（以标准差为单位） sep1 group 1 的 sep（None = 和 group 0 相同） 【返回】dict：x, y, group, p（每个人所在群体的 base rate）, sep（每个人的 sep） 【边界】p0 = p1 且 sep1 = None 时两个群体**同分布**：任何只依赖 x 的分数都同时满足三个公平准则。 """ group = np.repeat([0, 1], n_per_group) p = np.where(group == 1, p1, p0) s = np.where(group == 1, sep if sep1 is None else sep1, sep) y = (rng.random(2 * n_per_group) < p).astype(int) x = rng.normal(np.where(y == 1, s / 2, -s / 2), 1.0) return {"x": x, "y": y, "group": group, "p": p, "sep": s} def latent_scores(w): """ 在潜变量世界里构造三种分数，每一种"天生"满足一个公平准则。 【思路】 · blind：只看 x 的似然比分数 sigmoid(sep · x)。它对 y 的条件分布 P(s | y, a) 和群体无关， 所以 **separation** 精确成立（总体意义上）。但它没有用到 base rate，所以在 base rate 低的群体里偏高。 · calibrated：贝叶斯后验 P(y = 1 | x, a) = sigmoid(sep · x + logit(p_a))。它在两个群体内部都 calibrated， 所以 **sufficiency** 精确成立。代价是：同样的 x，base rate 低的群体得到更低的分数。 · rank：blind 分数在**本群体内部**的百分位（0~1）。两个群体的分数分布都是均匀分布， 所以 **independence** 精确成立（任何阈值下批准率都相等）。 为什么 sep · x 是似然比？log[N(x; s/2, 1) / N(x; −s/2, 1)] = s · x，正好是两个正态密度之比的对数。 【参数】w：make_latent_world 的输出 【返回】dict：name → 长度 2n 的分数数组（都在 [0, 1] 里） 【numpy 小知识】np.argsort(np.argsort(v)) 得到每个元素的"秩"（0 = 最小）；除以 (n − 1) 变成 0~1 的百分位。 """ lr = w["sep"] * w["x"] blind = 1.0 / (1.0 + np.exp(-lr)) logit_prior = np.log(w["p"] / (1.0 - w["p"])) calibrated = 1.0 / (1.0 + np.exp(-(lr + logit_prior))) rank = np.empty_like(blind) for g in (0, 1): m = w["group"] == g r = np.argsort(np.argsort(blind[m], kind="stable"), kind="stable") rank[m] = r / max(m.sum() - 1, 1) return {"blind (separation)": blind, "calibrated (sufficiency)": calibrated, "within-group rank (independence)": rank} if __name__ == "__main__": # 自检：zip/age 的分布、和第 1 讲数据的一致性、标签偏差的大小。 d = make_credit_plus(4000, np.random.default_rng(103)) d1 = make_credit(4000, np.random.default_rng(103), deployment=True) print("自检 self-check") print(f" 与第 1 讲同 seed 的部署数据逐位相同？ X: {np.array_equal(d['X'], d1['X'])}，" f"y: {np.array_equal(d['y'], d1['y'])}，group: {np.array_equal(d['group'], d1['group'])}") for g in (0, 1): m = d["group"] == g share_g1_zips = np.isin(d["zip"][m], GROUP1_ZIPS).mean() print(f" group {g}: n = {m.sum()}，住在街区 7–9 的比例 = {share_g1_zips:.3f}，平均年龄 = {d['age'][m].mean():.1f}，" f"年龄段人数 = {np.bincount(d['age_bucket'][m], minlength=3).tolist()}") rng = np.random.default_rng(0) yo_s = observed_label(d, rng, mode="shift", strength=1.0) yo_f = observed_label(d, rng, mode="flip", strength=0.25) for g in (0, 1): m = d["group"] == g print(f" group {g}: 真实还款率 = {d['y'][m].mean():.3f}，记录的还款率 shift(1.0) = {yo_s[m].mean():.3f}，" f"flip(0.25) = {yo_f[m].mean():.3f}") w = make_latent_world(5, np.random.default_rng(1), 0.6, 0.4) print(f" 潜变量世界（每组 5 人）：group = {w['group'].tolist()}，y = {w['y'].tolist()}")