统计推断、实验设计与回归

拟合系数或基准分数本身不建立不确定性、显著性或因果。推断问：在抽样和分配条件下，一个实验程序支持哪些结论。本课构造置信覆盖、有效原假设尾概率、配对比较、重抽样与回归目标，并保护最终评估免受模型选择影响。结果以前就写目标和协议，因为解释依赖这些选择。

检索准备：模块 12的最小二乘几何、模块 19的优化、模块 25的似然估计。全课应用模块 24 的抽样保证。实验 1、3 为标准库，实验 2 使用固定 NumPy 环境，全部 CPU、局部固定种子。

置信区间是固定参数 θ 的数据依赖集合 C(X)。1−α 水平保证指所声明模型内每个 θ 都有 P_θ(θ∈C(X))≥1−α，或者明确限定的渐近版本。抽样以前区间随机；数据实现后端点固定，没有先验时不由此给 θ 后验概率。一次区间或者涵盖真值、或者没有，概率保证属于反复使用构造程序。

IID 高斯、已知 σ²>0 时，x_bar 精确服从均值 μ、方差 σ²/n 的高斯。取标准高斯分位 z，使 P(|Z|≤z)=1−α，则 P_μ(|x_bar−μ|≤zσ/√n)=1−α。反解此事件得 C=[x_bar−zσ/√n,x_bar+zσ/√n]，从而证明精确覆盖。α=.05 时 z≈1.959964。任意观测中将已知方差换为估计量再用同常数，是另一个通常近似的程序。

IID [0,1] 观测共同均值 μ，Hoeffding 给半径 r=√[log(2/α)/(2n)] 的偏差概率≤α。因此闭区间 [x_bar−r,x_bar+r] 与 [0,1] 相交后，覆盖至少 1−α。严格/非严格事件差异只会加强闭区间的保守保证；裁到已知参数空间不删除真 μ。它在有界独立模型内无需指定具体分布，不适用于依赖行或无界回归残差。

常用 Bernoulli Wald 区间为 p_hat±z√[p_hat(1−p_hat)/n]，有限覆盖可能很差，尤其接近零一。裁剪不能修复零半径导致的漏覆盖。无成功时区间退化到零，而非零罕见率可高概率产生此数据。还有其他二项区间方法；这里直接比较透明的 Wald 近似与有效保守有界构造，不默称所有显示区间都精确。

覆盖探索器枚举合成 Bernoulli 的每个成功数，对每个计数构造区间，再将涵盖指定真 p 的样本概率相加。它是概率加权覆盖，而不是把各计数视为等可能。竖线标合成真值，水平区间颜色表示是否包含。仅因这是模型实验才能知道真率；现实未知 p 时，不能由一次数据算实际覆盖。

高斯 σ 未知时，IID 高斯条件下可用 n−1 自由度的 Student t 枢轴给精确区间。这是依赖均值与样本方差关系的命名分布结果，不能由任意经验离散度替换 σ 自动推出。一般有限方差 CLT 区间为渐近，可能在小样本失败。报告所用版本、独立单位、标准误和有意义单位的宽度；保守保证也不说每个实现区间都实用。

原假设 H₀ 指定一组数据规律或参数值。检验统计量 T 按预定方向概括证据。若每个原假设规律下拒绝概率≤α，检验大小至多 α。第一类错误是真原假设被拒绝，第二类错误是在指定备择下未拒绝；功效是该备择下拒绝概率。它们描述程序，不是一次结果“为错误”的概率。

简单完全指定原假设下，p 值是至少与所观测统计量同样极端的原假设尾概率，极端规则必须说明。有效 p 值满足 P₀(p_value≤u)≤u。连续原假设的适当 CDF 变换可给均匀分布，离散原子不能任意分割，常有保守大小。复合原假设需要全参数集合有效的构造，例如适当最坏尾部，而非随意选一个方便参数。

IID Bernoulli、固定 n 检验 p₀=.5，可选 T=|K−n/2|，精确双侧 p 值加总所有 T 不小于观测值的二项质量，含并列。离散双侧还有其他约定，数值可能不同，应命名规则。备择方向事先指定时可用单侧；看到哪边效应大才选方向会改变程序与错误控制。

