"""
=============================================================================
第 3 讲 · 代码 5：模型选择与诊断 —— 交叉验证、选择偏差、学习曲线、"下一步该做什么"
=============================================================================

对应讲义：lectures/lecture03_generalization_regularization.md 第 5、6 节

-----------------------------------------------------------------------------
【这个文件要回答的问题 What this file answers】
-----------------------------------------------------------------------------
    演示 1：用 k 折交叉验证 / LOOCV 选 ridge 的 λ。各种 k 选出的 λ 好不好？CV 的最小值能当作测试误差的估计吗？
    演示 2：选择偏差（selection bias）：在 M 个模型里挑 CV 最好的那个，它的 CV 分数是乐观的 ——
            纯噪声数据上也能"发现"规律。M 越大越严重；nested CV 给出诚实的估计。
    演示 3：学习曲线（learning curves）：高偏差和高方差的模型长什么样？各自该怎么办？
    演示 4：Andrew Ng 的"是优化算法的问题，还是目标函数的问题？"诊断。
    演示 5：Error analysis：流水线里哪个组件最值得改进？（把它的输出换成真值看看）
    演示 6：Ablative analysis：逐个去掉特征组，看哪些真正有用 —— 以及为什么去掉的顺序会影响结论。

-----------------------------------------------------------------------------
【运行方式】
-----------------------------------------------------------------------------
    python3 05_selection_diagnostics.py          # 约 4 秒

依赖本目录的 gen_tools.py。
"""

import numpy as np

from gen_tools import (add_intercept, cross_val_mse, fit_least_squares, fit_logistic_newton, fit_ridge,
                       kfold_indices, logistic_loss, make_sine, poly_features, print_header, ridge_loocv, sigmoid,
                       sine_target, to_unit_interval)

NOISE_VAR = 0.09


def poly_design(x, k):
    """多项式特征 [1, u, …, u^k]，u = 2x − 1（和代码 01 相同）。"""
    return poly_features(to_unit_interval(x), k)


# =============================================================================
# 演示 1：用交叉验证选 λ
# =============================================================================
def fold_mses(Phi, y, folds, lam):
    """
    每一折各自的验证 MSE（而不是合在一起的一个数）。

    【用途】one-standard-error rule 需要知道 CV 估计本身有多大的噪声：各折 MSE 的标准差 / √k 就是它的标准误（SE）。
    【返回】(k,) 数组
    """
    return np.array([np.mean((y[va] - Phi[va] @ fit_ridge(Phi[tr], y[tr], lam)) ** 2) for tr, va in folds])


def cv_curves(x, y, k_deg, lams, rng):
    """
    对一个数据集，算出每个 λ 的 2 折 / 5 折 / 10 折 CV MSE 和 LOOCV MSE，并给出 10 折 + 1-SE 规则选出的下标。

    【参数】x, y：训练数据；k_deg：多项式次数；lams：从小到大的 λ 列表；rng：划分折用的随机数生成器
    【返回】(curves, idx_1se)
        curves  ：dict，方法名 → (len(lams),) 的 CV MSE 数组
        idx_1se ：one-standard-error rule 选中的 λ 的下标 —— 在 10 折 CV 最小值的 1 个标准误以内，选**最大**（最简单）的 λ
    【边界】同一个方法对所有 λ 用**同一组**折（只划分一次）—— 否则不同 λ 之间的比较会混进"划分不同"的噪声。
        LOOCV 用 ridge_loocv 的封闭公式；它假设留一后惩罚的绝对大小 nλ 不变，和 k 折的约定（λ 乘以训练折的样本数）
        差一个 n/(n−1) 的因子，可以忽略。
    """
    Phi = poly_design(x, k_deg)
    out = {}
    for name, k in [("2-fold", 2), ("5-fold", 5)]:
        folds = kfold_indices(len(y), k, rng)
        out[name] = np.array([cross_val_mse(lambda A, b, C, lam=lam: C @ fit_ridge(A, b, lam), Phi, y, folds)
                              for lam in lams])
    folds = kfold_indices(len(y), 10, rng)
    per_fold = np.array([fold_mses(Phi, y, folds, lam) for lam in lams])        # (len(lams), 10)
    out["10-fold"] = per_fold.mean(axis=1)              # n = 30 能被 10 整除，各折一样大，等于合在一起的 MSE
    out["LOOCV"] = np.array([ridge_loocv(Phi, y, lam) for lam in lams])
    i_min = int(np.argmin(out["10-fold"]))
    se = per_fold[i_min].std(ddof=1) / np.sqrt(per_fold.shape[1])
    idx_1se = int(np.max(np.flatnonzero(out["10-fold"] <= out["10-fold"][i_min] + se)))
    return out, idx_1se


