3.4 模型评估与交叉验证

3.4 模型评估与交叉验证

本节摘要：模型不是拟合出来的，是验证出来的。RMSEC 训练误差、RMSECV 交叉验证误差、RMSEP 独立测试误差三本账各记各的，其中 RMSECV 最诚实；杠杆值抓"位置极端"的样本，残差抓"成分异常"的样本，两把探针互补。本节还有一个真实的翻车记录：忘掉浓度中心化，预测整体偏移一个均值。

3.3 节结尾那组漂亮的 RMSEP 有个前提：测试集是独立另配的。真实项目里标样金贵，独立测试集常常舍不得留，日常全靠交叉验证把关——而交叉验证用错了姿势，比不验证更危险（它会给你一种"验过了"的错觉）。本节把误差的会计学补全，它上承 3.3 的 PLS 实现，直接通往第 6 章的参数调优方法论。

一、三本误差账

RMSEC（校正均方根误差）：模型在自己吃过的数据上的误差。它天生乐观——参数就是照这批数据调的，过拟合的模型 RMSEC 照样漂亮。RMSECV（交叉验证误差）：把校正集切成几折，轮流留一折当"临时测试集"，误差在没参与训练的样本上累计。它反映的是"模型换个样本还行不行"，选成分数、选预处理全看它。RMSEP（预测误差）：真正独立的测试集，只在最终验收时动用一次。三者的大致关系健康时是 RMSEC 小于等于 RMSECV 小于等于 RMSEP；RMSECV 与 RMSEC 差距拉大就是过拟合的警报。

用 3.3 的同一套数据与 PLS 实现跑这三本账：

# 三本误差账：RMSEC 与五折 RMSECV
import numpy as np
rng = np.random.default_rng(42)
p = 200
x = np.linspace(1100, 1250, p)
s_ana = 0.8*np.exp(-0.5*((x-1180)/12)**2)
s_int = 0.5*np.exp(-0.5*((x-1210)/20)**2)
n = 40
c = rng.uniform(0, 10, n)
ci = rng.uniform(2, 6, n)
X = np.array([c[i]*s_ana + ci[i]*s_int + rng.normal(0, 0.02, p) for i in range(n)])

def pls_fit(X, y, ncomp):
    mu, my = X.mean(axis=0), y.mean()
    Xc, yc = X - mu, y - my
    Wl, Pl, Ql = [], [], []
    for _ in range(ncomp):
        w = Xc.T @ yc; w /= np.linalg.norm(w)
        t = Xc @ w
        pl = Xc.T @ t / (t @ t); q = yc @ t / (t @ t)
        Xc -= np.outer(t, pl); yc -= q*t
        Wl.append(w); Pl.append(pl); Ql.append(q)
    return mu, my, np.array(Wl).T, np.array(Pl).T, np.array(Ql)
def pls_pred(m, X):
    mu, my, W, P, Q = m
    return my + (X - mu) @ W @ np.linalg.inv(P.T @ W) @ Q
def rmse(a, b): return float(np.sqrt(np.mean((a-b)**2)))

def kfold_rmsecv(X, y, k=5, ncomp=2):
    idx = np.arange(len(y)); rng.shuffle(idx)
    folds = np.array_split(idx, k)
    errs = []
    for f in folds:
        tr = np.setdiff1d(idx, f)                  # 其余四折当训练集
        m = pls_fit(X[tr], y[tr], ncomp)
        errs.append((pls_pred(m, X[f]) - y[f])**2)
    return float(np.sqrt(np.concatenate(errs).mean()))

print(f"RMSEC(k=2) = {rmse(pls_pred(pls_fit(X, c, 2), X), c):.3f}")
for kk in [1, 2, 3, 4]:
    print(f"k={kk}: RMSECV={kfold_rmsecv(X, c, 5, kk):.4f}")
​

输出：

RMSEC(k=2) = 0.005
k=1: RMSECV=0.3643
k=2: RMSECV=0.0056
k=3: RMSECV=0.0060
k=4: RMSECV=0.0060
​

选成分数的标准答案就藏在这五行里：k=1 误差 0.364（少一个成分，干扰没讲完），k=2 落到 0.0056 的噪声水平，k=3、k=4 不再改善（0.0060 附近）——**RMSECV 曲线的拐点就是成分数的答案**，拐点之后再无改善的成分全是噪声。注意 RMSEC 0.005 与 RMSECV 0.0056 几乎贴着——这批模拟数据信噪比极高；真实数据两者差个两三倍是常态，差十倍就该怀疑过拟合了。

交叉验证的姿势也有讲究。五折够用、十折更稳，留一法（每折一个样本）在样本少时诱人，但它的高方差会让 RMSECV 抖动。最要命的是分层陷阱：如果四十个样本里前二十个是低浓度、后二十个是高浓度，顺序切折会让某些折全是单一浓度——交叉验证前先打乱或按浓度分层，否则验证结果是假的。