p 值不是 P(H₀|数据)，不是结果“由偶然造成”的概率，也不是精确数据集概率。后验赔率需要先验与备择似然。拒绝过简原假设也不说明哪个假设失败；依赖、混杂、选择和错误方差模型都可能使尾概率无效。应分开写统计原假设与科学目标。

未拒绝不证明相等。小或噪声大的实验可能无法检测重要效应。等效问题需要预指定有意义容差与为其设计的程序。除了检验，应报告估计、区间及实际尺度。一个 p 值略低 .05、另一个略高 .05，其显著标签不同并不自动构成已经检验的差异。

功效规划说明为何样本量属于协议。已知 σ 高斯单侧检验拒绝 x_bar>μ₀+z_(1−α)σ/√n，备择 μ=μ₀+Δ 下功效为 1−Φ(z_(1−α)−Δ√n/σ)。Δ=0 时大小 α，正 Δ 增功效，n 四倍使信号/标准误项两倍。它依赖高斯模型。看完结果选择的效应不是前瞻功效设定，事后功效也不替代规划。

有限简单原假设的有效性可直接证明。将不同统计量值按由极端到不极端排列，依序累积其原假设质量。各值的 p 值就是含同值并列的累计质量。对任意 u，p 值≤u 的事件是一个初始值集合，按构造其总质量≤u；若没有够小的累计量，事件为空。因此 P₀(p_value≤u)≤u。质量跳过阈值说明离散保守性。观察后再选极端顺序改变构造，需要新的有效性分析。

测量尺度也属于问题。.01 秒、.01 概率单位或 .01 个标准差有不同实际意义。按应用预指定材料效应或容差，再规划精度与功效。检验可检测理想原假设的微小偏离，却几乎不改变实际决策。

应指定被抽样或分配的单位。完全复制模型中，一人十行是一名独立受试者。一个提示的不同变体可能共享难度；按行拆分或重抽样会将近似相同信息放到不同组并低估不确定性。若主体、时间块或簇定义独立模型与目标，拆分和重抽样要保留它们。

随机分配按照已知设计使标签独立于分配前特征，以支持内部比较；随机抽样支持总体代表性。随机实验可能内部可靠却不具代表性，代表性观察样本仍可能混杂。随机化也不自动修复流失、单位干扰、不依从或测量错误。目标可为分配效应、额外假设下的治疗效应，或预测差异，应说清。

潜在结果 Y_i(1)、Y_i(0) 表示同一单位两种分配下的结果，只能观察实际分配的一项，所以个体差异不可直接得到。随机化在适当分配概率、无干扰与稳定处理含义下可支持平均比较。观察回归系数不会仅因小 p 值获得因果解释。调整已测变量也需要因果条件，不能保证消除未测混杂。

配对是在同一抽样单位上比较 A_i 与 B_i。独立单位的 D_i=A_i−B_i 均值估计方差为 Var(D_i)/n。Var(D)=VarA+VarB−2Cov，正的共同难度协方差可提升精度。不配对公式丢掉协方差；任意行顺序拼对并非有效设计。对须有意义，独立条件针对不同对。

观察结果前报告分配、排除、缺失处理、主要结果、方向、样本量与分析。多个训练种子描述固定数据下优化器波动，不替代独立总体受试者；一个训练种子下大量测试例也不测跨种子波动。基准可能需两层分析，但不能把每个种子—样本对全当联合独立。

非参数自助法（bootstrap）从经验分布有放回重抽样，并重新计算完整估计器，以近似抽样波动。IID 均值重样本仍为 n 项、索引独立抽取。重复是条件数值实验，不是新总体数据。增加重复数减少近似模拟噪声，不会补足未见罕见事件、修复偏样本或错误独立单位。

百分位区间取自助分布的 .025、.975 分位，在合适正则条件和足够数据时有用，但不自动精确。偏差、边界、离散与非光滑估计器可使覆盖差。单观测经验分布产生全部相同均值，总体却仍可能波动。配对数据重抽整个对或差值；簇应按独立簇及目标重抽，大小不等时还要分析权重。时间序列需另有依赖方法。