def demo_cv():
    print_header("演示 1：用 k 折交叉验证 / LOOCV 选 ridge 的 λ（k = 14 的多项式，n = 30）")
    rng = np.random.default_rng(61)
    k_deg, n = 14, 30
    lams = np.logspace(-7, 1, 17)
    x_te = rng.uniform(0.0, 1.0, size=2000)
    f_te = sine_target(x_te)
    Phi_te = poly_design(x_te, k_deg)

    def test_mse(Phi, y, lam):
        return float(np.mean((Phi_te @ fit_ridge(Phi, y, lam) - f_te) ** 2)) + NOISE_VAR

    # (a) 一个数据集上的曲线
    x, y = make_sine(n, rng)
    Phi = poly_design(x, k_deg)
    curves, i1 = cv_curves(x, y, k_deg, lams, rng)
    test = np.array([test_mse(Phi, y, lam) for lam in lams])
    print("(a) 一个数据集：每一列是一种 CV 估计的 MSE，最后一列是真实的测试误差（2000 个测试点 + 噪声方差）。")
    print(f"    * 标出每列的最小值；10 折 + 1-SE 规则选的是 λ = {lams[i1]:.0e}。\n")
    names = ["2-fold", "5-fold", "10-fold", "LOOCV"]
    print(f"{'lambda':>9}" + "".join(f"{nm:>11}" for nm in names) + f"{'test':>11}")
    for i, lam in enumerate(lams):
        cells = []
        for nm in names:
            star = "*" if i == int(np.argmin(curves[nm])) else " "
            cells.append(f"{curves[nm][i]:>10.4g}{star}")
        star = "*" if i == int(np.argmin(test)) else " "
        print(f"{lam:>9.0e}" + "".join(cells) + f"{test[i]:>10.4g}{star}")
    a_loo, a_test, a_i1 = curves["LOOCV"], test, i1              # 保存 (a) 的结果：下面 (b) 会复用这些变量名
    a_pick = int(np.argmin(curves["10-fold"]))

    # (b) 200 个数据集上的表现
    R = 200
    methods = names + ["10-fold+1SE"]
    chosen = {nm: [] for nm in methods}
    cvmin = {nm: [] for nm in names}
    loglam = {nm: [] for nm in methods}
    oracle = []
    for _ in range(R):
        x, y = make_sine(n, rng)
        Phi = poly_design(x, k_deg)
        curves, i1 = cv_curves(x, y, k_deg, lams, rng)
        test = np.array([test_mse(Phi, y, lam) for lam in lams])
        oracle.append(test.min())
        for nm in methods:
            i = i1 if nm == "10-fold+1SE" else int(np.argmin(curves[nm]))
            chosen[nm].append(test[i])
            loglam[nm].append(np.log10(lams[i]))
            if nm in cvmin:
                cvmin[nm].append(curves[nm][i] - test[i])
    oracle = np.array(oracle)
    print(f"\n(b) {R} 个数据集。oracle = 每个数据集上按测试误差挑 λ（不可实现的上限）：中位数 {np.median(oracle):.4f}，平均 {np.mean(oracle):.4f}。")
    print("    test@chosen = 选出的 λ 的真实测试误差；P(>2x oracle) = 选出的模型比 oracle 差一倍以上的比例（\"翻车\"）。\n")
    print(f"{'method':>12}{'median test@chosen':>20}{'mean test@chosen':>18}{'P(>2x oracle)':>15}"
          f"{'median(minCV-test)':>20}{'median log10 lam':>18}")
    summ = {}
    for nm in methods:
        c = np.array(chosen[nm])
        gap = f"{np.median(cvmin[nm]):>+20.4f}" if nm in cvmin else f"{'-':>20}"
        summ[nm] = (float(np.median(c)), float(np.mean(c)), float(np.mean(c > 2 * oracle)), float(np.median(loglam[nm])))
        print(f"{nm:>12}{summ[nm][0]:>20.4f}{summ[nm][1]:>18.4g}{summ[nm][2]:>15.3f}{gap}{summ[nm][3]:>18.1f}")
    print(f"""
观察 Observations：
  · (a) 四条 CV 曲线的谷底都在 λ = {lams[a_pick]:.0e}，形状和测试曲线一致：λ 太小过拟合、太大欠拟合。但单个数据集上的测试曲线本身
    在 1e−4 ～ 3e−3 之间起伏不平，CV 选的 λ 的测试误差是 {a_test[a_pick]:.3f}，1-SE 规则选的 λ = {lams[a_i1]:.0e} 是 {a_test[a_i1]:.3f}，最好的是 {a_test.min():.3f}。
  · 看 λ 很小的几行：LOOCV 在 λ = 3e−7 处只有 {a_loo[1]:.3f}，真实测试误差却是 {a_test[1]:.3g}。CV 只在**训练点所在的位置**检查预测；
    这里的灾难发生在 [0, 1] 两端没有训练点的空白处（多项式外推），留一法几乎"看不到"那里。
  · (b) 典型情况下（中位数）各种方法都选得不错：{summ['2-fold'][0]:.3f} / {summ['5-fold'][0]:.3f} / {summ['10-fold'][0]:.3f} / {summ['LOOCV'][0]:.3f}，oracle 是 {np.median(oracle):.3f}。
    但平均值差别很大（{summ['2-fold'][1]:.3g} / {summ['5-fold'][1]:.3g} / {summ['10-fold'][1]:.3g} / {summ['LOOCV'][1]:.3g}）：少数数据集上 CV 选了太小的 λ，测试误差爆炸。
  · 2 折每次只用一半数据训练，CV 曲线偏悲观、偏向更大的 λ（中位数 10^{summ['2-fold'][3]:.1f} vs LOOCV 的 10^{summ['LOOCV'][3]:.1f}）。
    它"翻车"的频率并不低（P(>2x oracle) = {summ['2-fold'][2]:.3f}），但翻得不重 —— 很少选到会爆炸的极小 λ，所以平均值最好之一。
    这和教科书的说法一致：k 越大，CV 的偏差越小，但估计本身的方差可能更大（LOOCV 的 n 个训练集几乎一模一样）。
  · one-standard-error rule（在最小值 1 个标准误以内选最简单的模型，ESL §7.10）：中位数 {summ['10-fold+1SE'][0]:.3f}（典型情况略差一点），
    翻车率 {summ['10-fold+1SE'][2]:.3f}，平均 {summ['10-fold+1SE'][1]:.3g}（所有方法里最好）。它承认"CV 曲线本身有噪声"，在差不多好的模型里偏向更稳的那个。
  · "minCV − test"（中位数）：2 折略悲观（{np.median(cvmin['2-fold']):+.4f}），5 折以上是乐观的（{np.median(cvmin['5-fold']):+.4f} ～ {np.median(cvmin['10-fold']):+.4f}）——
    被选中的模型的 CV 分数比它的真实测试误差好看，因为"取最小值"这个动作本身引入了偏差（演示 2）。
  · LOOCV 对 ridge 有封闭公式（gen_tools.ridge_loocv），和训练一次一样便宜；一般模型的 LOOCV 要训练 n 次。""")