二、翻车记录：忘掉浓度中心化

下面这个 bug 值得每位第一次写多元回归的人亲眼看一次。PCR 实现里光谱矩阵中心化了，浓度向量忘了：

# 一个真实的翻车：X 中心化了，y 没有
def pcr_fit_bug(X, y, ncomp):
    mu = X.mean(axis=0)
    U, S, Vt = np.linalg.svd(X - mu, full_matrices=False)   # X 中心化了
    T = U[:, :ncomp] * S[:ncomp]
    bt = np.linalg.lstsq(T, y, rcond=None)[0]      # y 忘了减均值！
    return mu, Vt[:ncomp], bt

mu, V, bt = pcr_fit_bug(X, c, 2)
pred_bug = (X - mu) @ V.T @ bt
print(f"buggy PCR RMSEC={rmse(pred_bug, c):.3f}")
print(f"pred mean={pred_bug.mean():.3f}  truth mean={c.mean():.3f}")
​

输出：

buggy PCR RMSEC=5.336
pred mean=0.000  truth mean=5.336
​

预测均值 0.000，真值均值 5.336——模型整体偏移了一个浓度均值，RMSEC 直接飙到 5.3。机理一句话：中心化后的得分 T 均值为零，回归里又没有截距项，输出被迫以零为中心。这个 bug 的阴险之处在于它有"半个体面"：预测与真值的相关系数可能依然很高，散点图上是一条漂亮的斜线，只是整体平移——看相关系数报告的人根本发现不了。防线只有一条：任何回归模型上线前，先对账预测均值与真值均值，两行代码的事。

三、两把探针：杠杆值与残差

模型验收还有最后一个问题：哪些样本在拖后腿。两把探针分工明确——杠杆值量样本在得分空间里的位置有多极端（是否外推），残差量模型解释不了的部分有多大（是否含异物）：

# 杠杆值与残差：两类离群样本的鉴别
s_odd = np.exp(-0.5*((x-1140)/8)**2)              # 校正集里没有的异物峰
X[7] = X[7] + 2.0*s_odd                           # 7 号样本混入异物
m3 = pls_fit(X, c, 2)
T = (X - m3[0]) @ m3[2]                           # 得分矩阵
H = T @ np.linalg.inv(T.T @ T) @ T.T              # 帽子矩阵
lev = np.diag(H)                                  # 每样本杠杆值
resid = c - pls_pred(m3, X)                       # 每样本残差

i_max = int(np.argmax(lev))
print(f"max leverage sample #{i_max}, lev={lev[i_max]:.3f}, resid={resid[i_max]:+.3f}")
print(f"sample #7: lev={lev[7]:.3f}, resid={resid[7]:+.3f}")
print(f"median |resid|={np.median(np.abs(resid)):.4f}  mean lev={lev.mean():.3f}")
​

输出：

max leverage sample #11, lev=0.136, resid=+0.013
sample #7: lev=0.041, resid=-0.105
median |resid|=0.0059  mean lev=0.050
​

读这组数的姿势：杠杆值均值 0.050（等于成分数除以样本数，2 除以 40，理论值），警戒线通常画三倍均值即 0.150——样本 11 的 0.136 已在警戒线边缘，它浓度或干扰恰好在分布的边上，属于位置极端：模型在这种样本上属于外推，预测要打折扣。样本 7 则是另一番景象：杠杆值平平（0.041），残差却冲到负 0.105——中位残差才 0.006，它放大了近二十倍——它的谱里混进了校正集没见过的异物峰，模型空间里根本没这个方向，这是成分异常。两把探针画出诊断图的四个象限：高杠杆低残差（外推但模型跟得上）、高杠杆高残差（既极端又异常，优先剔除复核）、低杠杆高残差（异物或浓度标错）、双低（好公民）。离群样本不是一律删除——先查它的来历，删除一个"真离群"等于销毁一条证据。

本节要点回顾

• RMSEC 乐观、RMSECV 诚实、RMSEP 一锤定音，三者分摊拉开即过拟合警报；
• 成分数看 RMSECV 拐点，实验里 k=2 之后曲线躺平，多出的成分全是噪声；
• 交叉验证先打乱或分层，顺序切折的验证结果是假的；
• 忘掉浓度中心化的 bug：预测整体偏移一个均值，相关系数却依旧漂亮——上线前先对账均值；
• 杠杆值抓位置极端，残差抓成分异常，四个象限分类处置，删除样本前先查来历。

鉴定台的四步到此走完。下一章进入高阶技术：重叠峰的解卷积与机器学习鉴定——但请带着本节的纪律进去，高阶方法是放大器，验证不严的模型会被它放大得更离谱。