简介这是一份围绕鸢尾花数据集展开的数理统计课程大作业面向高校统计、数据科学专业学生也适合希望了解马氏距离、高斯混合模型、主成分分析、线性判别分析等经典方法实际应用的读者。文档按引言、数据及符号说明、算法原理、数据处理分析、刀切法算法比对等章节组织系统展示从数据预处理、特征降维、聚类分类到模型评估的完整流程其中马氏距离用于度量样本相似度高斯混合模型实现聚类主成分分析与线性判别分析分别进行无监督与有监督降维并通过刀切法比较不同降维方法下的分类性能。资源包仅含一份Word文档大小约五百一十一KB内容包含摘要、目录、正文与结果展示。目前已有一千四百二十六人学习下载是一份内容具体、可复用的数理统计项目范本。读者可据此理解各统计模型在真实数据上的适用场景、参数估计与模型评估思路直接参考大作业的写作结构与实验设计适用于课程报告、答辩展示或复试复习。1. 鸢尾花数据集大作业一条能从头跑到尾的降维分类链路做数理统计课程大作业时很多同学的第一反应是拿鸢尾花数据集跑一遍 PCA画个二维散点图就交差。但这套大作业解法显然不是这个思路它把马氏距离、高斯混合模型、PCA、LDA 和刀切法串成了一条完整的实验链路从数据预处理、降维、聚类建模到留一法评估每一步都有可验证的中间结果。对于正在选课或者想把这套经典数据做出层次感的读者来说这份作业的价值在于它提供了「四种方法如何配合」的完整范式——适合当模板也适合当复现练习。下面按我拆这份作业时的顺序把关键实现和踩过的坑都铺开来讲。2. 算法选型与原理马氏距离、GMM、PCA 和 LDA 该怎么搭配2.1 马氏距离尺度无关的相似度度量先说马氏距离。欧氏距离在鸢尾花这种多变量数据上有个明显的毛病它把每个特征同等对待但四个特征的取值范围差异很大。花萼长度在 4.3 到 7.9 之间花瓣宽度只有 0.1 到 2.5。如果直接用欧氏距离波动范围大的特征会主导距离计算波动范围小的特征几乎不起作用。马氏距离通过协方差矩阵把特征之间的相关性和尺度差异一起考虑进去定义是对于一个均值为 μ、协方差为 Σ 的分布样本 x 到 μ 的马氏距离是 (x-μ)ᵀΣ⁻¹(x-μ) 再开根号。这个式子有两个关键性质。第一当 Σ 是单位矩阵时马氏距离退化为欧氏距离说明欧氏距离只是马氏距离的一个特例。第二马氏距离是尺度无关的也就是说你改变特征的单位马氏距离不变这让它天然适合鸢尾花这种原始特征和标准化特征并存的情况。在实际分类任务里我的习惯是计算样本到每个类别中心的马氏距离距离最小的那个类别作为预测结果。注意这里用的是每个类别自己的协方差矩阵不是全体样本的协方差。2.2 高斯混合模型用多个正态分布描述数据高斯混合模型的核心假设是整个数据分布可以由若干个高斯成分加权叠加而成。对鸢尾花来说每个类别假设对应一个高斯成分那么一个样本属于某个类别的概率就是这个成分的概率密度乘上该成分的先验权重。模型参数包括每个成分的均值向量、协方差矩阵和混合系数求解方式是极大似然估计。手工推导极大似然时你会发现高斯分布的参数估计公式很直观均值是样本均值协方差是样本协方差。但混合模型里每个样本不知道属于哪个成分所以不能直接套公式。常见做法是用 EM 算法迭代求解先根据当前参数计算每个样本属于每个成分的后验概率E 步再根据这些后验概率加权更新均值、协方差和混合系数M 步重复直到对数似然收敛。这块作业里没有手写 EM而是直接调用库实现我复现时也比较推荐先调库把流程跑通再回头手写 EM 验证。2.3 PCA无监督降维如何工作PCA 做的事情是把一组可能存在相关性的变量通过正交变换转换成一组线性不相关的变量这些新变量叫主成分。它不利用类别标签所以是一种无监督方法。目标很直接找到一组投影方向使得数据投影后的方差最大化因为方差大意味着信息保留多。求法也不复杂计算协方差矩阵或相关矩阵做特征值分解把特征值从大到小排列取前 k 个特征值对应的特征向量组成投影矩阵。特征值的大小代表该主成分解释的方差量特征值之和占全体特征值之和的比例就是累积方差贡献率。对鸢尾花这种四维数据前两个主成分通常能解释超过 95% 的方差这就是降维到二维的理论依据。PCA 的降维同时也是一个去噪过程去掉的主成分往往是噪声或线性冗余较大的方向。2.4 LDA有监督降维与 PCA 的关键差异LDA 是 Fisher 提出的一种线性判别方法和 PCA 最大的区别是它使用了类别标签。它的目标是找到一个投影方向 u让同类样本投影后尽可能聚集异类样本投影后尽可能分离。数学上表达为最大化类间散度与类内散度的比值这个优化问题的解是 Sw⁻¹Sb 的特征向量其中 Sw 是类内散度矩阵Sb 是类间散度矩阵。这里有个值得一提的约束LDA 最多能降到类别数减一的维度。鸢尾花有三个类别所以 LDA 最多降到二维。而 PCA 最多能降到特征维度四维。这是两者的一个本质差异——LDA 的降维上限由标签类别数决定不是由特征数决定。理论上 LDA 利用了标签信息效果应该好于 PCA但这要建立在数据近似服从高斯分布、各类协方差差异不大的前提下鸢尾花数据基本满足这个条件。2.5 刀切法为什么用它评估而不直接划分测试集刀切法也叫留一法每次从训练样本中剔除一个样本用剩下的样本建立模型再对剔除的那个样本做预测循环直到每个样本都被预测过一次最后统计判错个数除以总体本数得到判错率。对 150 个样本的鸢尾花数据这相当于做 150 次训练和预测。为什么不用简单的随机划分训练集和测试集因为鸢尾花数据集本身很小随机划分一次的结果受划分方式影响很大评估结果不稳定。刀切法的好处是每个样本都被当作测试样本一次评估更充分、更接近模型的真实泛化能力。代价是计算量比较大150 次循环每次都要重新做降维和建模。不过在现代计算机上这个规模也就是几秒到几十秒的事。3. 数据准备与预处理读取、标准化与分布检查3.1 数据读取与列顺序检查鸢尾花数据最常见的格式是 CSV每行一个样本包含花萼长度、花萼宽度、花瓣长度、花瓣宽度四个特征最后一列是类别标签。有些版本在第一列还有序号读取的时候要特别注意。我拿到这份作业配套的数据表时发现表格里列名顺序看起来是乱的序号、特征、类别的排列并不一致但代码里必须有一个统一的约定四列特征、一列标签、一行一个样本。import numpy as np import pandas as pd df pd.read_csv(iris.csv, headerNone) print(df.shape) # 期望是 (150, 5) print(df.head()) X df.iloc[:, :4].values.astype(float) # 前四列是特征 y df.iloc[:, 4].values # 最后一列是标签 # 检查类别分布 print(pd.Series(y).value_counts())这段代码做的事情很基础但很重要。headerNone是因为鸢尾花的原始 CSV 往往没有列名df.iloc[:, :4]取前四列作为特征矩阵df.iloc[:, 4]取最后一列作为标签。类别分布检查是为了确认三类样本各 50 条如果数量不均衡后面的 LDA 和 GMM 结果都会受影响。我一般先跑一遍这个检查确认数据格式没问题再进入下一步。3.2 标准化与相关矩阵两种等效的消除量纲手段鸢尾花四个特征的量纲都是厘米但取值范围差异明显。花萼长度几乎跨越 4 厘米以上花瓣宽度跨度只有 2 厘米多一点。直接用原始数据做 PCA方差大的特征会天然获得更高的主成分权重这并不代表它真的更重要。两种常见的解决办法一是做 z-score 标准化让每个特征均值为 0、方差为 1二是直接计算相关矩阵代替协方差矩阵。前者是改变数据本身后者是改变参与计算的矩阵本质上效果一致。# 方法一z-score 标准化 mu X.mean(axis0) std X.std(axis0, ddof0) Xs (X - mu) / std # 方法二用相关矩阵参与 PCA R np.corrcoef(X.T)注意np.std的ddof参数样本标准差用ddof1除以 n-1总体标准差用ddof0。数理统计课程里讲极大似然估计时用的通常是总体标准差作业里一般也延续这个约定。我复现时统一用ddof0保持和 GMM 的极大似然口径一致。标准化的另一个好处是让后面的马氏距离计算数值更稳定——协方差矩阵的条件数会变小。3.3 分布检查确认高斯假设是否成立GMM 的前提是每个类别大致服从高斯分布所以建模前要确认数据是否符合这个假设。最简单的检查是画分布图看每个特征的分布是否呈现钟形。更严格的检查可以做正态性检验但作为大作业散点图和直方图已经足够说明问题。import matplotlib.pyplot as plt fig, axes plt.subplots(2, 2, figsize(10, 8)) feature_names [sepal length, sepal width, petal length, petal width] for i, ax in enumerate(axes.flat): ax.hist(X[y setosa, i], bins10, alpha0.5, labelsetosa) ax.hist(X[y versicolor, i], bins10, alpha0.5, labelversicolor) ax.hist(X[y virginica, i], bins10, alpha0.5, labelvirginica) ax.set_title(feature_names[i]) ax.legend() plt.tight_layout() plt.show()从这里能直观看到两类信息每个特征在不同类别上的分布是否大致对称以及类别之间的重叠程度。鸢尾花数据里花瓣长度和花瓣宽度对类别的区分度明显高于花萼长度和花萼宽度这为后面降维方法的选择提供了一个直觉基础——数据里本身就存在和类别强相关的方向LDA 有希望找到它。4. PCA 与 LDA 降维实战从零手写核心代码4.1 PCA特征值分解与主成分选取不用sklearn从零实现 PCA 其实也就十几行代码。核心步骤是标准化数据 → 计算相关矩阵 → 特征值分解 → 按特征值降序排列 → 取前 k 个特征向量做投影矩阵。def pca_fit(X, n_components): # X 是已经标准化后的数据 R np.corrcoef(X.T) # 相关矩阵 eigvals, eigvecs np.linalg.eigh(R) # 对称矩阵用 eigh 更稳 order np.argsort(eigvals)[::-1] # 从大到小排列 eigvals eigvals[order] eigvecs eigvecs[:, order] W eigvecs[:, :n_components] # 取前 k 个特征向量 explained_ratio eigvals[:n_components] / eigvals.sum() return W, explained_ratio W, ratio pca_fit(Xs, 2) X_pca Xs W print(explained variance ratio:, ratio) print(total:, ratio.sum())这里用np.linalg.eigh而不是np.linalg.eig因为相关矩阵是对称矩阵eigh专门利用对称性数值稳定性更好速度也更快。explained_ratio是每个主成分解释的方差占比鸢尾花数据降到二维时这个占比通常能到 0.95 以上意味着丢掉两个维度只损失不到 5% 的方差。X_pca是降维后的数据后面 GMM 和刀切法都用它作为输入。初次跑 PCA 时建议打印一下完整特征值而不是只看前两个。如果第三个特征值明显小于前两个说明降到二维是合理的选择如果特征值之间差距不大那说明数据在四个方向上都有信息强行降到二维会丢失较多信息。4.2 LDA类内散度矩阵与广义特征值问题LDA 的从零实现比 PCA 稍绕因为需要分别计算类内散度矩阵和类间散度矩阵然后求解Sw⁻¹Sb的特征值问题。通用做法是直接解广义特征值问题代码上可以用scipy.linalg.eigh一步到位避免显式求逆带来的数值问题。from scipy.linalg import eigh def lda_fit(X, y, n_components): classes np.unique(y) n_features X.shape[1] Sw np.zeros((n_features, n_features)) Sb np.zeros((n_features, n_features)) overall_mean X.mean(axis0) for c in classes: Xc X[y c] n_c Xc.shape[0] mu_c Xc.mean(axis0) Sw (Xc - mu_c).T (Xc - mu_c) mean_diff (mu_c - overall_mean).reshape(-1, 1) Sb n_c * mean_diff mean_diff.T eigvals, eigvecs eigh(Sb, Sw) # 广义特征值问题 order np.argsort(eigvals)[::-1] W eigvecs[:, order[:n_components]] return W W_lda lda_fit(Xs, y, 2) X_lda Xs W_ldaeigh(Sb, Sw)传入两个矩阵时就变成了广义特征值问题它在内部处理了Sw⁻¹Sb的求解同时利用了 Sb 和 Sw 的对称性。这里有个关键点LDA 的降维维度不能超过类别数减一。鸢尾花三个类别n_components最大只能设 2如果你设 3eigh会返回不满秩的结果后面投影矩阵会出现数值异常。另外要注意W_lda的列不是单位正交的。PCA 的特征向量天然正交LDA 的投影方向受类内散度矩阵的影响向量之间一般不正交这是正常现象不要用 PCA 的思维去要求 LDA。4.3 降维维度的选择与投影可视化降维维度怎么定不能只看累积方差贡献率。我见过很多同学把 PCA 降到二维之后发现分类效果不如原始四维数据就怀疑算法有问题——其实不是降维本身就是去掉部分信息如果信息被去掉的部分恰好和分类相关分类效果自然会下降。所以维度选择要和后面的分类任务绑定看降维后 GMM 或马氏距离分类的正确率曲线而不是单独看方差贡献率。import matplotlib.pyplot as plt plt.figure(figsize(12, 5)) plt.subplot(1, 2, 1) scatter plt.scatter(X_pca[:, 0], X_pca[:, 1], cy, cmapviridis, s30) plt.colorbar(scatter) plt.title(PCA 2D projection) plt.subplot(1, 2, 2) scatter plt.scatter(X_lda[:, 0], X_lda[:, 1], cy, cmapviridis, s30) plt.colorbar(scatter) plt.title(LDA 2D projection) plt.tight_layout() plt.show()从两张图对比能直观感受到 PCA 和 LDA 的差别PCA 的投影方向是让总体方差最大三个类别在投影面上可能仍然部分重叠;而 LDA 的投影方向是让类间距离最大、类内距离最小通常能看到更清晰的类别分界。特别是 Setosa 这一类在 LDA 投影下几乎能被一条直线分开而 PCA 投影下虽然有分离趋势但边界更模糊。这就是有监督信息带来的效果差异。4.4 降维后还要注意什么降维不是终点而是后续建模的输入。用降维后的数据做 GMM 时等价于只在高方差方向PCA或高判别方向LDA上建模这相当于主动滤掉了那些不影响分类的噪声维度。不过要注意降维过程本身使用的数据包含了全部样本如果后面用刀切法做评估严格的流程应该是在每一折训练集内重新做降维防止测试样本的信息通过降维矩阵泄露到训练过程中。这个问题在下一章展开讲它是整个流程里最容易翻车的地方。5. GMM 聚类与马氏距离判别避坑信息泄露、标签错位、奇异矩阵5.1 GMM 参数拟合与聚类降维之后用 GMM 对数据建模。对鸢尾花这类数据通常设置三个高斯成分对应三个类别。核心参数除了成分数还有协方差类型和初始化次数。from sklearn.mixture import GaussianMixture gmm GaussianMixture( n_components3, covariance_typefull, max_iter300, n_init10, random_state42 ) gmm.fit(X_pca) labels gmm.predict(X_pca) # 每个样本的成分归属 print(converged:, gmm.converged_) print(log likelihood:, gmm.lower_bound_)covariance_typefull表示每个成分使用完整的协方差矩阵允许不同方向上的方差不同也允许维度之间存在相关性这是 GMM 相对于 k-means 的核心优势——k-means 实际上假设每个簇在所有方向上方差相同且互不相关这对鸢尾花这种特征相关性较强的数据来说太苛刻了。n_init10表示 EM 算法用 10 组不同的初始值各跑一遍选对数似然最大的那组结果避免初始值把迭代带进局部最优。random_state42固定随机种子保证结果可复现这是所有实验开始前必须做的一件事。5.2 马氏距离判别每个类别用各自的协方差GMM 建模之后可以用两类方式做判别直接用 GMM 的后验概率或者用马氏距离判断样本离哪个类中心更近。二者原理相通但马氏距离更直观也更好解释。def mahalanobis_distance(x, mu, cov): delta x - mu inv_cov np.linalg.pinv(cov) # 用伪逆防止矩阵奇异 return np.sqrt(delta inv_cov delta) def predict_by_mahalanobis(X, means, covs): n X.shape[0] preds np.zeros(n, dtypeint) for i in range(n): dists [mahalanobis_distance(X[i], means[c], covs[c]) for c in range(3)] preds[i] np.argmin(dists) return preds means gmm.means_ # 每个成分的均值 covs gmm.covariances_ # 每个成分的协方差矩阵 preds predict_by_mahalanobis(X_pca, means, covs)这里有一个细节容易忽略计算马氏距离时每个样本要和三个类别中心各算一次距离每个距离用对应的类别协方差矩阵。有些实现偷懒用了全局协方差那算出来的东西叫标准化欧氏距离不完全是马氏距离效果也会打折扣。np.linalg.pinv是伪逆它对奇异矩阵有很好的容错普通np.linalg.inv遇到奇异或接近奇异的协方差矩阵会直接报错或给出极大值。5.3 避坑一GMM 的聚类标签和真实标签不对齐现象用gmm.predict得到的标签是 0、1、2直接和真实标签比较算正确率结果惨不忍睹可能只有 40% 上下但看可视化图发现聚类分明、效果很好。原因GMM 的聚类标签是成分编号编号本身没有语义。初始值不同编号 0 可能对应 Setosa也可能对应 Virginica。直接把聚类编号和真实类别编号对比等于用一个随机排列的标签去匹配真实标签。解决做标签对齐让每个聚类编号映射到它最可能对应的真实类别。def align_labels(pred, true): mapping {} for t in np.unique(true): mask true t pred_sub pred[mask] if len(pred_sub) 0: mapping[pred_sub[np.argmax(np.bincount(pred_sub))]] t return mapping mapping align_labels(labels, y) aligned np.array([mapping.get(l, -1) for l in labels]) accuracy (aligned y).mean() print(aligned accuracy:, accuracy)这段代码对每个真实类别找出该类别样本中占比最高的聚类编号建立映射关系。本质上是找两个集合之间最可能的一一对应。如果某个聚类编号没出现在任何真实类别中映射不到就标记为 -1避免崩掉。5.4 避坑二降维信息泄露导致刀切法结果虚高现象刀切法算出来的正确率接近 100%比正常预期高出一大截明显不合理。原因在留一循环之前先用全部 150 个样本做了 PCA 或 LDA求出了投影矩阵。然后循环中每次留出测试样本但投影矩阵是全量样本计算的测试样本的信息已经通过投影矩阵混进了训练过程。这相当于考试时把答案夹在草稿纸里正确率当然虚高。解决降维必须在每一折训练集内部重新计算。def jackknife_with_pca(X, y, n_components2): n len(y) correct 0 for i in range(n): mask np.ones(n, dtypebool) mask[i] False X_train, y_train X[mask], y[mask] X_test X[i:i1] mu X_train.mean(axis0) std X_train.std(axis0, ddof0) Xs_train (X_train - mu) / std Xs_test (X_test - mu) / std W, _ pca_fit(Xs_train, n_components) # 只用训练集计算投影矩阵 Z_train Xs_train W Z_test Xs_test W gmm GaussianMixture(n_components3, random_state42, n_init5) gmm.fit(Z_train) pred gmm.predict(Z_test)[0] # 标签对齐后比较 ... return correct / n这段代码的关键是第 7~9 行标准化的均值和标准差、PCA 的投影矩阵全部来自训练集。测试集只有标准化和投影这一步。这一步错了整个实验的结论都没有意义。你可能觉得 150 个样本的信息泄露一点点没关系但留一法评估的逻辑就是用极端方式刻画泛化误差任何一点泄露都会让结果失去参考价值。5.5 避坑三协方差矩阵奇异马氏距离算出 NaN现象计算马氏距离时提示矩阵奇异或者除法遇到无效值距离结果出现 inf 或 NaN。原因协方差矩阵是半正定的理论上不会奇异但当特征之间完全线性相关或者样本数少于特征数时计算机会因为数值精度问题报奇异。降维后如果再选择了一个和类别数相同的维度某些成分的协方差矩阵可能退化。解决用伪逆np.linalg.pinv替代普通逆或者给协方差矩阵加一个小的正则化项例如cov 1e-6 * np.eye(d)。前者更通用后者更贴近 GMM 本身的处理方式。GaussianMixture的reg_covar参数就是做这件事的默认是 1e-6如果数据量小或者维度高可以适当调大。5.6 避坑四PCA 不标准化等于白做现象PCA 第一主成分的载荷几乎全集中在某一个特征上其他特征权重接近零降维图上看不出类别区分。原因原始特征尺度差异大协方差矩阵由方差最大的那个特征主导PCA 自然优先保留它。这不是数据本身的规律而是量纲造成的假象。解决先标准化再 PCA或者直接用相关矩阵代替协方差矩阵。一句话标准化是 PCA 的前置条件不是可选项。6. 用刀切法横向对比降维方法、维度与分类效果的取舍6.1 刀切法评估脚本把整个流程串起来刀切法评估的完整思路是留一循环 → 训练集内降维 → GMM 建模 → 标签对齐 → 预测留出样本 → 统计判错数。计算判错率的公式是判错数除以总体本数。把第 5 章的jackknife_with_pca扩展一下同时支持 PCA 和 LDA 两种降维方式就能得到一组可对比的准确率。def jackknife_eval(X, y, methodpca, n_components2): n len(y) correct 0 for i in range(n): mask np.ones(n, dtypebool) mask[i] False X_train, y_train X[mask], y[mask] X_test X[i:i1] mu X_train.mean(axis0) std X_train.std(axis0, ddof0) Xs_train (X_train - mu) / std Xs_test (X_test - mu) / std if method pca: W, _ pca_fit(Xs_train, n_components) else: W lda_fit(Xs_train, y_train, n_components) Z_train Xs_train W Z_test Xs_test W gmm GaussianMixture(n_components3, random_state42, n_init5) gmm.fit(Z_train) pred_train gmm.predict(Z_train) mapping align_labels(pred_train, y_train) pred_test gmm.predict(Z_test)[0] pred_label mapping.get(pred_test, -1) correct (pred_label y_test) return correct / n这里有两点值得说明。第一GMM 的标签对齐映射是在每一折训练集上重新计算的因为每次训练集不同聚类编号和真实类别的对应关系可能变化。第二n_components对 LDA 有硬约束最大只能取 2而 PCA 可以取到 4对比实验时要分别设置合理的取值区间。6.2 典型结果与参数参考在鸢尾花数据集上几组配置的结果通常落在下面的区间降维方法降维维度留一法正确率典型区间说明不降维40.95 左右特征全保留GMM 直接拟合PCA20.90~0.94信息损失较多类别边界有重叠PCA30.95~0.97保留大部分有效信息LDA20.97~0.98利用标签信息投影方向更具判别力这个表说明两个问题一是降维程度的选取不能贪低PCA 降到二维虽然损失方差占比不到 5%但损失的这部分恰好对分类有帮助所以正确率明显低于三维二是 LDA 在同样维度下通常优于 PCA因为它的投影方向是冲着分类去的。这些数字在不同随机种子下会有小幅波动但整体趋势非常稳定。做实验报告时建议多跑几个随机种子把均值和标准差都写上。6.3 一个值得保留的进阶做法如果想让结果更细可以把刀切法里的 GMM 换成不同协方差类型跑一组对比。full允许各方向方差不同、维度间相关diag假设维度独立但方差不同spherical假设各方向方差相同。鸢尾花数据用full效果最好这个结论本身也可以写进报告。从那以后我跑这类小样本实验都强制走一遍固定随机种子、每一折重新降维、标签对齐三步检查凡是这三步有一处漏掉结果我都当无效处理。做实验最怕的不是算法不 work而是流程里有个不起眼的环节悄悄污染了结论。希望帮到你。本文还有配套的精品资源点击获取