置换检验通过原假设下保持数据律不变的变换，或随机设计允许的分配，生成原假设参考分布。两组 IID 且完整分布相等时，给定混合观测的标签可交换。仅均值相等，在方差或形状不同下不保证该性质。随机实验中，所有单位无效应的尖锐原假设可按已知分配设计检验；变换须保留组大小、配对及分配限制。

实验 1 也独立重抽两组给有限百分位区间。相同观测回答不同问题：近似估计量波动与可交换原假设尾。围绕观测效应的自助分布不能直接当置换零效应参考。读图以前说明目标，不能只选得到喜欢显著标签的方法。

Monte Carlo 置换常使用包含原统计量的 (b+1)/(B+1)，其中 b 是至少同样极端的随机零变换数，B 是抽取数。这在正确可交换随机化方案下避免有限抽样给出无依据零 p 值；有效性仍依赖变换抽取方式及观察分配是否加入可交换集合。完整枚举则用极端数/所有允许数。

构成估计器的预处理若属于重抽程序，每次 bootstrap 都应重新拟合。保持拟合变换固定，则目标是条件程序，遗漏其拟合波动。固定最终模型上的独立主体重抽可度量条件于该模型的评估波动，不自动纳入重新训练或超参数选择。应说明什么固定，与模块 25 的条件区别一致。

固定 n×d 设计 X 的线性模型 y=Xβ+ε，在独立同方差高斯、已知 σ² 时，负 logL 为 ||y−Xβ||²/(2σ²) 加常数，所以最小二乘是 MLE。X 满列秩时正规方程 XᵀXβ_hat=Xᵀy 给唯一系数；实际应使用合适求解或 QR/SVD，而非显式逆 Gram。模块 12 给投影解释。

条件于 X，误差零均值、Cov=σ²I 时，代入系数规则可推 E[β_hat|X]=β，Cov[β_hat|X]=σ²(XᵀX)⁻¹。n>d、满秩时，残差平方和/(n−d) 无偏：投影在期望中留下 n−d 个残差方向。高斯条件另支持精确 t 系数区间。误差 Cov 不同时，即使系数仍能算，方差分析也需调整。

截距列使精确无约束最小二乘残差和零，Xᵀr=0 可检查正规方程。拟合后残差一般相关，不能把原观测误差独立条件错误用于每个残差。查看残差与拟合值、特征、组、时间的关系，找均值、尺度或依赖错配。小训练残差不证明正确噪声、区间校准或未来准确。

固定新 x₀ 的拟合条件均值方差为 σ²x₀ᵀ(XᵀX)⁻¹x₀。独立新结果增加 σ²，预测方差为 σ²[1+x₀ᵀ(XᵀX)⁻¹x₀]。所以均值区间不同于结果预测区间。范围外外推可放大误差与模型错设；公式能延伸不证明真实关系线性。系数仅在另有因果条件时才能解释干预。

二值响应的 logistic 模型为 P(Y_i=1|x_i)=s(z_i)，z_i=x_iᵀβ，s=1/(1+exp(−z))。条件独立 Bernoulli 给负 logL=Σ[log(1+exp z_i)−y_i z_i]。梯度 Xᵀ(s−y)，Hessian XᵀWX、w_i=s_i(1−s_i)≥0，所以凸。完全分离可能没有有限无惩罚 MLE，凸不保证达到最小值。稳定 softplus/对数计算防止大 logit 溢出。

实验 2 在训练上拟合合成线性和 logistic 模型，报告保留数据平方误差、log 损失、Brier 与粗校准分箱。校准问总体频率是否匹配预测，不等于分类准确率高。小箱噪声大，有限表不能证明精确总体校准。线性同方差系数区间也不能未经推导搬到 logistic。

训练拟合参数，验证选模型、惩罚和停止，最终测试在所有选择冻结后用独立目标样本评估。用测试结果调参会把它变选择数据。改名文件或事后计算标准误不能撤销所选规则与结果的依赖。真正最终评估需要未使用样本或正确选择感知程序。

