简介：这份数理统计（孙海燕）大作业文档面向学习数理统计、模式识别与机器学习基础的高校学生及课程教师，围绕经典鸢尾花数据集展开完整分析，可作为课程作业撰写、统计方法复习与实验报告参考。文档为1个docx文件，压缩包约511KB，体量轻便，内容涵盖摘要、引言、数据及符号说明、算法原理分析、数据处理分析与结果讨论等部分，并串联马氏距离、高斯混合模型、主成分分析、线性判别分析和刀切法等知识点。读者可据此了解如何对鸢尾花四维属性做预处理、降维去噪、聚类分类与模型评估，特别是主成分分析与线性判别分析的差异、混合高斯模型在聚类中的优势，以及用刀切法比较降维程度和降维方法。目前已有1429人学习，适合需要一份结构完整、方法链条清晰的大作业范例来查漏补缺的读者。
1. 降维的收益不在维度，而在判错率
四个特征压到两个，信息必然有损，判错率却可能下降——这不是玄学。降维扔掉的往往是三类东西：高度共线的冗余方向、方差小到接近噪声的方向、以及让协方差矩阵求逆变得病态的方向。鸢尾花数据集是验证这件事的最小样本：150 条记录、3 个类别、4 个数值特征，花瓣长度与花萼长度的相关系数高到 0.87 附近，最后一个主成分的方差小到几乎可以忽略。一份数理统计课程的大作业把这条链路完整走了一遍：用主成分分析和线性判别分析降维，用混合高斯模型在极大似然意义下建模，用马氏距离做判别，最后用刀切法比较不同降维方法与不同降维程度的判错率差异。下面拆的是每一步为什么这么设计、以及手写实现时那些藏在协方差矩阵里的细节。
2. 马氏距离：为什么直接量欧氏距离会在鸢尾花上骗你
2.1 欧氏距离默认了三件不成立的事
欧氏距离把每个特征当成一根等长、互相垂直、互不影响的坐标轴。这三条默认在鸢尾花上都不成立。四个特征量纲一致（都是厘米），尺度这一条勉强过关；但方差差得远——花瓣长度的样本标准差接近 1.8 厘米，花萼宽度只有 0.4 厘米左右，同样的偏移量在不同轴上代表的「异常程度」完全不同。更麻烦的是相关性：花萼长度和花瓣长度的相关系数在 0.87 附近，花瓣越长的样本，花萼往往也越长，这两根坐标轴记录的很可能是同一件事。
后果很具体。拿两个花瓣长度差 1 厘米的样本算距离，欧氏距离会把这份差异在两根高度相关的轴上各计一次，于是这个方向的差异被系统性放大；反过来，两个花瓣长度接近、只在花萼宽度上差一点点的样本，距离显得很小。分类器看到的「远近」和真实的「像不像」脱了钩。这类问题的表现不是模型不收敛，而是收敛到一个被相关性扭曲的度量上，判错率卡在一个不上不下的水平，怎么调参数都下不去。
只做标准化解决不了这件事。标准化是把每一维除以自己的标准差，让方差都变成 1，它处理的是尺度，不处理方向之间的耦合。要把「花瓣长度和花萼长度说的是一回事」这个事实编码进距离里，必须动协方差矩阵。
2.2 拆开 Σ⁻¹：马氏距离其实是「先转轴、再缩放」
对均值 μ、协方差 Σ 的多变量向量，马氏距离的平方写成 d²(x, μ) = (x − μ)ᵀ Σ⁻¹ (x − μ)。把这个式子拆开看，它的几何含义就藏不住了。协方差矩阵是对称正定的，可以做特征分解 Σ = UΛUᵀ，于是 Σ^(−1/2) = UΛ^(−1/2)Uᵀ，代入得：
d²(x, μ) = ‖Λ^(−1/2) Uᵀ (x − μ)‖²
读法很直接：Uᵀ 先把样本旋转到协方差矩阵的主轴上，Λ^(−1/2) 再把每个主轴按各自的标准差归一，最后在归一后的空间里量一把普通欧氏距离。旋转那一步在消除相关性，缩放那一步在消除方差差异。所以马氏距离不是什么新奇的度量，它就是把数据中心化、摆正、按各自尺度拉伸之后，再量欧氏距离。
三条常见性质的来源也就清楚了。当 Σ = I 时 U 和 Λ 都是单位阵，退化为欧氏距离；尺度无关，因为量纲被 1/σ 除掉了；考虑特征间联系，因为旋转把相关方向对齐了。
这里有个容易忽略的连接：白化变换 UΛ^(−1/2)Uᵀ 和主成分分析用的是同一个分解。PCA 取 Λ 最大的那几个方向并保留原尺度，马氏距离则用上全部方向并按 √λ 缩放。两者共享一套特征结构，区别只是「保留哪些方向」和「是否缩放」。
再往前走一步是判别。若各类共享同一个协方差 Σ、先验相同，判别规则 argmin_c (x − μ_c)ᵀ Σ⁻¹ (x − μ_c) 展开后是 xᵀΣ⁻¹x − 2μ_cᵀΣ⁻¹x + μ_cᵀΣ⁻¹μ_c，第一项与类别无关，可整体约掉，剩下的判别函数对 x 是线性的。这和线性判别分析在协方差相同假设下得到的形式是同一件事。若每类用自己的 Σ_c，二次项消不掉，判别边界变成二次曲面，那就是二次判别分析。马氏距离、LDA、QDA 三者其实是同一条推导链上的三个位置。
2.3 手写实现：为什么不该用 inv 求逆
理解公式和让代码不出现 NaN，中间隔着一层数值稳定性。下面这段是手写版本里最核心的一小块：
import numpy as np def mahalanobis_sq(X, mu, Sigma, ridge=1e-6): """X: (n, p) 样本矩阵；mu: (p,) 类均值；Sigma: (p, p) 协方差 返回 (n,) 每个样本到 mu 的平方马氏距离""" p = X.shape[1] # 对角扰动：协方差接近奇异时，直接求逆会被最小特征值方向主导 S = Sigma + ridge * np.eye(p) d = X - mu # 用 solve 解 S z = d^T，而不是 inv(S) @ d： # 少一次矩阵乘法，且不会显式构造条件数极差的逆矩阵 z = np.linalg.solve(S, d.T) # (p, n) return np.einsum('ij,ji->i', d, z) # 逐样本 (x-mu)^T S^{-1} (x-mu)ridge是一个很小的对角扰动，量级取相对于协方差特征值足够小的值，标准化后的数据上 1e-3 也是安全的；它的作用不是提高精度，而是防止 Σ 奇异或近奇异时solve返回一堆毫无意义的大数。np.linalg.solve与np.linalg.inv的复杂度同阶，但前者少一次完整的矩阵乘法，且在条件数偏大时误差更小。最后的einsum做的是逐行内积，等价于 (d @ S⁻¹ * d).sum(axis=1)，省掉构造完整矩阵的中间步骤。
还有一层比数值稳定性更根本的隐患：协方差矩阵本身是估出来的。4 个特征、每类 50 个样本，Σ 有 10 个自由参数，估计还算稳。维度一旦升到几十，参数数量是 p(p+1)/2 的平方级增长，样本根本喂不饱，Σ⁻¹ 会被最小特征值方向主导，距离被噪声放大。这解释了为什么降维和马氏距离总是绑在一起出现——降维不只是为了算得快，是为了让 Σ⁻¹ 可估。后面用主成分分析和线性判别分析先降到一维或二维，再求马氏距离，本质上是在保证这个逆矩阵有意义。
提示：如果想确认实现对不对，拿 Σ = I 跑一遍，结果应当与欧氏距离逐元素相等；再拿 Σ 为对角阵跑一遍，结果应当等于各维差值除以标准差后的平方和。这两个退化情形是最省事的自检手段。
3. PCA 与 LDA：有没有标签，决定了你在找方差还是找可分性
3.1 主成分分析：把「信息量大」翻译成「方差大」
主成分分析要解决的是一个纯优化问题：找一个单位方向 u，让数据投影上去之后的方差最大。投影后的方差是 uᵀΣu，约束是 ‖u‖ = 1，写出拉格朗日函数求导，得到 Σu = λu。也就是说，最优方向就是协方差矩阵最大特征值对应的特征向量，最大方差就是那个特征值。要找第二个与第一方向正交且方差最大的方向，答案就是第二大特征值对应的特征向量，依此类推。
这条推导解释了 PCA 的全部行为：它只看协方差矩阵，只关心数据在哪个方向上散得开，完全不知道标签的存在。方差大等于信息量大这个假设，在多数场景下成立，但它并非普适——如果某个方向方差极大却和类别无关（比如测量仪器的系统性漂移），PCA 会优先保留这个方向，等于主动放大噪声。
鸢尾花上还有一个预处理细节值得强调：原作业明确用相关矩阵代替协方差矩阵。四个特征虽然都是厘米，但方差差一个量级，不标准化的话第一主成分几乎就是花瓣长度单变量。用相关矩阵 R（即标准化后再求协方差）等价于先做 z-score 再分解，这样每一维在起跑线上权重相同。代价是把「本来方差大」这一信息也抹掉了，如果特征量纲一致且希望保留绝对尺度，就不该这么做。
3.2 线性判别分析：一个比值同时管两件事
线性判别分析的思路朴素到可以一句话说完：把样本投影到一条直线上，让同类投影点尽量挤在一起，异类投影点尽量分开。为了实现这句话，定义两个散度矩阵——类内散度 S_w = Σ_c Σ_{x∈c} (x − μ_c)(x − μ_c)ᵀ，类间散度 S_b = Σ_c n_c (μ_c − μ)(μ_c − μ)ᵀ，其中 n_c 是第 c 类样本数、μ 是全局均值。目标是最大化
J(u) = uᵀ S_b u / uᵀ S_w u
分子是投影后的类间方差，分母是投影后的类内方差，比值越大，投影后越可分。这个比值对 u 的缩放不变，所以约束 ‖u‖ = 1 或者 uᵀS_w u = 1 都行。对后者用拉格朗日乘子法，求导得到 S_b u = λ S_w u，即广义特征值问题；两边左乘 S_w⁻¹ 得到 S_w⁻¹S_b u = λu，取最大特征值对应的特征向量即可。
二分类时有个更直观的闭式解：u ∝ S_w⁻¹(μ₁ − μ₂)。注意这个形式——它在「两类均值差」外面套了一个 S_w⁻¹。这个 S_w⁻¹ 正是 2.2 节里出现过的那个东西，所以 LDA 的方向可以读成「在马氏度量下，把两类均值拉开的方向」。类内散度大的方向会被压低权重，类内散度小的方向会被放大权重，这正是「兼顾放大类间差异、缩小类内差异」在数学上的落地方式。
和 PCA 的根本差别也在这里：PCA 的目标函数里没有 y，LDA 的 S_b 里全是类均值。有标签的时候，LDA 用了更多信息，理论上限更高；没标签的时候 LDA 根本无从下手。这不是孰优孰劣，是信息量的问题。
3.3 三分类最多降到二维，这是秩决定的硬约束
一个经常踩的坑：对三类鸢尾花调用 LDA 并设置降维到 3 维，程序报错，或者返回一个全零方向。原因不在实现，在秩。
S_b 是 C 个秩 1 矩阵 n_c (μ_c − μ)(μ_c − μ)ᵀ 之和，而这 C 个向量 (μ_c − μ) 带权重之后线性相关（它们的加权和为零），所以它们张成的子空间维数最多是 C − 1。于是 rank(S_b) ≤ C − 1，S_w⁻¹S_b 的非零广义特征值最多只有 C − 1 个。三分类的鸢尾花，判别方向最多两个，第三维投影过去是恒等于零的退化方向。
这条约束对所有判别式降维都成立：想降到 k 维，类别数必须满足 C ≥ k + 1。做降维维度对比实验时，这一条决定了实验表里 LDA 那一行最多只有 1 维和 2 维两列，再往上没有意义。反过来说，正因为 LDA 的上限是 2，用 4 维原始特征做留一法对比时，LDA 路线天然就是「降到 2 维」，不存在「LDA 降到 3 维」这种配置。
3.4 同一份数据上的两条投影路径
把两种方法都降到二维，投影出来的形状不一样。PCA 的第一主成分几乎沿着「整体尺寸」方向走，因为花瓣长度和花萼长度同向变化，这个方向的方差最大；第二主成分是剩下方差最大的正交方向，通常捕捉花瓣宽度与花萼宽度之间的对比。LDA 的第一判别方向则明显偏向于分开后两类，setosa 因为离得太远，几乎在任何方向上都能被一条直线切走。
这个差异在留一法结果里会变成一个很好看的现象：几乎所有判错都集中在第二类和第三类之间，第一类的判错率为零。原因不是算法偏心，而是数据结构如此——setosa 在四个特征上和另外两类完全没有重叠，剩下的两类在花瓣宽度、花瓣长度上有一段实打实的重叠区。理解了这一点，看到判错率 3% 左右就不会觉得是实现有问题，也不会盲目去调参。
下面是两个核心类的骨架，实现了「用相关矩阵分解」和「用 Cholesky 白化后分解」两条路径：
import numpy as np class PCA: def __init__(self, n_components=2): self.n_components = n_components def fit(self, X): Xc = X - X.mean(axis=0) Xs = Xc / Xc.std(axis=0, ddof=1) # 标准化，等价于改用相关矩阵 R = np.corrcoef(Xs, rowvar=False) w, V = np.linalg.eigh(R) # 对称矩阵专用，比 eig 快且稳 idx = np.argsort(w)[::-1] self.explained_ratio_ = w[idx] / w.sum() # 每个主成分解释的方差比例 self.components_ = V[:, idx[:self.n_components]].T return self def transform(self, X): Xc = X - X.mean(axis=0) Xs = Xc / Xc.std(axis=0, ddof=1) return Xs @ self.components_.T def fit_lda(X, y, n_components=2, eps=1e-6): n, p = X.shape classes = np.unique(y) mu_all = X.mean(axis=0) S_w = np.zeros((p, p)) S_b = np.zeros((p, p)) for c in classes: Xc = X[y == c] mu_c = Xc.mean(axis=0) S_w += (Xc - mu_c).T @ (Xc - mu_c) # 类内散度累加 diff = (mu_c - mu_all).reshape(-1, 1) S_b += len(Xc) * (diff @ diff.T) # 类间散度按类样本数加权 # 直接对 S_w^{-1} S_b 求特征值会得到非对称矩阵，可能出复数 # 先做 Cholesky 白化，把它转成对称矩阵再分解 L = np.linalg.cholesky(S_w + eps * np.eye(p)) Linv = np.linalg.inv(L) M = Linv @ S_b @ Linv.T # 对称，且特征值就是广义特征值 w, V = np.linalg.eigh(M) idx = np.argsort(w)[::-1][:n_components] return (Linv.T @ V[:, idx]).T # 变回原特征空间的方向PCA.fit里用np.linalg.eigh而不是eig，因为相关矩阵对称，eigh只读下三角、保证返回实数特征值并按升序排列，最后用[::-1]倒序取最大。explained_ratio_保留全部四个主成分的比例，方便判断降维到几维够用——鸢尾花上通常前两个加起来就超过九成。
fit_lda里最关键的两行是 Cholesky 白化。S_w正定时可以分解成 LLᵀ，令 M = L⁻¹ S_b L⁻ᵀ，则 M 对称，其特征值与 S_w⁻¹S_b 完全相同，特征向量经 L⁻ᵀ 变换回去就是原空间的判别方向。这样做避免了非对称矩阵特征分解可能出现复数的问题，也让eigh能正常使用。eps的作用是兜底：当某一维在类内几乎没有变化、S_w 接近奇异时，不加扰动 Cholesky 会直接失败。
注意：
eps不是越大越安全。取值超过 S_w 最小特征值的量级，判别方向会被扰动带偏，表现为降维后的类间距离莫名缩小。判断标准是看 S_w 的条件数，超过 1e6 就该怀疑数据本身是否适合做判别分析。
4. 高斯混合模型：极大似然为什么解不出闭式，EM 在补什么
4.1 对数似然里那个 log-sum，让闭式解消失
单个高斯成分的极大似然估计有漂亮的闭式解：均值就是样本均值，协方差就是样本协方差，把似然函数取对数、对参数求导令其为零，两步就能推出来。这也是为什么原作业里能直接写出 μ 和 Σ 的估计式。
混合模型就没这个运气了。K 个成分的密度是 p(x) = Σ_k π_k N(x | μ_k, Σ_k)，整个数据集的似然取对数之后是
log p(X) = Σ_i log ( Σ_k π_k N(x_i | μ_k, Σ_k) )
问题出在最外层那个 log 套在求和号外面。对 μ_k 求偏导时，求导算子要穿过 log，再穿过对 k 的求和，得到的形式里每一项都同时含有所有成分的参数。三个参数组（π、μ、Σ）互相耦合在同一个方程里，没有任何一组能单独解出来。这不是推导技巧不够，是结构上就没有闭式解。
顺便澄清一个常见混淆：单个成分的 MLE 形式（样本均值、样本协方差）在混合模型里依然会出现，但它只在「已知每个样本属于哪个成分」的前提下成立。而归属恰恰是未知的，这正是死结所在。
4.2 E 步与 M 步：责任度就是后验概率
EM 的思路是给这个死结补一个中间变量。既然不知道样本属于哪个成分，那就先猜一个概率分布，用猜的分布去算参数的期望，再拿新参数去更新猜的分布，来回迭代。
E 步算的是每个样本对每个成分的「责任度」，形式上是后验概率：
γ_ik = π_k N(x_i | μ_k, Σ_k) / Σ_j π_j N(x_i | μ_j, Σ_j)
分母是样本 i 在整个混合分布下的密度，分子是它单独落在第 k 个成分下的加权密度。γ_ik 越大，说明这个样本越像是由第 k 个成分生成的。所有 γ_ik 对 k 求和为 1，所以它是一组软分配权重，和 k-means 里非 0 即 1 的硬分配形成对照。
M 步拿这组权重去更新参数，形式与单高斯的 MLE 几乎一样，只是每个样本都带了权重：
π_k = N_k / n，其中 N_k = Σ_i γ_ik μ_k = (1/N_k) Σ_i γ_ik x_i Σ_k = (1/N_k) Σ_i γ_ik (x_i − μ_k)(x_i − μ_k)ᵀ
N_k 可以理解成「软计数」——不是数有多少个样本属于第 k 类，而是把所有样本的归属概率加起来。这三行更新保证了每轮迭代后对数似然单调不减，这是 EM 最重要的性质，也是判断实现是否正确的最直接依据：如果跑出来的对数似然出现下降，那一定是 E 步或 M 步写错了，不存在「正常波动」。
4.3 协方差形状：GMM 比 k-means 强在哪
k-means 的目标是最小化各点到簇心的欧氏距离平方和，这个目标隐含了三个假设：每个簇是球形的（各方向方差相同）、各簇大小相同、簇与簇的边界是垂直平分面。这三个假设在鸢尾花上都不成立。后两类的分布是斜着叠在一起的，花瓣长度和花瓣宽度的联合分布呈现明显的椭圆形状，长轴方向一致但中心错开。k-means 会沿着一条垂直于中心连线的直线切开，切掉的那部分里混着另一类的样本。
全协方差的 GMM 允许每个成分有自己的椭球形状——不同方向不同方差、椭球可以倾斜、可以大小不一。判别边界从超平面变成二次曲面，能够顺着椭球的边界弯过去。这就是「考虑了不同方向上方差不同，也考虑了不同维度之间相关性」这句话在几何上的具体含义。
有一个等价关系能把两者彻底串起来：k-means 是 GMM 在 Σ_k = σ²I（各向同性、各成分相同）、且 σ² → 0 时的极限，此时责任度 γ_ik 退化为 0/1 的硬分配，M 步的均值更新退化为簇内均值。所以 GMM 不是「另一种聚类算法」，它是 k-means 放松假设后的推广；反过来，k-means 是 GMM 的一个特例，速度快但形状假设强。
4.4 实现细节：初始化、正则与收敛判定
EM 只保证收敛到局部最优，不保证全局最优。初始化的好坏直接决定最后停在哪，所以实现里必须做两件事：固定随机种子保证可复现，以及多次重启取对数似然最大的那一组。下面这段是核心循环：
import numpy as np def gmm_em(X, K, n_iter=200, tol=1e-6, seed=0, reg=1e-6): rng = np.random.default_rng(seed) n, p = X.shape idx = rng.choice(n, K, replace=False) mu = X[idx].copy() # 随机挑 K 个点做初值 Sigma = np.array([np.cov(X, rowvar=False) + reg * np.eye(p)] * K) pi = np.full(K, 1.0 / K) loglik_old = -np.inf for _ in range(n_iter): # E 步：在对数域算加权似然，避免小概率下溢成 0 logp = np.zeros((n, K)) for k in range(K): d = X - mu[k] z = np.linalg.solve(Sigma[k], d.T) maha = np.einsum('ij,ji->i', d, z) logp[:, k] = (np.log(pi[k]) - 0.5 * maha - 0.5 * np.log(np.linalg.det(Sigma[k])) - 0.5 * p * np.log(2 * np.pi)) m = logp.max(axis=1, keepdims=True) # log-sum-exp 稳定化 gamma = np.exp(logp - m) gamma /= gamma.sum(axis=1, keepdims=True) loglik = float((m + np.log(np.exp(logp - m).sum(axis=1, keepdims=True))).sum()) # M 步：软计数加权更新 Nk = gamma.sum(axis=0) mu = (gamma.T @ X) / Nk[:, None] for k in range(K): d = X - mu[k] Sigma[k] = (gamma[:, k][:, None] * d).T @ d / Nk[k] + reg * np.eye(p) pi = Nk / n if abs(loglik - loglik_old) < tol: break loglik_old = loglik return pi, mu, Sigma, gammalogp在对数域累加，是因为密度值在高维下会小到 double 下溢成 0，一旦变成 0 再取对数就是负无穷。m = logp.max(...)是 log-sum-exp 技巧：先把每行的最大值提出来，剩下的部分取指数就不会溢出，最后再把最大值加回去。这个技巧在 GMM、隐马尔可夫、变分推断里反复出现，属于通用做法。
reg * np.eye(p)在协方差更新里也必须加。当某个成分只分到极少样本时（N_k 很小），样本协方差可能接近奇异，下一轮求逆就会炸，进而让这个成分吸走更多样本或者彻底塌缩成一个点。加对角扰动相当于给协方差一个先验，代价是略微高估方差，但在小样本下这个代价远小于崩溃。
tol的判断对象是平均对数似然的绝对变化量，不是相对变化量。如果数据量很大，tol=1e-6可能跑满 200 轮都不收敛，此时应该看的是参数变化或者直接改成相对判据。
5. 刀切法的排查清单：留一法为什么更贵却更可信
5.1 留一法在估计什么
刀切法的流程可以写成五步：从某一类的 n₁ 个训练样本里剔除一个，用剩下的样本建立判别规则；用这条规则去判别刚剔除的那个样本；重复直到这一类每个样本都被剔除并判别过一次，记下判错个数；对每个类都做一遍；最后把各类判错个数加起来除以总数，得到判错率。
严格来说，刀切法和留一交叉验证在术语上有细微差别：前者一般指每次删掉一个样本后重新计算统计量，用来估计这个统计量的偏差和方差；后者是用同一套循环来估计模型的泛化误差。在这份作业的场景里两者是同一套循环——每次留出一个样本、用其余样本重训、预测留出样本——所以可以互换着叫。真正需要区分的是「回报法」：把训练集本身拿去做预测再统计判错率，这样得到的数字一定偏乐观，因为它衡量的是模型在见过的数据上的表现，而模型参数正是被这些数据调出来的。
留一法几乎用上了全部 n−1 个样本训练，且每个样本都恰好被当作一次「没见过的数据」，所以它对泛化误差的估计近似无偏。代价是训练次数等于样本数，150 个样本就要重训 150 次。这正是「运算量大」的来源，也是它在小数据集上还能接受、在大数据集上必须换成 k 折的原因。
5.2 分层留一：三类各留一个更均衡
鸢尾花三类各 50 个样本，整体留一时每轮剔掉的样本类别分布基本均匀，看不出问题。但只要类别不均衡，整体留一就会出状况：某一轮剔除的刚好是稀有类里唯一的一个样本，训练集里这个类直接消失，判别规则根本没法学，这一轮的判错是必然的，但它反映的不是模型能力，而是划分方式的问题。
分层留一的做法是每一轮从每个类里各剔除一个样本，这样训练集的类别比例始终和原始数据一致。对三类鸢尾花而言，一轮同时留出三个样本、重训一次，训练次数从 150 次降到 50 次，计算量直接压到三分之一，估计的仍然是「模型在没见过的新样本上的表现」。这个改动在实现上只需把循环单位从「一个样本」换成「一组样本」。
下面这段是带分层和预处理内嵌的框架：
import numpy as np def stratified_loo_error(X, y, build_model): """build_model 接收训练集，返回一个带 fit/predict 的对象。 关键点：预处理必须放在 build_model 内部，用训练折 fit。""" classes = np.unique(y) folds = [np.where(y == c)[0] for c in classes] wrong = 0 total = 0 for i in range(max(len(f) for f in folds)): test_idx = [f[i] for f in folds if i < len(f)] mask = np.ones(len(y), dtype=bool) mask[test_idx] = False model = build_model() model.fit(X[mask], y[mask]) # 标准化、PCA/LDA 都在这里面完成 pred = model.predict(X[test_idx]) wrong += int((pred != y[test_idx]).sum()) total += len(test_idx) return wrong / total # 判错率框架本身只有十几行，所有风险都集中在build_model的实现里。如果build_model内部调用的标准化器是在传入的训练集上 fit 的，那就正确；如果是在调用之前对全量数据做了一次标准化再切片传入，那就泄漏了。这个区别在代码上只差一行，在结果上能把判错率压低好几个百分点。
5.3 五条踩坑记录
第一条：判错率异常低，低到接近零。现象是留一法跑出来判错率 0.7%，比预期好太多。原因几乎总是数据泄漏——标准化、主成分分析、判别分析的投影方向，如果是在全量数据上拟合的，留出样本的信息已经以均值、方差、投影方向的形式渗进了训练过程。解决方式是把整个预处理链路包进build_model，每一折用训练部分重新 fit，留出样本只做 transform。
第二条：每次跑的判错率都不一样，波动几个百分点。原因是 GMM 的随机初始化和 EM 的局部最优，换一个种子就落到另一个极值点。解决办法是固定种子并做多次重启，取对数似然最大的那组参数作为最终模型；实验报告里应该说明重启次数，否则数字不可复现。三分类、每类 50 个样本的规模下，5 次重启基本能稳定下来。
第三条：LDA 降到 3 维报错，或者投影结果有一列全零。原因是 3.3 节说的秩约束，rank(S_b) ≤ C − 1，三分类最多两个有效判别方向。解决办法是把n_components限制在min(C-1, p)以内，并在代码里加断言，别等到特征值出现数值噪声才反应过来。
第四条：马氏距离判别出现 NaN，或者判错率突然飙到 50% 以上。原因是某类的协方差矩阵接近奇异——留一法每一折都会少一个样本，如果恰好让某一维在类内完全不变（比如花瓣宽度只有 0.1 厘米的几档取值），方差估计变成零，求逆就崩了。解决方式是加对角扰动、用伪逆、或者干脆在这一折里降维后再算距离。判断信号是看 Σ 的条件数，超过 1e6 就要警惕。
第五条：留一法跑得极慢，一次实验十几分钟。原因是每一折都重新做
本文还有配套的精品资源，点击获取