# =============================================================================
# 演示 2：选择偏差与 nested CV
# =============================================================================
def single_feature_cv(X, y, folds):
    """
    同时对 M + 1 个模型算 CV MSE：模型 0 = 只有截距；模型 j（1..M）= 截距 + 第 j 个特征的简单线性回归。

    【思路】一元线性回归有封闭解：slope = Σ(x − x̄)(y − ȳ) / Σ(x − x̄)²，intercept = ȳ − slope·x̄。
        用 broadcasting 一次对 M 个特征都算，不写循环。
    【参数】X：(n, M)；y：(n,)；folds：kfold_indices 的结果
    【返回】(M + 1,) 每个模型的 CV MSE
    """
    sq = np.zeros(X.shape[1] + 1)
    for tr, va in folds:
        xm, ym = X[tr].mean(axis=0), y[tr].mean()
        xc = X[tr] - xm
        slope = xc.T @ (y[tr] - ym) / np.sum(xc ** 2, axis=0)            # (M,)
        pred = ym + (X[va] - xm) * slope                                 # (n_va, M)
        sq[0] += np.sum((y[va] - ym) ** 2)
        sq[1:] += np.sum((y[va][:, None] - pred) ** 2, axis=0)
    return sq / len(y)


def fit_single(X, y, j):
    """拟合模型 j（0 = 只有截距），返回 (intercept, slope)。"""
    if j == 0:
        return y.mean(), 0.0
    x = X[:, j - 1]
    slope = float((x - x.mean()) @ (y - y.mean()) / np.sum((x - x.mean()) ** 2))
    return y.mean() - slope * x.mean(), slope