标准化、插补、特征选择及学习变换应在每个训练部分拟合，再用于验证/测试。用测试标签选特征直接泄漏；用其特征拟合预处理也改变所声称归纳协议。某些转导任务刻意用无标签目标输入，但应明确另一个任务。交叉验证需保留组/时间结构，并在每个训练折重新拟合预处理。

交叉验证平均保留折分数以支持模型选择，但训练集重叠，折分数一般依赖。其经验标准差除以 √折数不自动成为有效独立重复标准误。嵌套交叉验证将内调参与外评估分开，目标是完整选择流程；仍依赖总体、单位与拆分机制。更多折不消除分布漂移或主体泄漏。

m 个预指定有效 α/m 检验，通过并集界控制至少一个假拒绝概率≤α，称 Bonferroni。无需检验间独立，但每项有效且事件族受覆盖。m=20、α=.05 则阈值 .0025。未记录自适应搜索、切换结果和可选停止不自动属于固定族。错误发现率是另一目标，需另有程序。

类似地，每个 m 区间失败≤α/m，则全覆盖概率≥1−α。有界均值半径使用 log(2m/α)，由模块 24 推导。宽保守界是充分声明，不是精确实际样本量需求。固定模型保留准确率、选择模型验证结果和指定族同时保证是不同对象；模块 30 将从同逻辑推学习界。

实验报告列问题、目标量、独立单位、抽样分配、固定主要分析、排除、多重比较计划、效应、区间/检验解释与限制。数值复现不能暗示种子保证总体有效。观察关联、统计显著、预测校准与因果效应是不同结论，应准确写证据建立哪一项及其条件。

实验 1 · 自助法与精确置换

下载 lab1_bootstrap_and_permutation.py

"""IID two-group percentile bootstrap and an exact exchangeability test."""
import itertools
import math
import random


def mean(xs):
    return math.fsum(xs)/len(xs)


def quantile(xs, u):
    ordered = sorted(xs)
    index = (len(xs)-1)*u
    lo = int(index)
    hi = min(lo+1, len(xs)-1)
    return ordered[lo]+(index-lo)*(ordered[hi]-ordered[lo])


if __name__ == "__main__":
    a, b = [2., 3., 4., 5.], [0., 1., 2., 3.]
    observed = mean(a)-mean(b)
    rng, replicates = random.Random(26026), 5000
    boot = [mean([rng.choice(a) for _ in a])-mean([rng.choice(b) for _ in b]) for _ in range(replicates)]
    lo, hi = quantile(boot, .025), quantile(boot, .975)
    print(f"IID groups: A={a}, B={b}; difference={observed:.6f}")
    print(f"seed=26026; percentile bootstrap B={replicates}, approximate interval=[{lo:.6f}, {hi:.6f}]")
    combined = a+b
    stats = []
    for indices in itertools.combinations(range(8), 4):
        chosen = set(indices)
        stats.append(mean([combined[i] for i in chosen])-mean([combined[i] for i in range(8) if i not in chosen]))
    extreme = sum(abs(t) >= abs(observed)-1e-12 for t in stats)
    print(f"Exact label permutations={len(stats)}; inclusive two-sided extreme={extreme}; p={extreme/len(stats):.6f}")
    print("Null: the two groups share an exchangeable distribution, or random assignment with a sharp no-effect null.")
    print("Equal means alone with unequal distributions do not justify unrestricted label exchangeability.")
    print("A small percentile bootstrap sample has no automatic exact 95% coverage guarantee.")

IID groups: A=[2.0, 3.0, 4.0, 5.0], B=[0.0, 1.0, 2.0, 3.0]; difference=2.000000
seed=26026; percentile bootstrap B=5000, approximate interval=[0.500000, 3.500000]
Exact label permutations=70; inclusive two-sided extreme=10; p=0.142857
Null: the two groups share an exchangeable distribution, or random assignment with a sharp no-effect null.
Equal means alone with unequal distributions do not justify unrestricted label exchangeability.
A small percentile bootstrap sample has no automatic exact 95% coverage guarantee.

