鸢尾花数理统计大作业：马氏距离、GMM、PCA、LDA与刀切法全链路解析

简介这份数理统计大作业文档面向高校统计学、机器学习课程的学习者围绕经典鸢尾花数据集展开完整分析帮助读者理解多重变量数据的建模流程。内容涵盖马氏距离度量样本相似度、混合高斯模型聚类、主成分分析降维去噪、线性判别分析兼顾类间与类内差异以及刀切法评估降维方法与模型泛化能力并配有数据预处理与结果讨论可作为课程作业参考或方法对照范本。资源包共1个docx文件约511KB结构完整、便于直接查阅与二次整理。目前已有1429人学习下载适合需要系统梳理数理统计知识点、完成同类大作业或复习降维与聚类方法的读者参考借鉴。1. 一份把鸢尾花拆到骨头的数理统计大作业值不值得下鸢尾花数据集大概是每个学统计、学机器学习的人绕不开的第一课但绝大多数人拿它跑个sklearn的load_iris()、画个散点图、套个分类器看准确率就收工了。这份数理统计大作业不一样它把鸢尾花当成一个完整的分析对象从马氏距离、高斯混合模型、主成分分析、线性判别分析一路推到刀切法留一法做模型比对是一条从距离度量到降维、再到概率建模、最后用重采样方法验证的完整链路。适合正在做课程设计、需要一份能照着复现的统计建模范例的读者也适合想把 PCA 和 LDA 的区别真正讲清楚、而不是背结论的从业者。它解决的核心问题是面对一组有标签的高维数值数据怎么选降维方式、怎么建模、怎么用一套站得住脚的评估方法把不同方案比出高下。2. 先把数据和方法论立住鸢尾花的四个属性到底怎么用2.1 数据集结构与符号约定鸢尾花数据集共 150 个样本3 个类别各 50 个每个样本 4 个属性花萼长度、花萼宽度、花瓣长度、花瓣宽度单位都是厘米。类别标签是 Setosa、Versicolour、Virginica 三类。这份作业在符号上做了统一约定输入向量记作 x样本数 n四个属性分别对应花萼长度、花萼宽度、花瓣长度、花瓣宽度因变量用 0/1/2 编码三个类别样本均值记 μ协方差矩阵记 Σ相关矩阵记 R线性变换方向记 u。这套符号看着琐碎但后面推导马氏距离、PCA 特征分解、LDA 的类间类内散度时全靠它符号不统一公式推到一半就会乱。数据来源是经典的公开数据集作业里明确写了从标准机器学习数据仓库获取150 条样本、4 个数值型属性、3 个类别这个规模小到可以在笔记本上几秒钟跑完又足够体现多重变量分析的全部要素。四个属性之间存在明显的线性相关尤其是花瓣长度和花瓣宽度、花萼长度和花瓣长度之间这正是后面能做降维的前提——如果四个属性彼此独立PCA 降维就没有意义了。2.2 为什么选这四个方法而不是直接上分类器很多人第一反应是鸢尾花这么简单直接逻辑回归或者随机森林不就完了这份作业的价值恰恰在于它不走这条路。它要回答的是统计建模层面的问题数据的分布假设成不成立、降维到什么程度合适、有标签和无标签的降维差多少、评估方法本身可不可靠。选马氏距离是因为欧氏距离在四个量纲相同但方差不同的属性上会失真马氏距离用协方差矩阵做白化把属性间的相关性和尺度差异都吸收进去是后面用高斯均值做判别的距离基础。选高斯混合模型是因为鸢尾花的四个属性基本符合高斯假设用多个高斯成分去拟合每个类别的分布比硬套一个球形聚类比如 k-均值更贴合真实数据的椭圆分布。选 PCA 和 LDA 做对比是因为一个无监督、一个有监督正好能体现“用不用标签”对降维结果的影响。选刀切法是因为样本只有 150 条留一法能把每个样本都当一次测试集评估结果比单次划分稳定得多。提示这套方法组合不是随便凑的每一步的输出都是下一步的输入——马氏距离依赖协方差估计协方差估计依赖降维后的数据降维效果又靠刀切法来评判链条是闭合的。2.3 数据预处理与标准化的必要性四个属性虽然都是厘米量纲一致但方差差异不小花瓣长度的方差明显大于花萼宽度。PCA 如果直接在协方差矩阵上做特征分解方差大的属性会主导主成分方向这不是我们想要的。作业里的处理方式是用相关矩阵代替协方差矩阵等价于先对每个属性做标准化再分解这样能消除量纲和方差量级的影响。具体操作上标准化就是每个属性减去均值、除以标准差import numpy as np def standardize(X): # 按列属性计算均值和标准差 mu X.mean(axis0) sigma X.std(axis0, ddof1) # 样本标准差无偏估计 # 避免除零标准差为 0 的属性直接置 1 sigma[sigma 0] 1.0 return (X - mu) / sigma, mu, sigma这段代码的关键在ddof1样本标准差要用无偏估计和后面协方差矩阵的估计口径保持一致否则相关矩阵的对角线不严格等于 1特征值分解出来的主成分方差贡献率会对不上。sigma 0的保护是防御性写法鸢尾花数据里不会触发但换成别的数据集就可能踩到。标准化之后的矩阵再算相关矩阵 RR 的对角线全是 1非对角线是属性间的皮尔逊相关系数这一步是后面 PCA 和 LDA 共同的输入基础。3. 四个核心算法的实现从马氏距离到刀切法3.1 马氏距离为什么它比欧氏距离更适合做判别马氏距离的定义是对均值 μ、协方差矩阵 Σ 的多变量向量样本 x 到分布中心的距离为 sqrt((x-μ)ᵀ Σ⁻¹ (x-μ))。它有三个关键性质当 Σ 是单位矩阵时退化成欧氏距离本质上是标准化后的欧氏距离不受量纲变化影响是无量纲量。这三点决定了它在鸢尾花问题里的地位——四个属性方差不同、彼此相关欧氏距离会把方差大的方向算得过重马氏距离用 Σ⁻¹ 把数据“拉回”各向同性的空间再量距离。实现上要注意协方差矩阵可能接近奇异尤其是降维到一维或两维、样本又少的时候直接求逆会数值不稳。常见做法是加一个极小的对角正则项或者用伪逆def mahalanobis(x, mu, cov, reg1e-6): # cov 为协方差矩阵reg 为正则项防止奇异 cov_reg cov reg * np.eye(cov.shape[0]) diff x - mu inv_cov np.linalg.inv(cov_reg) return np.sqrt(diff.T inv_cov diff)reg这个参数看着不起眼但它是血泪经验——降维到一维时协方差矩阵是 1×1样本少的时候方差估计可能极小不加正则项求逆会炸出天文数字判别结果直接翻车。正则项取 1e-6 到 1e-4 之间一般够用太大又会把马氏距离拉回欧氏距离失去意义。判别规则就是算测试样本到每个类别高斯均值的马氏距离取最小者所属类别。3.2 高斯混合模型极大似然估计下的参数求解高斯混合模型把每个类别的分布看成若干个高斯成分的加权和概率密度是各成分密度乘以混合权重的求和。作业里的思路是对降维后的数据在每个类别内假设高斯分布用极大似然估计求均值和协方差然后对新样本算它落在各类别分布下的似然取最大者。极大似然估计下均值的估计就是样本均值协方差的估计是样本协方差除以 n 还是 n-1 取决于是否有偏作业里用无偏估计。如果类别内成分数大于 1就需要 EM 算法迭代求混合权重、均值和协方差。鸢尾花每个类别内部相对集中单成分高斯已经能拟合得不错但代码结构上留出成分数的接口更通用from sklearn.mixture import GaussianMixture def fit_gmm(X, n_components1, reg_covar1e-6): # reg_covar 加到协方差对角线防止某成分方差塌缩 gmm GaussianMixture( n_componentsn_components, covariance_typefull, # 允许各方向方差不同且相关 reg_covarreg_covar, random_state42 ) gmm.fit(X) return gmmcovariance_typefull是关键选择它允许每个高斯成分有完整的协方差矩阵能刻画不同方向方差不同、维度间相关的情况这正是作业里强调的“比 k-均值更好”的地方。reg_covar是另一个后悔药某个成分如果只分到一两个样本协方差会接近奇异不加正则项 EM 会发散。random_state固定是为了结果可复现EM 对初值敏感不固定种子每次跑出来的聚类标签可能不一样。3.3 主成分分析相关矩阵特征分解与降维程度选择PCA 的数学核心是对相关矩阵或协方差矩阵做特征分解特征值从大到小排序取前 k 个特征值对应的单位特征向量作为投影方向把原数据投影到这 k 个方向上得到主成分。特征值的大小代表该主成分承载的方差也就是信息量。作业里设计了一个 PCA 类fit函数以训练样本为输入求相关矩阵、排序特征值、选最大的 n_components 个特征向量transform做投影。def pca_fit(X, n_components): # X 已标准化相关矩阵等于协方差矩阵 R np.corrcoef(X, rowvarFalse) eigvals, eigvecs np.linalg.eigh(R) # 对称矩阵用 eigh 更稳 # eigh 返回升序反转成降序 idx np.argsort(eigvals)[::-1] eigvals eigvals[idx] eigvecs eigvecs[:, idx] # 取前 n_components 个方向 components eigvecs[:, :n_components] explained eigvals[:n_components] / eigvals.sum() return components, explained def pca_transform(X, components): return X components用np.linalg.eigh而不是eig是因为相关矩阵是对称正定的eigh专门处理对称矩阵数值更稳定、返回的特征向量正交性更好。explained是各主成分的方差贡献率选几个主成分就看累计贡献率到多少鸢尾花通常前两个主成分就能解释 95% 以上的方差。降维程度不是拍脑袋定的作业里正是用刀切法去比较不同降维维度下的判别效果把“降几维”变成一个可以量化评估的选择而不是凭经验。3.4 线性判别分析用标签信息放大类间、压缩类内LDA 和 PCA 最大的区别是它用了标签。它的目标是找一个投影方向 u使得投影后同类样本尽量聚拢、异类样本尽量分开。数学上就是最大化类间散度与类内散度的比值解出来是类内散度矩阵逆乘以类间散度矩阵的特征值问题取最大特征值对应的特征向量。二分类时这个方向就是 Fisher 判别方向多分类时取前 C-1 个特征向量C 是类别数鸢尾花三分类最多降到两维。def lda_fit(X, y): classes np.unique(y) n_features X.shape[1] mean_overall X.mean(axis0) S_w np.zeros((n_features, n_features)) # 类内散度 S_b np.zeros((n_features, n_features)) # 类间散度 for c in classes: Xc X[y c] mean_c Xc.mean(axis0) # 类内各样本到本类均值的散布 S_w (Xc - mean_c).T (Xc - mean_c) # 类间类均值到总均值的散布乘样本数加权 diff (mean_c - mean_overall).reshape(-1, 1) S_b len(Xc) * (diff diff.T) # 解 S_w^{-1} S_b 的特征值问题 eigvals, eigvecs np.linalg.eig(np.linalg.inv(S_w) S_b) idx np.argsort(eigvals.real)[::-1] return eigvecs[:, idx].real这里S_w累加的是每个类内样本相对本类均值的散布S_b是类均值相对总均值的散布并按类样本数加权两者相除就是 Fisher 准则。np.linalg.eig返回可能是复数取.real是因为理论上这个广义特征值问题解是实的数值误差会带出极小虚部。LDA 在有标签时通常比 PCA 判别效果好因为它直接把“分类”这个目标编码进了投影方向而 PCA 只关心方差最大可能把判别信息丰富的低方差方向丢掉——这是作业里反复强调的一点也是很多人做特征降维时最容易忽略的边界。3.5 刀切法留一法怎么做模型比对刀切法留一法的流程是从某类 n 个训练样本中剔除一个用剩下 n-1 个建立判别规则再对剔除的那个样本做判别重复直到每个样本都被判别一次统计判错个数对所有类别做完后算总判错率。它的优点是每个样本都当过一次测试集评估结果几乎无偏缺点是运算量随样本数线性增长150 个样本就要训练 150 次。def loo_accuracy(X, y, fit_predict): n len(y) errors 0 for i in range(n): # 留出第 i 个样本 mask np.ones(n, dtypebool) mask[i] False X_train, y_train X[mask], y[mask] X_test, y_test X[i:i1], y[i] pred fit_predict(X_train, y_train, X_test) if pred ! y_test: errors 1 return 1 - errors / nfit_predict是一个回调把“训练预测”封装进去这样同一套留一循环可以复用于 PCAGMM、LDAGMM 等不同组合比对时只换回调不换评估框架。这个设计是整份作业里最值得抄的部分——评估逻辑和模型逻辑解耦换降维方法、换维度数评估代码一行不用改。作业里 4.4 节的算法比对就是靠这个框架把 PCA 降到不同维度和 LDA 降到不同维度的判错率放在一起看得出降维程度和方法差异的结论。4. 复现时最容易翻车的几个地方4.1 现象PCA 和 LDA 降维后判别准确率忽高忽低原因标准化没做或者做错了。PCA 如果直接在原始协方差矩阵上分解方差大的属性主导方向LDA 如果没标准化类内散度矩阵的条件数会很差求逆不稳定。解决所有降维方法统一在标准化后的数据上做且训练集和测试集的均值、标准差必须用训练集算不能全量数据一起算否则留一法里就发生了信息泄漏。4.2 现象高斯混合模型训练报协方差奇异或 EM 不收敛原因某个高斯成分分到的样本太少协方差矩阵接近奇异或者初值选得差导致 EM 卡在局部极值。解决加reg_covar正则项1e-6 起步固定random_state保证可复现成分数从 1 开始试别一上来就设 3 个成分——鸢尾花每个类别内部单成分就够成分数多了反而过拟合。4.3 现象马氏距离算出来是 NaN 或极大值原因协方差矩阵求逆时接近奇异尤其是降维到一维、样本又少的情况。解决加对角正则项再求逆正则项量级要试1e-6 到 1e-4 之间太大马氏距离退化成欧氏距离太小压不住奇异。另一个检查点是降维后的数据维度不能超过样本数否则协方差矩阵秩亏。4.4 现象留一法跑得极慢或者结果和单次划分差很多原因留一法本身就是 n 次训练如果每次训练里还嵌套了交叉验证或者 EM 多次重启复杂度会爆炸结果差很多往往是因为单次划分的测试集太小评估本身不可靠。解决留一循环里不要嵌套重型调参超参数提前定好如果确实慢可以改用 k 折交叉验证做粗筛再用留一法在候选方案上精评作业里选留一法正是因为 150 个样本规模扛得住。4.5 现象LDA 投影维度设成 4 维报错或结果异常原因LDA 最多只能降到 C-1 维C 是类别数三分类最多两维设成 4 维时类间散度矩阵秩不够特征向量里混入零空间方向。解决投影维度上限卡在 min(n_features, C-1)鸢尾花就是 2超过这个数的维度没有判别信息纯属噪声。5. 把评估框架用活一个可复用的降维比对技巧这份作业最值得带走的不是某个算法实现而是那套“降维方法 × 降维维度 × 留一评估”的比对框架。我一般会把它改成一个二维扫描横轴是降维维度PCA 从 1 到 4LDA 从 1 到 2纵轴是留一判错率把两条曲线画在同一张图上一眼就能看出 PCA 降到几维追平 LDA、LDA 在低维时优势有多大。这个扫描的代码骨架如下import matplotlib.pyplot as plt results {PCA: [], LDA: []} for k in range(1, 5): acc_pca loo_accuracy(X_std, y, make_pca_gmm_predictor(k)) results[PCA].append(acc_pca) for k in range(1, 3): # LDA 最多 C-12 维 acc_lda loo_accuracy(X_std, y, make_lda_gmm_predictor(k)) results[LDA].append(acc_lda) for method, accs in results.items(): plt.plot(range(1, len(accs)1), [1-a for a in accs], markero, labelmethod) plt.xlabel(降维维度) plt.ylabel(留一判错率) plt.legend() plt.show()make_pca_gmm_predictor(k)和make_lda_gmm_predictor(k)是两个工厂函数返回符合fit_predict签名的回调内部先降维再拟合 GMM 再按马氏距离判别。这样换维度只改 k换方法只换工厂函数评估循环完全不动。跑出来的典型形态是PCA 一维判错率明显偏高二维迅速下降三维四维基本持平LDA 一维就已经很低二维略降。这个形态本身就是结论——有标签时一维 LDA 就能接近无标签 PCA 两维以上的效果和作业里“LDA 兼顾类间类内、效果更好”的判断对得上。几个实操上的边界要留意。第一降维维度不能超过样本数留一法里训练集是 149 个样本四维完全没问题但换成小样本数据集就要卡这个上限。第二GMM 的成分数在扫描时固定住别每个维度都重新调否则比对的是“调参后的最好结果”而不是“方法本身的能力”。第三判错率曲线要配合混淆矩阵看鸢尾花里 Setosa 和另外两类的分离度极高判错基本都发生在 Versicolour 和 Virginica 之间只看总判错率会掩盖这个结构。从那以后我每次做降维相关的比对都强制先把评估框架和模型逻辑拆开、把标准化口径统一、把留一循环里的信息泄漏堵死再去看曲线。这套习惯比记住任何一个公式都管用因为翻车往往不在算法本身而在评估链条上某个不起眼的预处理步骤。希望这份拆解能帮到你把鸢尾花这组老数据真正用出统计建模的味道。本文还有配套的精品资源点击获取