def demo_selection_bias():
    print_header("演示 2：在 M 个模型里挑 CV 最好的 —— 选择偏差（selection bias）与 nested CV")
    rng = np.random.default_rng(62)
    n, R = 50, 200
    print(f"纯噪声数据：n = {n}，y ~ N(0, 1)，M 个特征 x_j ~ N(0, 1) 全部和 y 无关。候选模型 = {{只有截距}} ∪ {{截距 + 一个特征 x_j}}。")
    print("用 5 折 CV 挑最好的模型，再用全部数据重新拟合它。它的真实测试 MSE 可以精确算出：1 + θ₀² + θ₁²（≥ 1）。")
    print(f"nested CV：外层 5 折，每个外层训练集上**重新做一遍**内层 5 折选择 + 拟合，在外层验证折上评估。每行 {R} 个数据集的平均。\n")
    print(f"{'M':>6}{'CV score of winner':>20}{'true MSE of winner':>20}{'optimism':>10}{'nested CV':>11}"
          f"{'P(pick intercept)':>19}")
    rows = {}
    for M in [1, 10, 100, 1000]:
        cv_best, true_mse, nested, pick0 = [], [], [], []
        for _ in range(R):
            X = rng.normal(size=(n, M))
            y = rng.normal(size=n)
            scores = single_feature_cv(X, y, kfold_indices(n, 5, rng))
            j = int(np.argmin(scores))
            b0, b1 = fit_single(X, y, j)
            cv_best.append(scores[j])
            true_mse.append(1.0 + b0 ** 2 + b1 ** 2)
            pick0.append(j == 0)
            # nested CV：整个"选模型 + 拟合"的流程被当作一个算法，放进外层 CV 里评估
            sq = 0.0
            for tr, va in kfold_indices(n, 5, rng):
                inner = single_feature_cv(X[tr], y[tr], kfold_indices(len(tr), 5, rng))
                jj = int(np.argmin(inner))
                c0, c1 = fit_single(X[tr], y[tr], jj)
                pred = c0 + (c1 * X[va, jj - 1] if jj > 0 else 0.0)
                sq += float(np.sum((y[va] - pred) ** 2))
            nested.append(sq / n)
        rows[M] = (float(np.mean(cv_best)), float(np.mean(true_mse)), float(np.mean(nested)), float(np.mean(pick0)))
        r = rows[M]
        print(f"{M:>6}{r[0]:>20.4f}{r[1]:>20.4f}{r[0] - r[1]:>+10.4f}{r[2]:>11.4f}{r[3]:>19.3f}")
    print(f"""
观察 Observations：
  · 所有特征都是噪声，最好的模型是"只有截距"，真实 MSE ≈ 1。可是 M = 1000 时，被选中的模型的 CV 分数平均只有 {rows[1000][0]:.3f}，
    看起来"解释了 {1 - rows[1000][0]:.0%} 的方差"；它的真实 MSE 却是 {rows[1000][1]:.3f}（比只用截距还差）。
  · 乐观偏差随 M 增大（M = 10：{rows[10][0] - rows[10][1]:+.3f}；M = 1000：{rows[1000][0] - rows[1000][1]:+.3f}）：M 个带噪声的 CV 分数里取最小值，
    最小值大约比均值低 σ_CV·√(2 log M)。这正是 §3 的 union bound 的另一面：验证集是一个"有限假设类 + m 个样本"的学习问题，
    log M 越大，越容易被某个假设的好运气骗到。
  · 单独的 CV 分数对**选出来的模型**不再是无偏估计：数据已经被"用来选择"了一次。nested CV 把"选择"也放进外层循环，
    给出诚实的估计（M = 1000：{rows[1000][2]:.3f}，略高于真实值 {rows[1000][1]:.3f}，因为外层每次只用 4/5 的数据）。
  · 实践教训：Kaggle 的 public leaderboard 过拟合、论文里"在测试集上调参"、在很多特征里挑和结果相关的那几个
    （p-hacking），都是同一件事。测试集只能看一次；看了多少次，就要为多少次付出 log M 的代价。""")