枚举 70 个分配，保留数值相同但索引不同的项。比较原假设分布与自助估计分布，说明各自条件。

实验 2 · 回归残差与概率校准

下载 lab2_regression_and_calibration.py

"""OLS and logistic fitting on specified independent synthetic data; NumPy 1.26.4."""
import numpy as np


def sigmoid(z):
    z = np.asarray(z)
    out = np.empty_like(z, dtype=float)
    positive = z >= 0
    out[positive] = 1/(1+np.exp(-z[positive]))
    exp_z = np.exp(z[~positive])
    out[~positive] = exp_z/(1+exp_z)
    return out


if __name__ == "__main__":
    rng = np.random.default_rng(26027)
    x = rng.uniform(-2, 2, 240)
    design = np.column_stack([np.ones(len(x)), x])
    y = 1+2*x+rng.normal(0, .5, len(x))
    train, test = np.arange(160), np.arange(160, 240)
    fitted, _, rank, _ = np.linalg.lstsq(design[train], y[train], rcond=None)
    residual = y[train]-design[train]@fitted
    sigma2 = residual@residual/(len(train)-2)
    covariance = sigma2*np.linalg.solve(design[train].T@design[train], np.eye(2))
    print("NumPy baseline1.26.4; seed26027; fixed160 training /80 test observations")
    print(f"OLS rank={rank}; intercept={fitted[0]:.6f}; slope={fitted[1]:.6f}")
    print(f"estimated coefficient SEs={np.sqrt(np.diag(covariance)).round(6).tolist()}; residual variance={sigma2:.6f}")
    print(f"normal equation residual max={np.max(np.abs(design[train].T@residual)):.3e}; held-out MSE={np.mean((y[test]-design[test]@fitted)**2):.6f}")
    probabilities = sigmoid(-.3+1.2*x)
    labels = (rng.random(len(x)) < probabilities).astype(float)
    weights = np.zeros(2)
    for _ in range(2500):
        p = sigmoid(design[train]@weights)
        gradient = design[train].T@(p-labels[train])/len(train)
        weights -= .2*gradient
    logits = design[test]@weights
    predictions = sigmoid(logits)
    nll = np.mean(np.logaddexp(0, logits)-labels[test]*logits)
    brier = np.mean((predictions-labels[test])**2)
    print(f"Logistic intercept={weights[0]:.6f}; slope={weights[1]:.6f}; final train gradient norm={np.linalg.norm(gradient):.3e}")
    print(f"Held-out NLL={nll:.6f}; Brier score={brier:.6f}")
    for lo, hi in ((0., .25), (.25, .5), (.5, .75), (.75, 1.)):
        mask = (predictions >= lo) & (predictions < hi)
        if np.any(mask):
            print(f"Calibration bin[{lo:.2f},{hi:.2f}): n={mask.sum()}, mean prediction={predictions[mask].mean():.6f}, success fraction={labels[test][mask].mean():.6f}")
    print("Bins are finite diagnostics, not proof of population calibration; OLS SE formula assumes the stated independent homoscedastic model.")

NumPy baseline1.26.4; seed26027; fixed160 training /80 test observations
OLS rank=2; intercept=0.996048; slope=1.997925
estimated coefficient SEs=[0.040504, 0.032941]; residual variance=0.262094
normal equation residual max=4.863e-14; held-out MSE=0.217459
Logistic intercept=0.036903; slope=1.254244; final train gradient norm=5.212e-16
Held-out NLL=0.560591; Brier score=0.187537
Calibration bin[0.00,0.25): n=20, mean prediction=0.153584, success fraction=0.050000
Calibration bin[0.25,0.50): n=15, mean prediction=0.368301, success fraction=0.466667
Calibration bin[0.50,0.75): n=17, mean prediction=0.624004, success fraction=0.705882
Calibration bin[0.75,1.00): n=28, mean prediction=0.862986, success fraction=0.714286
Bins are finite diagnostics, not proof of population calibration; OLS SE formula assumes the stated independent homoscedastic model.

先推目标和梯度，检查正规方程残差，列系数标准误条件。有限校准箱噪声不能证明精确校准。

