简介：一份围绕鸢尾花数据集展开的数理统计大作业完整文档，适合正在学习多元统计、模式识别或机器学习入门课程的高校学生参考。资源包含1个docx文档，压缩包约511KB，内容系统演示了马氏距离、混合高斯模型、主成分分析、线性判别分析和刀切法在花萼、花瓣共四个属性数据上的综合应用。文档从符号说明与样本数据选取整理入手，依次给出各类算法的原理推导与公式说明，并通过PCA与LDA降维、GMM聚类以及留一法模型评估，形成完整的数据处理分析和结果对比；其中还讨论了不同降维程度与降维方式的差异，能够帮助读者理解统计模型在实际数据上的建模、评估与选择思路。目前已有1426人学习下载，可作为撰写数理统计实验报告、复现经典分类聚类方法或准备课程答辩时的参考样板。
1. 五十年经典数据集遇上四把统计利器：这套组合拳到底在验证什么
一份数理统计大作业把鸢尾花数据集作为试验台，同时用上马氏距离、高斯混合模型、主成分分析、线性判别分析和刀切法，很容易被看作在罗列算法。但如果只盯着算法清单，就会错过真正值得搞懂的东西：这五样东西按「先降维、再建模、后判类、最后验证」的顺序排列，组成了一条完整的统计决策链。任何一个环节换掉，结果的含义都会变。
鸢尾花数据集很小，150 个样本、4 个数值属性、3 个类别，却是数理统计里最合适的最小实验场。它属性之间存在线性相关，量纲又不一致，于是有了降维的必要；属性基本符合高斯分布假设，于是高斯混合模型能派上用场；类别之间又不完全可分，所以需要一个不偏不倚的评价方式来比较不同降维方案。读懂这份作业，等于把「为什么数据决定方法，方法又如何在数据上兑现」这条主线完整看了一遍，比单独记住某个算法的公式更有复用价值。
2. 数据与马氏距离：为什么欧氏距离在四维特征上直接失灵
2.1 鸢尾花数据的真实结构：相关性与量纲如何影响距离
鸢尾花数据集包含 150 个观察值，分为 3 类，每类 50 个，每个观察值有花萼长度、花萼宽度、花瓣长度、花瓣宽度 4 个属性。数据从公开机器学习数据库中下载，列的顺序是花萼长度、花萼宽度、花瓣长度、花瓣宽度，最后一列是类别标签。符号上，称自变量为 x，它是 4 维输入向量；n 表示样本个数，取 150；μ 表示均值向量，Σ 表示协方差矩阵，R 表示相关矩阵。
这份数据的结构藏着一个陷阱：花萼长度分布在 4.3 到 7.9 厘米之间，花瓣宽度则只有 0.1 到 2.5 厘米的量级。如果用欧氏距离去度量两个样本的相似度，数值范围更大的变量会在平方和中占绝对主导，花萼长度和花瓣长度几乎决定距离，而花瓣宽度的差异被淹没。更隐蔽的是属性之间的相关性。花萼长度和花瓣长度、花瓣宽度之间存在明显的线性关联，这些变量反映的信息有重叠。欧氏距离假设各维度相互独立、等权重，这个假设在鸢尾花数据上从头到尾都不成立。
从类别结构看，Setosa 这一类和另外两类在花瓣属性上有非常清晰的边界，几乎可以线性切分；Versicolour 与 Virginica 则在花瓣长度和花瓣宽度上有较大重叠。这意味着后续方法要解决的核心问题不是「能不能分」，而是「在不完全可分的地方，如何度量相似度、如何建模密度、如何评价降维后信息保留了多少」。这三层问题环环相扣，第一环就落在距离度量上。
2.2 马氏距离的本质：旋转坐标系再按方差缩放
马氏距离的定义是一个均值为 μ、协方差矩阵为 Σ 的多变量向量 x 到该分布中心的距离，形式为 d(x, μ) = sqrt((x − μ)^T Σ^(−1) (x − μ))。表面看只是把欧氏距离里的单位矩阵换成了 Σ^(−1)，但这个替换改变了度量的几何含义：它先把坐标系按协方差矩阵旋转到各主轴上，再按每个方向上的方差做缩放，最后才计算欧氏距离。用一个类比来说，欧氏距离看的是「直线距离」，马氏距离看的是「在这个分布里有多异常」。
马氏距离有三个性质值得展开。第一，当 Σ = I 时，各项独立且方差相同，马氏距离退化为欧氏距离，所以欧氏距离只是马氏距离的特例。第二，马氏距离等价于对 x 和中心 μ 做标准化之后再求欧氏距离，标准化的尺度由协方差矩阵决定，因此它天然消除了量纲影响。第三，它对变量的线性变换具有不变性，也就是尺度无关，无论把厘米换成英寸还是把数据整体做线性旋转，距离排序不变。这三个性质合起来说明一件事：当类别分布不是正圆而是旋转椭球时，欧氏距离会同时忽略方向和展布，马氏距离则把椭球的形状直接编码进度量里。
在后续的判别流程里，高斯混合模型会对每个类别拟合出一个高斯成分，每个成分有自己的 μ 和 Σ。此时样本到底属于哪个类别，最合理的度量就是到该成分中心的马氏距离，因为 Σ 就是这个椭圆簇的形状。把马氏距离换回欧氏距离，等于强行把所有簇当成正圆，重叠区域里的样本会被系统性地误判。
2.3 落成代码：一个最小可用的马氏距离函数
下面是实现马氏距离的 Python 函数，输入一组样本、一个中心向量和一个协方差矩阵：
import numpy as np def mahalanobis_distance(x, mu, cov): # x: (n, d) 每个样本一行；mu: (d,) 中心；cov: (d, d) 协方差 # 用伪逆而不是逆，防止协方差矩阵奇异时直接报错 cov_inv = np.linalg.pinv(cov) diff = x - mu.reshape(1, -1) left = diff @ cov_inv return np.sqrt(np.einsum('ij,ij->i', left, diff))参数说明：x 的形状是 (n, d)，鸢尾花场景里通常是 150 行 4 列，但降维后可能变成 150 行 2 列；mu 是某个高斯成分的均值向量长度 d；cov 是该成分的协方差矩阵形状 (d, d)。einsum 这一步计算的是 left 和 diff 逐行求内积，相当于对角化地算出行向量 diff @ cov_inv @ diff.T 的结果，避免用矩阵乘法再取对角线的写法。返回的是每个样本到该中心的马氏距离数组。
提示：用伪逆 pinv 只是一种兜底手段。当协方差矩阵奇异时，它给出最小二乘意义下的距离，避免程序崩溃，但也会掩盖协方差估计本身退化的问题。真正要紧的是在建模阶段通过正则化防止协方差矩阵病态，而不是在距离计算阶段补救。这套思路在后面 GMM 的 EM 迭代里会再次出现。
3. 主成分分析与线性判别分析：降维的两条路线差在标签用没用
3.1 两种降维的分水岭：方差最大化和可分性最大化
主成分分析（PCA）和线性判别分析（LDA）都是线性降维方法，形式上都是找一个线性变换把 p 维数据压到低维空间，但目标函数完全不同。PCA 是无监督方法，不理会类别标签，只要求投影后的数据方差尽量大，因为方差被认为是信息量的代理。LDA 是有监督方法，由分类任务驱动，目标是把同类样本的投影点拉近、把异类样本的投影点推远，也就是同时放大类间离差、缩小类内离差。
这个差异直接解释了为什么在有标签的情况下 LDA 理论上比 PCA 效果好：PCA 找到的方差最大方向不一定能区分类别，可能两个类别在该方向上大致重合，投影后反而更难分；LDA 的优化目标天然就是「判别友好」的，它找到的方向直接服务于后续分类。在鸢尾花数据上，三个类别的形状类似，各类协方差结构差异不大，正好落在 LDA 的假设范围内，所以 LDA 的降维结果在分类任务上更有优势。
需要强调 LDA 并不是在所有场景下都优于 PCA。它的监督信息来自标签，如果标签噪声大、或类别数不准确，LDA 会把标签里的错误也学进去；PCA 则完全不依赖标签，在无监督探索、数据可视化、去噪场景下依然是更稳健的起点。两者的选择本质上是「要不要为分类目标牺牲无偏性」的取舍。
3.2 PCA 的源码实现：用相关矩阵等价于先标准化再分解
下面是 PCA 类的最小实现，接口设计是 fit 学习投影方向、transform 执行投影：
class PCA: def __init__(self, n_components=2): self.n_components = n_components def fit(self, X): # 先做标准化：(X - mean) / std self.mean_ = X.mean(axis=0) self.std_ = X.std(axis=0) X_std = (X - self.mean_) / self.std_ # 标准化后的协方差矩阵就是相关矩阵 cov = np.cov(X_std, rowvar=False) # 形状 (d, d) eigvals, eigvecs = np.linalg.eigh(cov) idx = np.argsort(eigvals)[::-1] # 特征值降序 self.components_ = eigvecs[:, idx[:self.n_components]] return self def transform(self, X): X_std = (X - self.mean_) / self.std_ return X_std @ self.components_代码里的关键逻辑有两点。第一，fit 里先标准化再求协方差矩阵，得到的其实就是相关矩阵，这避免花萼宽度（0.1 到 2.5 厘米）和花瓣长度（1.0 到 6.9 厘米）之间的量纲差异主导主成分方向。第二，np.linalg.eigh 返回升序排列的特征值和特征向量，所以要用 argsort 反转索引，再取前 n_components 列构成投影矩阵。投影后的方差正好等于对应特征值，因此特征值大小直接反映该主成分保留了多少信息量。
提示：PCA 的标准化步骤不能省。如果直接对原始协方差矩阵做特征分解，数值范围大的变量会主导第一个主成分方向，低方差但有用的变量可能被完全忽略。这不是实现细节问题，而是 PCA 定义的一部分——用相关矩阵就是默认「每个原始变量先被压缩到单位方差再竞争信息量」。这份大作业里专门强调用相关矩阵代替协方差矩阵来消除量纲影响，原因就在这。
3.3 LDA 的源码实现：广义特征值问题里藏着类别数上限
LDA 的求解逻辑比 PCA 多了一步：构造类内散度矩阵 Sw 和类间散度矩阵 Sb，然后解广义特征值问题。Sw 衡量每个类别内部样本的离散程度，Sb 衡量各类中心相对全局中心的离散程度。目标是让投影后的类间散度尽量大、类内散度尽量小，数学上等价于最大化瑞利商，其极大值由 S_w^(−1) S_b 的最大特征值给出。实现代码如下：
class LDA: def __init__(self, n_components=2): self.n_components = n_components def fit(self, X, y): classes = np.unique(y) mu = X.mean(axis=0) Sw = np.zeros((X.shape[1], X.shape[1])) Sb = np.zeros((X.shape[1], X.shape[1])) for c in classes: Xc = X[y == c] muc = Xc.mean(axis=0) # 类内散度：各类样本到类中心的离差平方和 Sw += (Xc - muc).T @ (Xc - muc) n_c = len(Xc) diff = (muc - mu).reshape(-1, 1) # 类间散度：类中心到全局中心，按样本数加权 Sb += n_c * diff @ diff.T # 广义特征值：Sb w = lambda Sw w，用伪逆转成普通特征问题 eigvals, eigvecs = np.linalg.eigh(np.linalg.pinv(Sw) @ Sb) idx = np.argsort(eigvals)[::-1] self.components_ = eigvecs[:, idx[:self.n_components]] return self def transform(self, X): return X @ self.components_代码中 Sw 累加每个类别的 (Xc − muc)^T (Xc − muc)，它本质上是每个类别协方差矩阵的按样本数加权和。Sb 的构造用到了每个类中心 muc 与全局中心 mu 的差的外积，再乘上该类的样本数 n_c。这里用 np.linalg.pinv(Sw) 而不是直接求逆，目的和上一节一致：Sw 在特征维度较高时可能不可逆，伪逆让程序至少能跑通。
但 LDA 有一个结构性的限制：类别数为 C 时，Sb 的秩最多是 C−1，因此 S_w^(−1) S_b 的非零特征值也最多只有 C−1 个。对鸢尾花数据集，C=3，LDA 最多给出 2 个有效投影方向。如果强行设置 n_components=3，第三个方向对应的是接近零的特征值，方向由数值噪声决定，投影结果没有任何统计意义。这也是后续分析里 LDA 降维只取 1 维或 2 维的原因——不是不想保留更多信息，而是方法本身不允许。相比之下 PCA 可以降到 1 维、2 维或 3 维，于是刀切法就有了对比的空间：同样降到 1 维、2 维时，PCA 和 LDA 哪个误判率更低。
4. 高斯混合模型与极大似然：软聚类如何把方差和相关性织进密度
4.1 从 K 均值到 GMM：硬聚类和软聚类的本质差别
k-均值是最常见的聚类方法，它把每个样本硬性归到距离最近的中心，更新中心时只取该类内样本的均值。这个做法的局限在于它只关心中心位置，完全忽略簇的形状。当不同类别在各个方向上的方差不一样、维度之间存在相关时，k-均值会把旋转椭球形的簇当成球形来切分，边界自然不准。
高斯混合模型（GMM）解决了这个问题。它假设数据由 K 个高斯成分混合生成，密度函数写为 p(x) = Σπ_k N(x|μ_k, Σ_k)，其中 π_k 是第 k 个成分的混合权重、满足 Σπ_k = 1，N(x|μ_k, Σ_k) 是均值为 μ_k、协方差矩阵为 Σ_k 的多元高斯密度。每个成分不仅学自己的中心，还学一个完整的协方差矩阵，因此可以表达不同方向上展布不同、维度间有相关的椭圆簇。在鸢尾花数据上，4 个属性是连续数值、近似正态，用高斯成分描述每个类别非常自然。
与 k-均值最本质的差别是软分类：GMM 不给每个样本只贴一个标签，而是计算样本属于每个成分的后验概率，称为责任。当数据点落在两个簇的重叠区域时，它同时以一定概率归属两边。这种软输出不仅在聚类上更贴近真实分布，也为后面的判别提供了后验概率，可以和马氏距离无缝衔接。
4.2 EM 的 E 步与 M 步：参数更新的底层机制
GMM 的参数没有解析解，因为每个样本究竟由哪个成分生成是隐变量，标准做法是极大似然意义下的 EM 迭代。每次迭代分两步：E 步在固定参数时计算每个样本属于每个成分的责任；M 步用责任作为权重，重新估计 π、μ、Σ。具体代码形式如下：
def gmm_em(X, K=3, max_iter=200, tol=1e-4, seed=42): rng = np.random.default_rng(seed) n, d = X.shape pi = np.ones(K) / K # 随机选 K 个样本作为初始中心；更稳的做法是用 k-means 初始化 mu = X[rng.choice(n, K, replace=False)] cov = np.tile(np.eye(d), (K, 1, 1)) + 1e-6 * np.eye(d) for _ in range(max_iter): resp = np.zeros((n, K)) # E 步：计算责任，即第 k 个成分对样本 x_i 的后验概率 for k in range(K): resp[:, k] = pi[k] * multivariate_normal_pdf(X, mu[k], cov[k]) resp /= resp.sum(axis=1, keepdims=True) # M 步：按责任加权更新参数 Nk = resp.sum(axis=0) # 每个成分的等效样本数 pi = Nk / n for k in range(K): mu[k] = (resp[:, k:k+1] * X).sum(axis=0) / Nk[k] diff = X - mu[k] cov[k] = (resp[:, k:k+1] * diff).T @ diff / Nk[k] cov[k] += 1e-6 * np.eye(d) # 正则，防止协方差奇异 return pi, mu, covE 步的关键在责任矩阵 resp 的归一化：先按每个成分的似然加权算出原始值，再对每行求和归一，保证每个样本对各个成分的归属概率加起来等于 1。M 步里 Nk 是第 k 个成分累积的责任，直观理解就是「有多少个等效样本被这个成分认领」。μ 的更新是责任加权的均值，Σ 的更新是责任加权的外积和除以 Nk，相当于对协方差矩阵做加权最大似然估计。收敛判据通常用相邻两次迭代的对数似然增量，小于 tol 就停止。
注意：协方差矩阵加 1e-6 的对角项不是可有可无的细节。当某个成分的等效样本数很少，或者数据在高维空间里稀疏时，Σ 会趋向半正定甚至奇异，这时候马氏距离里的 Σ^(−1) 会出现极大值，分类结果瞬间失真。正则项相当于给协方差设了一个下限，让数值计算保持稳定。
4.3 降维、GMM、马氏距离如何构成一条完整预测管线
整个分析流程是：原始 4 维数据先经过 PCA 或 LDA 降到 1～2 维，然后用 GMM 在这个低维空间拟合 3 个高斯成分，每个成分对应一个类别。对新样本做预测时，把它投影到同样的低维空间，分别计算它到 3 个成分中心的马氏距离，距离最近的那个成分就是预测类别。这个判断方式和直接用后验概率 p(k|x) 等价，因为多元高斯密度的对数恰恰是负的马氏距离平方的一半，再减去一个与方向有关的常数项；对同一个样本比较不同成分时，常数项虽有差异，但主要排序由马氏距离决定。
降维在这里不是可有可无的预处理，而是直接影响 GMM 的训练质量。4 维时每个高斯成分的协方差矩阵是 4×4，包含 4 个方差和 6 个协方差参数，3 个成分一共 30 个参数，要靠 150 条样本估计，方差很大。降到 2 维后每个协方差矩阵只有 5 个独立参数，3 个成分共 15 个参数，估计稳定性明显提升。同时，降维丢弃的维度往往对应噪声和冗余信息，GMM 在这些方向上不必再拟合无意义的分布。这就是「降维同时去噪」在概率建模层面的真实含义。
5. 刀切法与避坑手册：留一法如何拆穿降维的自我感觉良好
5.1 刀切法的机制：每个样本都有一次「没见过」的机会
刀切法在这里的本质就是留一法交叉验证。具体步骤是：从某一类 Gi 的 n_i 个训练样本中剔除一个，用剩下 n_i−1 个样本重新建立判别规则，然后用这条规则对剔除的那一个样本判类；如此重复，直到该类每个样本都被判过一次，记录错判个数 k_i；所有类别执行相同操作后，总误判率等于总错判数除以总样本数。
为什么不用回报法？回报法用全部样本训练、再对全部样本测试，训练集和测试集完全重叠，相当于「考自己出过的题」。模型只是记住了数据，误判率低得没有意义。刀切法让每个样本都经历一次「模型没在训练时见过我」的考验，得到的误判率是对泛化能力更合理的估计。它相比 K 折交叉验证的差别在于每次只留一个样本，训练集与完整数据集只差 1 条，模型学到的信息损失最小，适合 150 条样本这样的小数据集；代价是运算量相对大，因为整个「降维 + GMM + 距离判类」流程要重新执行 150 次。
5.2 刀切法循环里最容易被忽视的顺序问题
刀切法代码结构看似简单，但顺序问题会直接毁掉结论。正确做法是把降维、GMM 拟合全部放进循环里，每次只用训练集的样本完成这些拟合，测试样本只参与最后的距离计算：
def jackknife_evaluate(fit_predict, X, y): n = len(y) errors = 0 for i in range(n): train_idx = np.delete(np.arange(n), i) test_idx = i pred = fit_predict( X[train_idx], y[train_idx], X[test_idx:test_idx+1] ) errors += (pred[0] != y[i]) return errors / nfit_predict 内部依次完成降维拟合、GMM 拟合、马氏距离判类，所有拟合参数都只从训练样本中得到。如果为了省时间，先在整个数据集上做 PCA，再用降维后的结果循环留一，属于信息泄漏——降维阶段已经看过测试样本的分布，测试误差会被系统性低估。刀切法比较不同降维方案时，必须保证每种方案都在同等的信息条件下被评估，否则「PCA 误判率比 LDA 低」这类结论可能是偷看数据造成的假象。
5.3 我在鸢尾花任务上踩过的五个坑
第一个坑：直接用欧氏距离判类，错判率明显偏高。现象是 Setosa 与两类边界清晰的样本没问题，但 Versicolour 与 Virginica 重叠区域的样本大量分错。原因是欧氏距离忽略特征相关性和量纲，把椭圆簇当成球簇。解决方法是改用马氏距离，并配合 GMM 学出的协方差矩阵。
第二个坑：EM 结果每次跑都不一样。现象是多次运行后，有时两个高斯成分挤在同一簇里，有时某个类别的样本没有任何成分认领。原因是随机初始化让小概率的局部最优反复出现。解决方法是固定随机种子，最好用 k-均值的结果初始化 μ，再进入 EM 迭代；同时保留协方差正则项，防止单个成分的等效样本数趋近零。
第三个坑：不降维直接用 GMM 建模，误判率反而波动大。现象是 4 维下交叉验证的结果时好时坏。原因是 4 维协方差矩阵参数多、150 条样本估计不稳，过拟合风险高。解决方法是先用 PCA 或 LDA 降到 2 维，让每个成分的参数降到 5 个。
第四个坑：LDA 降维设置 n_components=3。现象是第三个投影方向上的特征值几乎为零，降维后 GMM 在这条轴上疯狂拟合噪声。原因是三类的类间散度矩阵秩最多为 2，不可能存在第三个有效方向。解决方法是把 LDA 的降维维度上限设为类别数减一，即 2。
第五个坑：在整个数据集上先做 PCA 再做刀切法。现象是误判率低得惊人，但换成新样本后表现明显变差。原因是降维阶段的均值、标准差和特征向量偷看了所有样本的信息，测试样本的信息已经参与训练。解决方法是在刀切法循环内部重新执行完整的降维拟合，确保每个测试样本在训练阶段完全不可见。
6. 设计取舍的最后一道思考：四种方法怎么咬合成一个自洽框架
把这四种方法拆开看，各自解决不同层面的问题：马氏距离是度量，决定「怎么算相似」；PCA 和 LDA 是特征变换，决定「在哪个空间里算相似」；GMM 是密度建模，决定「用什么样的分布假设描述每个类别」；刀切法是验证策略，决定「凭什么相信前面的选择」。它们不在同一个层面上，却恰好覆盖了统计分析从描述到推断的完整路径。真正有价值的不是某一算法的公式，而是这条流程的先后关系：先根据数据相关性选度量，再用降维控制模型复杂度，然后建概率模型，最后用留一法检验泛化能力。
这个框架也有明显的边界。PCA 不考虑标签，对分类任务可能保留无用方向；LDA 假设各类协方差相近，如果类别形状差异很大，线性判别会失真；GMM 假设每个簇可被高斯刻画，数据严重偏态时拟合效果有限；刀切法在小样本上偏差小但方差大，且运算成本高。进阶路线也清晰可见：用 SVD 直接做 PCA 可以绕开协方差矩阵的条件数问题；若怀疑各类协方差不同，可以改用二次判别分析（QDA），它相当于让每个类别保留自己的协方差矩阵，这和 GMM 的思想同源；若要决定 GMM 的成分个数，可以用 BIC 对模型复杂度做权衡，而不是固定 3 个成分。把这些变体对照着读，原本孤立的公式会连成一张网。
从那以后，我每次拿到一份新的数值型数据，都会强制走一遍「看相关性 → 定度量 → 选降维 → 建概率模型 → 留一验证」这个流程，先不管算法多高级，把每个环节的假设和数据对齐，再谈调参。踩过的坑少了，结论也更能站住脚。希望帮到你。
本文还有配套的精品资源，点击获取