# =============================================================================
# 演示 3：学习曲线
# =============================================================================
def demo_learning_curves():
    print_header("演示 3：学习曲线 —— 训练 / 测试误差随训练集大小 n 的变化")
    rng = np.random.default_rng(63)
    x_te = rng.uniform(0.0, 1.0, size=2000)
    f_te = sine_target(x_te)
    models = {"k=1 (line)": (1, 0.0), "k=12, no reg": (12, 0.0), "k=12, ridge 1e-3": (12, 1e-3)}
    ns = [20, 40, 80, 160, 320, 640]
    R = 200
    print("数据：y = sin(2πx) + N(0, 0.3²)。噪声下限 0.09。假设目标是测试 MSE ≤ 0.12（\"desired performance\"）。")
    print(f"每个 n 重复 {R} 个训练集，表中是**中位数**（高次多项式的测试误差是重尾的，见代码 01 演示 1）。\n")
    head = "".join(f"{nm + ' tr/te':>26}" for nm in models)
    print(f"{'n':>5}{head}")
    res = {nm: {} for nm in models}
    for n in ns:
        cells = []
        for nm, (k, lam) in models.items():
            tr, te = [], []
            Phi_te = poly_design(x_te, k)
            for _ in range(R):
                x, y = make_sine(n, rng)
                Phi = poly_design(x, k)
                th = fit_ridge(Phi, y, lam)
                tr.append(np.mean((Phi @ th - y) ** 2))
                te.append(np.mean((Phi_te @ th - f_te) ** 2) + NOISE_VAR)
            res[nm][n] = (float(np.median(tr)), float(np.median(te)))
            cells.append(f"{res[nm][n][0]:.4f} / {res[nm][n][1]:.4f}")
        print(f"{n:>5}" + "".join(f"{c:>26}" for c in cells))
    lin, hv, rg = res["k=1 (line)"], res["k=12, no reg"], res["k=12, ridge 1e-3"]
    print(f"""
观察 Observations：
  · 直线（k = 1）= 高偏差（high bias）：训练和测试误差很快就靠在一起（n = 640：{lin[640][0]:.3f} / {lin[640][1]:.3f}），
    但都停在 0.3 附近，远高于目标 0.12。**更多数据没用**（n 从 80 到 640，测试误差 {lin[80][1]:.3f} → {lin[640][1]:.3f}）。
    该做的是：更多 / 更好的特征、更大的模型、更小的正则化。
  · 12 次多项式不加正则 = 高方差（high variance）：小 n 时训练误差远低于测试误差（n = 20：{hv[20][0]:.3f} vs {hv[20][1]:.3f}），
    差距随 n 缩小（n = 640：{hv[640][0]:.4f} vs {hv[640][1]:.4f}），两条曲线还在靠近 —— 外推一下，**更多数据有用**。
    该做的是：更多数据、更少的特征、更强的正则化。
  · 加了 ridge（λ = 1e−3）的同一个模型：小 n 时测试误差低得多（n = 20：{rg[20][1]:.3f} vs {hv[20][1]:.3f}），
    n 大时两者几乎一样（数据多了，正则化就不那么重要）。正则化相当于"用先验换数据"。
  · 读学习曲线的要点：看**最终的差距**（gap，诊断 variance）和**训练误差本身离目标多远**（诊断 bias）。
    训练误差都达不到目标时，再多数据也救不了 —— 训练误差通常只会随 n 上升。""")


# =============================================================================
# 演示 4：是优化算法的问题，还是目标函数的问题？
# =============================================================================
def weighted_cost(theta, X, y, fp_cost=10.0, fn_cost=1.0):
    """
    我们真正关心的指标：平均错误代价。把非垃圾邮件（y = 0）判成垃圾（误删正常邮件）代价 10，漏掉垃圾邮件代价 1。

    【返回】(平均代价, 误删率 FP rate, 漏判率 FN rate)
    """
    pred = (X @ theta >= 0).astype(int)
    fp = np.mean(pred[y == 0] == 1)
    fn = np.mean(pred[y == 1] == 0)
    cost = (fp_cost * np.sum((pred == 1) & (y == 0)) + fn_cost * np.sum((pred == 0) & (y == 1))) / len(y)
    return float(cost), float(fp), float(fn)


def make_spam(n, rng):
    """
    玩具"垃圾邮件"数据：两个特征，y = 1 是垃圾邮件。第二个特征的尺度是第一个的 100 倍（模拟未标准化的计数特征）。
    两类有重叠，不可分。
    """
    n1 = n // 2
    X1 = rng.normal([1.0, 1.0], 1.0, size=(n1, 2))
    X0 = rng.normal([-1.0, -1.0], 1.0, size=(n - n1, 2))
    X = np.vstack([X1, X0]) * np.array([1.0, 100.0])
    y = np.concatenate([np.ones(n1, dtype=int), np.zeros(n - n1, dtype=int)])
    return X, y


def gd_logistic(X, y, steps):
    """逻辑回归的 batch GD（固定步长 1/L，L = λ_max(XᵀX/n)/4），只跑 steps 步 —— 故意不跑到收敛。"""
    lr = 1.0 / (np.linalg.eigvalsh(X.T @ X / len(y))[-1] / 4.0)
    theta = np.zeros(X.shape[1])
    for _ in range(steps):
        theta -= lr * X.T @ (sigmoid(X @ theta) - y) / len(y)
    return theta