实验 3 · 覆盖、选择与独立单位

下载 lab3_coverage_selection_and_units.py

"""Coverage enumeration and deliberately invalid selection/unit shortcuts."""
import math
import random


def coverage(n, p, method):
    total = 0.
    for k in range(n+1):
        estimate = k/n
        radius = math.sqrt(math.log(40)/(2*n)) if method == "Hoeffding" else 1.959963984540054*math.sqrt(estimate*(1-estimate)/n)
        if estimate-radius <= p <= estimate+radius:
            total += math.comb(n, k)*p**k*(1-p)**(n-k)
    return total


if __name__ == "__main__":
    for p in (.5, .01):
        for method in ("Wald plug-in", "Hoeffding"):
            print(f"Exact model coverage n=20, p={p:.2f}, {method} nominal95%: {coverage(20,p,method):.6f}")
    rng = random.Random(26028)
    candidates, n_val, n_test = 50, 100, 2000
    val = [sum(rng.random() < .5 for _ in range(n_val))/n_val for _ in range(candidates)]
    winner = max(range(candidates), key=val.__getitem__)
    test = sum(rng.random() < .5 for _ in range(n_test))/n_test
    print(f"50 identical population-accuracy .5 candidates: selected validation={val[winner]:.6f}; independent untouched test={test:.6f}")
    print("Calling selected validation a final test is invalid. Repeatedly selecting on the final test creates the same problem.")
    clusters, copies, p = 100, 10, .3
    print(f"Benchmark with {clusters} independent Bernoulli subjects, {copies} copied rows each:")
    print(f"wrong row SE={math.sqrt(p*(1-p)/(clusters*copies)):.6f}; correct subject SE={math.sqrt(p*(1-p)/clusters):.6f}")
    print("Repair: keep subject groups intact in splitting/resampling, use the subject as the independent unit, and predefine the comparison.")
    observed_difference, n, sigma = .01, 1000000, 1.
    z = observed_difference*math.sqrt(n)/sigma
    print(f"Known-sigma Gaussian mean test: effect={observed_difference:.3f}, n={n}, z={z:.1f}, two-sided p={math.erfc(abs(z)/math.sqrt(2)):.3e}")
    print("A small p-value does not make this .01-unit effect practically large or establish causality.")

Exact model coverage n=20, p=0.50, Wald plug-in nominal95%: 0.958611
Exact model coverage n=20, p=0.50, Hoeffding nominal95%: 0.997423
Exact model coverage n=20, p=0.01, Wald plug-in nominal95%: 0.182050
Exact model coverage n=20, p=0.01, Hoeffding nominal95%: 1.000000
50 identical population-accuracy .5 candidates: selected validation=0.630000; independent untouched test=0.493500
Calling selected validation a final test is invalid. Repeatedly selecting on the final test creates the same problem.
Benchmark with 100 independent Bernoulli subjects, 10 copied rows each:
wrong row SE=0.014491; correct subject SE=0.045826
Repair: keep subject groups intact in splitting/resampling, use the subject as the independent unit, and predefine the comparison.
Known-sigma Gaussian mean test: effect=0.010, n=1000000, z=10.0, two-sided p=1.524e-23
A small p-value does not make this .01-unit effect practically large or establish causality.

解释罕见率覆盖崩塌、选择乐观性和副本标准误。舍入到 1.000000 的覆盖只是显示，不证明永远不漏。

练习 1–12 必做，选学 13–14 在十二小时以外加 35 分钟。

读 MIT 经典推断和作者网站 Introduction to Statistical Learning 的回归、重抽样及多重检验章。本课推导初等覆盖、配对方差和 logistic 导数。

反演事件成区间，定义有效尾，推配对方差和回归目标，写不泄漏的比较协议，解释重抽条件与不能建立的因果结论。

结业任务：修复罕见率 Wald 确定性、副本基准与测试调停止。解释显著性前先报效应、带条件程序与独立数。

进入下一课的标准：能为统计评估程序与限制辩护。下一课连接概率、编码与对数损失。返回总览。