""" ============================================================================= 第 2 讲 · 代码 1：GDA vs 逻辑回归（generative vs discriminative） ============================================================================= 对应讲义：lectures/lecture02_generative_kernels_svm.md 第 1、2 节 ----------------------------------------------------------------------------- 【这个文件要回答的问题】 ----------------------------------------------------------------------------- 演示 1：GDA 只用两个高斯密度和 Bayes 公式，它的后验 p(y=1|x) 却**逐位等于**一个 sigmoid(θᵀx)。 这是第 1 讲 §3.2 预告的“sigmoid 的第三个来源”。反过来不成立：Poisson 的类条件分布 也给出逻辑函数形式的后验 —— 逻辑回归的假设比 GDA 弱。 演示 2：GDA 的假设**成立**时，谁更省数据？（学习曲线：每类 10 ~ 1000 个样本，各重复 200 次） 演示 3：GDA 的假设**不成立**时（正类有一条“单侧长尾”），谁更稳健？GDA 的错误能不能靠更多数据消除？ ----------------------------------------------------------------------------- 【运行方式】 ----------------------------------------------------------------------------- python3 01_gda_vs_logistic.py # 约 1 秒 依赖本目录的 generative_models.py，以及第 1 讲的 linear_models.py（自动加入搜索路径）。 """ import numpy as np from generative_models import (bayes_error_contaminated, bayes_error_gaussian, fit_gda, gda_predict_proba, gda_to_logistic, linear_rule_error, make_contaminated_classes, make_gaussian_classes, make_gda_truth) from linear_models import add_intercept, fit_logistic_newton, print_header, sigmoid def fmt_vec(v, digits=3): """把向量打印成 [a, b, c] 的形式（固定小数位，方便对齐比较）。""" return "[" + ", ".join(f"{x:.{digits}f}" for x in v) + "]" # ============================================================================= # 演示 1：GDA 的后验 = 逻辑回归 # ============================================================================= def demo_posterior_is_logistic(): print_header("演示 1：GDA 的后验恰好是一个逻辑回归 —— sigmoid 的第三个来源") rng = np.random.default_rng(11) truth = make_gda_truth(d=2, delta=2.0, rho=0.5, phi=0.5) X, y = make_gaussian_classes(200, truth, rng) params = fit_gda(X, y) print("数据：x | y ~ N(μ_y, Σ)，d = 2，两类共享 Σ（相关系数 0.5），马氏距离 Δ = 2，每类 200 个。") print("GDA 的 MLE 是封闭解（数一数、求平均），不需要迭代：\n") print(f"{'param':>8}{'estimate':>36}{'truth':>36}") print(f"{'phi':>8}{params['phi']:>36.3f}{truth['phi']:>36.3f}") print(f"{'mu0':>8}{fmt_vec(params['mu0']):>36}{fmt_vec(truth['mu0']):>36}") print(f"{'mu1':>8}{fmt_vec(params['mu1']):>36}{fmt_vec(truth['mu1']):>36}") print(f"{'Sigma':>8}{fmt_vec(params['Sigma'].ravel()):>36}{fmt_vec(truth['Sigma'].ravel()):>36}") # 用 Bayes 公式算后验（两个高斯密度之比），再和 sigmoid(θᵀx) 比较 X_test, _ = make_gaussian_classes(5000, truth, rng) th_gda = gda_to_logistic(params) p_bayes = gda_predict_proba(params, X_test) p_sig = sigmoid(add_intercept(X_test) @ th_gda) print("\n把 GDA 参数换算成 θ = Σ⁻¹(μ₁ − μ₀)、θ₀ = −½μ₁ᵀΣ⁻¹μ₁ + ½μ₀ᵀΣ⁻¹μ₀ + log(φ/(1−φ))，") print(f"在 10000 个新点上比较 “Bayes 公式的后验” 和 “sigmoid(θ₀ + θᵀx)”：max |diff| = " f"{np.max(np.abs(p_bayes - p_sig)):.1e}") th_lr, it_lr = fit_logistic_newton(add_intercept(X), y) th_true = gda_to_logistic(truth) print("\n同一份训练数据上，三个“逻辑回归”的 θ = [θ₀, θ₁, θ₂] 和它们在总体上的错误率（精确值）：\n") print(f"{'model':<40}{'theta':>26}{'error rate':>12}") rows = [("GDA (closed form, then convert)", th_gda), (f"logistic regression ({it_lr} Newton steps)", th_lr), ("truth (Bayes-optimal rule)", th_true)] for name, th in rows: print(f"{name:<40}{fmt_vec(th):>26}{linear_rule_error(th, truth):>12.4f}") # (b) Poisson 的类条件分布：后验也是逻辑函数 print("\n(b) 换一个生成模型：x_j | y ~ Poisson(λ_{y,j})，各坐标独立，φ = p(y=1) = 0.3。") lam0, lam1, phi = np.array([1.0, 5.0, 3.0]), np.array([3.0, 6.0, 3.0]), 0.3 rng_b = np.random.default_rng(12) yb = (rng_b.uniform(size=5000) < phi).astype(int) Xb = rng_b.poisson(np.where(yb[:, None] == 1, lam1, lam0)) # log x! 的查表：logfact[k] = log(k!) = Σ_{i=1}^{k} log i logfact = np.concatenate([[0.0], np.cumsum(np.log(np.arange(1, Xb.max() + 1)))]) lp0 = np.sum(Xb * np.log(lam0) - lam0 - logfact[Xb], axis=1) + np.log(1 - phi) lp1 = np.sum(Xb * np.log(lam1) - lam1 - logfact[Xb], axis=1) + np.log(phi) post = np.exp(lp1 - np.logaddexp(lp0, lp1)) th_pois = np.concatenate([[np.log(phi / (1 - phi)) - np.sum(lam1 - lam0)], np.log(lam1 / lam0)]) print(f" Bayes 后验 vs sigmoid(θᵀx)，θ = [log(φ/(1−φ)) − Σ(λ₁ − λ₀), log(λ₁/λ₀)] = {fmt_vec(th_pois)}") print(f" 在 5000 个样本上 max |diff| = {np.max(np.abs(post - sigmoid(add_intercept(Xb) @ th_pois))):.1e}") print(f""" 观察 Observations： · GDA 的后验是用两个高斯密度和 Bayes 公式算出来的，代码里完全没有 sigmoid；它却和 sigmoid(θ₀ + θᵀx) 一致到 {np.max(np.abs(p_bayes - p_sig)):.0e}（舍入误差）。原因：共享 Σ 时 log-odds 里的二次项 −½xᵀΣ⁻¹x 两边抵消， 只剩 x 的线性函数。这就是第 1 讲 §3.2 说的“sigmoid 的第三个来源”。 · GDA 换算出的 θ 和直接拟合的逻辑回归 θ 不一样（两个不同的估计量），但都在真值附近； 总体错误率 {linear_rule_error(th_gda, truth):.4f} vs {linear_rule_error(th_lr, truth):.4f}， Bayes 最优是 {linear_rule_error(th_true, truth):.4f}（= Φ(−Δ/2) = {bayes_error_gaussian(2.0):.4f}）。 · (b) 类条件分布换成 Poisson，后验**也**是 sigmoid(θᵀx)（只是 θ 的公式不同）。任何“同一个指数族、 只有自然参数随类别变化”的类条件分布都是这样。所以：GDA ⇒ 逻辑回归，反过来不成立。 逻辑回归的假设更弱 —— 它不在乎 p(x|y) 长什么样，只要 log-odds 是线性的。""") # ============================================================================= # 演示 2、3：学习曲线 # ============================================================================= def learning_curve(truth, sizes, reps, rng, frac=0.0, shift=0.0, l2=1e-3): """ 对每个训练集大小重复 reps 次：生成数据 → 拟合 GDA 和逻辑回归 → 算两者的**精确**总体错误率。 【参数】truth：真实参数；sizes：每类样本数的列表；reps：每个大小重复几次；rng； frac、shift：0 表示 GDA 假设成立（make_gaussian_classes），否则用 make_contaminated_classes； l2：逻辑回归的 L2 系数。只加一点点（1e−3）：样本很少时数据常常线性可分，λ = 0 的 MLE 不存在（第 1 讲 §3.5）。 【返回】(gda_err, lr_err)：两个形状 (len(sizes),) 的数组，每个元素是 reps 次的平均错误率 """ gda_err, lr_err = [], [] for npc in sizes: eg, el = [], [] for _ in range(reps): if frac > 0: X, y = make_contaminated_classes(npc, truth, rng, frac=frac, shift=shift) else: X, y = make_gaussian_classes(npc, truth, rng) eg.append(linear_rule_error(gda_to_logistic(fit_gda(X, y)), truth, frac, shift)) th, _ = fit_logistic_newton(add_intercept(X), y, l2=l2) el.append(linear_rule_error(th, truth, frac, shift)) gda_err.append(np.mean(eg)) lr_err.append(np.mean(el)) return np.array(gda_err), np.array(lr_err) SIZES = [10, 15, 25, 50, 100, 250, 1000] REPS = 200 def demo_assumptions_hold(): print_header("演示 2：GDA 的假设成立时 —— 谁更省数据？") rng = np.random.default_rng(21) truth = make_gda_truth(d=10, delta=2.5, rho=0.5) bayes = bayes_error_gaussian(2.5) print("x | y ~ N(μ_y, Σ)，d = 10，Σ_jk = 0.5^|j−k|，马氏距离 Δ = 2.5。GDA 的模型完全正确，逻辑回归也正确（后验确实是 sigmoid）。") print(f"每个训练集大小重复 {REPS} 次；错误率是**精确的总体错误率**（不是测试集估计），Bayes 最优 = Φ(−1.25) = {bayes:.4f}。\n") g, l = learning_curve(truth, SIZES, REPS, rng) print(f"{'n per class':>12}{'GDA error':>12}{'LR error':>12}{'GDA excess':>12}{'LR excess':>12}{'LR/GDA excess':>15}") for k, npc in enumerate(SIZES): print(f"{npc:>12}{g[k]:>12.4f}{l[k]:>12.4f}{g[k] - bayes:>12.4f}{l[k] - bayes:>12.4f}" f"{(l[k] - bayes) / (g[k] - bayes):>15.2f}") win = [npc for k, npc in enumerate(SIZES) if g[k] < l[k]] lose = [npc for k, npc in enumerate(SIZES) if g[k] >= l[k]] ratio = [(l[k] - bayes) / (g[k] - bayes) for k, npc in enumerate(SIZES) if npc >= 25] print(f""" 观察 Observations： · “excess” = 错误率 − Bayes 最优，也就是“因为数据有限而多犯的错”。两者都随 n 趋于 0：假设都对，都是一致的估计。 · 每类样本数为 {'、'.join(map(str, win))} 时 GDA 更好；每类 ≥ 25 个样本时，逻辑回归多犯的错是 GDA 的 {min(ratio):.2f}～{max(ratio):.2f} 倍，而且这个比值不随 n 趋于 1。这就是 CS229 讲义说的 “模型正确时 GDA 是 asymptotically efficient（渐近有效）的”：它用上了“x 是高斯的”这条额外信息， 同样多的数据估得更准。 · 但在最小的样本量（每类 {'、'.join(map(str, lose))} 个，总共 {2 * lose[0] if lose else 0} 个点、d = 10）GDA 反而更差： 它要估计 Σ 的 55 个元素，样本太少时 Σ̂ 很不准，Σ̂⁻¹ 会放大误差。“渐近”二字很重要 —— 样本极少时，参数少的模型（或者加正则化）才是关键。""" if lose else f""" 观察 Observations： · “excess” = 错误率 − Bayes 最优，也就是“因为数据有限而多犯的错”。两者都随 n 趋于 0：假设都对，都是一致的估计。 · 所有样本量上 GDA 都更好：它用上了“x 是高斯的”这条额外信息。""") def demo_assumptions_violated(): print_header("演示 3：GDA 的假设不成立时 —— 正类有一条“单侧长尾”") rng = np.random.default_rng(31) truth = make_gda_truth(d=10, delta=2.5, rho=0.5) frac, shift = 0.1, 5.0 bayes = bayes_error_contaminated(truth, frac, shift, np.random.default_rng(32)) print(f"负类仍是 N(μ₀, Σ)；正类中 {frac:.0%} 的样本被额外推远 {shift:g}·(μ₁ − μ₀)（落在离边界很远、毫无争议的一侧，") print("相当于第 1 讲 §3.1 的“特别大的肿瘤”）。p(x|y=1) 不再是高斯，但边界附近的后验几乎还是 sigmoid。") print(f"Bayes 最优错误率（Monte Carlo 估计）= {bayes:.4f}。\n") g, l = learning_curve(truth, SIZES, REPS, rng, frac=frac, shift=shift) print(f"{'n per class':>12}{'GDA error':>12}{'LR error':>12}") for k, npc in enumerate(SIZES): print(f"{npc:>12}{g[k]:>12.4f}{l[k]:>12.4f}") # “n = ∞”：在一个巨大的样本上拟合一次，看两个估计量各自收敛到哪里 Xbig, ybig = make_contaminated_classes(100000, truth, np.random.default_rng(33), frac=frac, shift=shift) p_big = fit_gda(Xbig, ybig) th_g = gda_to_logistic(p_big) th_l, _ = fit_logistic_newton(add_intercept(Xbig), ybig, l2=1e-3) eg, el = linear_rule_error(th_g, truth, frac, shift), linear_rule_error(th_l, truth, frac, shift) print(f"{'100000':>12}{eg:>12.4f}{el:>12.4f} ← 近似 n = ∞") g0, g1 = linear_rule_error(th_g, truth, frac, shift, per_class=True) l0, l1 = linear_rule_error(th_l, truth, frac, shift, per_class=True) print(f"\n每类 100000 个样本时，分开看两类的错误率：GDA P(err|y=0) = {g0:.4f}, P(err|y=1) = {g1:.4f}") print(f"{'':>39}LR P(err|y=0) = {l0:.4f}, P(err|y=1) = {l1:.4f}") shift_mu = np.linalg.norm(p_big["mu1"] - truth["mu1"]) / np.linalg.norm(truth["mu1"] - truth["mu0"]) print(f"GDA 估计的 μ₁ 偏离“高斯部分”的中心 {shift_mu:.2f}·‖μ₁ − μ₀‖（= frac × shift = {frac * shift:.2f}，长尾把均值拉走了）") print(f""" 观察 Observations： · 每个样本量上逻辑回归都更好，而且差距**不随 n 消失**：GDA 的错误率停在 {eg:.4f}， 比 Bayes 最优（{bayes:.4f}）高 {eg - bayes:.4f}；逻辑回归收敛到 {el:.4f}，只比 Bayes 高 {el - bayes:.4f}。 这是模型错设（misspecification）造成的**偏差**，数据再多也修不好；演示 2 里 GDA 的优势是**方差**上的， 两者的绝对差距随 n 趋于 0。 · 分开看两类：GDA 在正类上犯错多（{g1:.4f} vs 负类 {g0:.4f}）—— 长尾把 μ̂₁ 往远处拉，也把 Σ̂ 沿分类方向撑大， 边界被推向正类那一侧。逻辑回归几乎不理会那些离边界很远、早已分对的点（log loss 在 margin 很大时趋于 0）， 和第 1 讲 §3.1 的“大肿瘤”是同一个机制。 · Ng & Jordan (2002) 的结论正是这两张表的组合：生成式模型更快地接近**它自己的**渐近错误率， 判别式模型的渐近错误率更低（或相同）。数据少、假设大致对 → 生成式；数据多、假设可疑 → 判别式。""") def main(): demo_posterior_is_logistic() demo_assumptions_hold() demo_assumptions_violated() if __name__ == "__main__": main()