def demo_optimizer_vs_objective():
    print_header("演示 4：Ng 的诊断 —— 是优化算法（optimizer）的问题，还是目标函数（objective）的问题？")
    rng = np.random.default_rng(64)
    X, y = make_spam(2000, rng)
    X_te, y_te = make_spam(20000, rng)
    A, A_te = add_intercept(X), add_intercept(X_te)
    # 方法 B：加权逻辑回归 —— 把每个 y = 0 的样本复制 10 份（等价于给它们 10 倍的权重），用 Newton 跑到收敛。
    #   它扮演 Ng 例子里的"SVM"：另一个方法，在我们真正关心的指标上表现更好。
    w_idx = np.concatenate([np.flatnonzero(y == 1)] + [np.flatnonzero(y == 0)] * 10)
    th_B, _ = fit_logistic_newton(A[w_idx], y[w_idx])
    J_B, a_B = logistic_loss(th_B, A, y), weighted_cost(th_B, A_te, y_te)
    print("问题：误删一封正常邮件的代价是 10，漏掉一封垃圾邮件的代价是 1。我们关心 a(θ) = 测试集上的平均代价（越小越好）。")
    print("方法 A = 普通逻辑回归，目标 J(θ) = 不加权的平均 log loss，用 GD 训练（特征没标准化，第二个特征的尺度是第一个的 100 倍）。")
    print("方法 B = 把正常邮件的权重设为 10 的逻辑回归（Newton 收敛）。它在 a 上比 A 好 —— 问题是：A 差在哪？")
    print("诊断：拿 A 自己的目标 J 去评价 B 的参数（J 是损失，越小越好）。")
    print("  J(θ_B) < J(θ_A)：B 在 A 自己的目标上都赢了 A ⇒ A 的优化器没把 J 优化好 → optimizer 的问题")
    print("  J(θ_B) ≥ J(θ_A)：J 认为 A 更好，a 却认为 B 更好 ⇒ J 和 a 不一致 → objective 的问题\n")
    print(f"{'parameters':<24}{'J (A objective)':>16}{'a = cost':>10}{'FP rate':>9}{'FN rate':>9}   diagnosis")
    print(f"{'B: weighted LR':<24}{J_B:>16.4f}{a_B[0]:>10.4f}{a_B[1]:>9.3f}{a_B[2]:>9.3f}   (reference)")
    rows = {}
    for steps in [5, 20, 50, 200, 1000, 5000]:
        th = gd_logistic(A, y, steps)
        rows[steps] = (logistic_loss(th, A, y), *weighted_cost(th, A_te, y_te))
    th_N, _ = fit_logistic_newton(A, y)
    rows["Newton"] = (logistic_loss(th_N, A, y), *weighted_cost(th_N, A_te, y_te))
    for key, r in rows.items():
        name = f"A: {key} GD steps" if key != "Newton" else "A: Newton (converged)"
        diag = "optimizer" if J_B < r[0] else "objective"
        print(f"{name:<24}{r[0]:>16.4f}{r[1]:>10.4f}{r[2]:>9.3f}{r[3]:>9.3f}   {diag}")
    opt_rows = [k for k, r in rows.items() if J_B < r[0]]
    obj_rows = [k for k, r in rows.items() if J_B >= r[0]]
    print(f"""
观察 Observations：
  · 只跑 {', '.join(str(k) for k in opt_rows)} 步时 J(θ_A) > J(θ_B) = {J_B:.4f}：B 的参数在 A 自己的目标上都比 A 好，
    A 的优化器显然没完成任务 —— 这时改目标函数、加特征都是白费力气，先修优化器（多跑几步、标准化特征、换 Newton）。
  · 从 {obj_rows[0]} 步开始 J(θ_A) ≤ J(θ_B)，诊断变成"objective"：J 认为 A 比 B 好，而我们关心的代价 a 认为 B 好得多
    （Newton 收敛后 a = {rows['Newton'][1]:.4f} vs B 的 {a_B[0]:.4f}）。不加权的 log loss 认为两种错误一样糟，而我们不这么认为 ——
    该改的是目标（加权损失，或者调整判决阈值），不是优化器。
  · 注意两个结论**不对称**：50 步时 GD 还远没收敛（J = {rows[50][0]:.4f}，最优是 {rows['Newton'][0]:.4f}），诊断却已经说"objective"。
    "J(θ_B) ≥ J(θ_A)"只证明了 J 和 a 不一致（目标确实有问题），并不证明优化器没问题 —— 两个问题可以同时存在。
    反过来，"J(θ_B) < J(θ_A)"则确凿地说明优化器有问题。所以修好一个之后要再诊断一次。
  · 这个诊断的价值在于它把两种完全不同的"下一步"分开。同样的思路在强化学习里也用
    （第 6 讲：是模拟器、RL 算法还是奖励函数的问题？）。""")


# =============================================================================
# 演示 5：Error analysis —— 流水线里哪个组件最值得改进？
# =============================================================================
def make_pipeline_data(n, rng, noise_a=1.0, noise_b=0.3):
    """
    一个两阶段的玩具流水线：真实的潜变量 s1、s2 ~ N(0, 1)，标签 y ~ Bernoulli(σ(3 s1 + 1.5 s2))。
    组件 A 输出 s1 的带噪估计 ŝ1 = s1 + N(0, noise_a²)（一个不太好的组件）；
    组件 B 输出 s2 的带噪估计 ŝ2 = s2 + N(0, noise_b²)（一个不错的组件）。
    最后一层逻辑回归只能看到 [ŝ1, ŝ2]。

    【返回】dict：s1, s2, s1_hat, s2_hat, y
    """
    s1, s2 = rng.normal(size=n), rng.normal(size=n)
    y = (rng.uniform(size=n) < sigmoid(3.0 * s1 + 1.5 * s2)).astype(int)
    return {"s1": s1, "s2": s2, "s1_hat": s1 + rng.normal(0.0, noise_a, size=n),
            "s2_hat": s2 + rng.normal(0.0, noise_b, size=n), "y": y}


def pipeline_accuracy(train, test, use_true_a, use_true_b):
    """把组件 A / B 的输出换成真值（"perfect component"），重新训练最后的逻辑回归，返回测试准确率。"""
    def feats(d):
        return add_intercept(np.column_stack([d["s1"] if use_true_a else d["s1_hat"],
                                              d["s2"] if use_true_b else d["s2_hat"]]))
    th, _ = fit_logistic_newton(feats(train), train["y"], l2=1e-4)
    return float(np.mean((feats(test) @ th >= 0) == test["y"]))


def demo_error_analysis():
    print_header("演示 5：Error analysis —— 依次把每个组件换成\"完美\"的，看整体准确率涨多少")
    rng = np.random.default_rng(65)
    train, test = make_pipeline_data(2000, rng), make_pipeline_data(50000, rng)
    print("流水线：输入 → 组件 A（估计 s1，噪声标准差 1.0）→ 组件 B（估计 s2，噪声标准差 0.3）→ 逻辑回归。")
    print("真实标签 y ~ Bernoulli(σ(3 s1 + 1.5 s2))。准确率在 50000 个测试样本上算。\n")
    configs = [("overall system", False, False), ("A perfect", True, False), ("B perfect", False, True),
               ("A + B perfect", True, True)]
    acc = {name: pipeline_accuracy(train, test, a, b) for name, a, b in configs}
    print(f"{'component(s) replaced by ground truth':<40}{'test accuracy':>14}{'gain vs overall':>17}")
    for name, _, _ in configs:
        print(f"{name:<40}{acc[name]:>14.4f}{acc[name] - acc['overall system']:>+17.4f}")
    print(f"""
观察 Observations：
  · 把 A 换成完美组件，整体准确率涨 {acc['A perfect'] - acc['overall system']:+.4f}；把 B 换成完美组件只涨 {acc['B perfect'] - acc['overall system']:+.4f}。
    所以下一步应该花时间改进 A，即使 B 看起来"也不完美"。
  · 两个都完美时是 {acc['A + B perfect']:.4f}：这是这条流水线的天花板（剩下的错误来自标签本身的随机性）。
    Ng 的人脸识别例子用的就是这张表：沿着流水线逐个"插入真值"，看每个组件能贡献多少上限。
  · 为什么不直接看"哪个组件的误差最大"？因为组件的误差对最终任务的影响不一样大：
    这里 s1 的系数是 s2 的两倍，A 的噪声又大，所以 A 的误差被放大了。error analysis 直接量的是**对最终指标的影响**。""")


# =============================================================================
# 演示 6：Ablative analysis —— 去掉哪个特征组损失最大？
# =============================================================================
def demo_ablation():
    print_header("演示 6：Ablative analysis —— 逐个去掉特征组，以及为什么顺序会影响结论")
    rng = np.random.default_rng(66)
    n, n_te = 2000, 50000

    def make(n):
        s1, s2 = rng.normal(size=n), rng.normal(size=n)
        y = (rng.uniform(size=n) < sigmoid(3.0 * s1 + 1.5 * s2)).astype(int)
        groups = {
            "A  (s1 sensor 1)": (s1 + rng.normal(0.0, 0.7, size=n))[:, None],
            "A' (s1 sensor 2)": (s1 + rng.normal(0.0, 0.7, size=n))[:, None],   # 和 A 测的是同一个东西
            "B  (s2 sensor)": (s2 + rng.normal(0.0, 0.3, size=n))[:, None],
            "junk (10 noise feats)": rng.normal(size=(n, 10)),
        }
        return groups, y

    G_tr, y_tr = make(n)
    G_te, y_te = make(n_te)

    def acc(names):
        """只用 names 里的特征组训练逻辑回归，返回测试准确率（names 为空时只有截距）。"""
        if names:
            F_tr = add_intercept(np.hstack([G_tr[g] for g in names]))
            F_te = add_intercept(np.hstack([G_te[g] for g in names]))
        else:
            F_tr, F_te = np.ones((n, 1)), np.ones((n_te, 1))
        th, _ = fit_logistic_newton(F_tr, y_tr, l2=1e-4)
        return float(np.mean((F_te @ th >= 0) == y_te))

    allg = list(G_tr)
    full = acc(allg)
    print("4 个特征组：A 和 A' 是同一个潜变量 s1 的两个独立的带噪测量（冗余）；B 测 s2；junk 是 10 个纯噪声特征。")
    print(f"全部特征的测试准确率 = {full:.4f}（{n} 个训练样本，{n_te} 个测试样本）。\n")
    print("(a) 一次只去掉一个组（leave-one-group-out）：")
    print(f"{'removed':<26}{'accuracy':>10}{'drop':>9}")
    single = {}
    for g in allg:
        single[g] = acc([h for h in allg if h != g])
        print(f"{g:<26}{single[g]:>10.4f}{single[g] - full:>+9.4f}")
    both = acc([h for h in allg if h not in ("A  (s1 sensor 1)", "A' (s1 sensor 2)")])
    label = "A + A' together"                  # 【Python 小知识】含单引号的文字先放进变量，避免在 f-string 里嵌套引号
    print(f"{label:<26}{both:>10.4f}{both - full:>+9.4f}")
    print("\n(b) Ng 式的累积消融（从完整系统开始，一次去掉一个组，按下面的顺序）：")
    order = ["junk (10 noise feats)", "A' (s1 sensor 2)", "B  (s2 sensor)", "A  (s1 sensor 1)"]
    cur = list(allg)
    print(f"{'system':<40}{'accuracy':>10}")
    print(f"{'full':<40}{full:>10.4f}")
    cum = {}
    for g in order:
        cur = [h for h in cur if h != g]
        cum[g] = acc(cur)
        print(f"{'  − ' + g:<40}{cum[g]:>10.4f}")
    a_drop = single["A  (s1 sensor 1)"] - full
    a2_drop = single["A' (s1 sensor 2)"] - full
    print(f"""
观察 Observations：
  · (a) 单独去掉 A 只损失 {a_drop:+.4f}，单独去掉 A' 只损失 {a2_drop:+.4f}，都比去掉 B（{single['B  (s2 sensor)'] - full:+.4f}）小 ——
    看起来 B 最重要。可 s1 在真实模型里的系数是 s2 的两倍：同时去掉 A 和 A' 损失 {both - full:+.4f}。
    A 和 A' 是冗余的，互相顶替。leave-one-out 式的消融会系统性地低估冗余组件的价值。
  · 去掉 junk 的变化是 {single['junk (10 noise feats)'] - full:+.4f}：10 个噪声特征在 2000 个样本上几乎无害（估计方差很小）；
    样本少、噪声特征多时它们才会明显拖后腿。
  · (b) 累积消融的结论取决于顺序：按这个顺序，A' 被去掉时只损失 {cum["A' (s1 sensor 2)"] - cum['junk (10 noise feats)']:+.4f}（因为 A 还在），
    最后去掉 A 时损失 {cum['A  (s1 sensor 1)'] - cum['B  (s2 sensor)']:+.4f}（因为 A' 已经不在了）。如果把 A 和 A' 的顺序对调，两者的"功劳"也会大致对调（A 和 A' 的地位是对称的）。
  · 实践建议：写论文 / 做汇报时的消融表，要说明顺序；遇到可能冗余的组件，同时报告"成对去掉"的结果。
    Error analysis（演示 5）问"把哪个组件做到完美收益最大"，ablative analysis 问"现有的哪个组件在起作用"，两者互补。""")


def main():
    demo_cv()
    demo_selection_bias()
    demo_learning_curves()
    demo_optimizer_vs_objective()
    demo_error_analysis()
    demo_ablation()


if __name__ == "__main__":
    main()