简介这份资源是数理统计课程大作业的完整文档面向正在学习数理统计、多元统计或机器学习基础的高校学生帮助其理解如何将课堂理论落到真实数据集分析中。文档以经典鸢尾花数据集为对象串联马氏距离、混合高斯模型、主成分分析、线性判别分析、刀切法等知识点从数据预处理、算法原理分析到降维与聚类结果比对形成一条较完整的分析链路适合作为课程作业参考或复习提纲。压缩包内共1个docx文件约511KB内容涵盖摘要、符号说明、算法原理、数据处理与结论等章节结构清晰便于按模块查阅。目前已有1426人学习下载说明其在同类课程作业中具有一定参考价值。读者可从中获得各方法的适用场景说明、降维与聚类流程的组织方式以及刀切法评估模型差异的思路对完成类似数据分析作业或梳理统计方法体系均有帮助。1. 鸢尾花数据集上的降维与聚类一份数理统计大作业的完整复现路径拿到一份数理统计大作业最怕的不是公式看不懂而是代码跑出来和报告里的数字对不上。这份围绕鸢尾花数据集展开的作业把马氏距离、高斯混合模型、主成分分析、线性判别分析、刀切法串成了一条完整的分析链路先用 PCA 和 LDA 把四维特征压到低维再用 GMM 在极大似然意义下拟合分布最后用留一法交叉验证去比较不同降维方式的泛化差异。它适合正在做统计学习课程设计的学生也适合想补一遍经典多元统计方法实操的工程师——因为这条链路里的每一步在真实项目里都能找到对应场景比如特征降噪、无监督分群、模型选型评估。下面按“数据怎么进、算法怎么落、结果怎么验”的顺序拆开讲。2. 数据准备与马氏距离从四列数值到可计算的样本空间2.1 鸢尾花数据的结构与本作业的符号约定鸢尾花数据集共 150 条样本三个类别各 50 条每条样本有四个数值型属性花萼长度、花萼宽度、花瓣长度、花瓣宽度单位都是厘米。作业里把输入向量记作 x样本数记作 n类别标签记作 0、1、2 分别对应三个品种。这个数据集之所以被反复使用是因为它的四个属性之间存在明显线性相关——花瓣长度和花萼长度高度正相关这正是 PCA 能起作用的前提同时三个类别在特征空间里并非完全线性可分第二类和第三类有重叠这又给分类和聚类留出了区分度。在动手写代码前先把数据整理成两个数组特征矩阵 X 形状为 (150, 4)标签向量 y 形状为 (150,)。常见做法是用 pandas 读入后直接转 numpy注意列顺序要和报告里的符号说明一致否则后面算出来的特征向量方向会对不上。作业里明确采用 UCI 的原始版本没有做额外的缺失值处理因为这份数据本身没有缺失。import numpy as np import pandas as pd # 列名按作业符号说明的顺序花萼长度、花萼宽度、花瓣长度、花瓣宽度 col_names [sepal_length, sepal_width, petal_length, petal_width, label] df pd.read_csv(iris.data, headerNone, namescol_names) # 标签映射为 0/1/2与作业符号说明保持一致 label_map {Iris-setosa: 0, Iris-versicolor: 1, Iris-virginica: 2} df[label] df[label].map(label_map) X df[[sepal_length, sepal_width, petal_length, petal_width]].values y df[label].values print(X.shape, y.shape) # (150, 4) (150,)这段代码的关键在于列顺序和标签映射。列顺序一旦和报告里的符号定义错位后面 PCA 载荷矩阵的解释就会整体偏移标签映射则决定了 LDA 的类间散度矩阵按什么顺序累加虽然不影响最终分类正确率但会影响中间量的可读性。参数上没有什么可调的唯一要注意的是读入时 headerNone因为原始文件没有表头。2.2 马氏距离的计算与它在分类中的角色马氏距离是这份作业里第一个需要落地的公式。它和欧氏距离的区别在于考虑了协方差结构当各维方差不同、维度之间相关时欧氏距离会被量纲和相关性带偏而马氏距离通过协方差矩阵的逆做白化把样本投影到一个各向同性的空间里再量距离。作业里给出的定义是 d(x, μ) sqrt((x-μ)^T Σ^{-1} (x-μ))其中 μ 是均值向量Σ 是协方差矩阵。当 Σ 为单位阵时退化为欧氏距离这个性质在调试时很有用——可以先令 ΣI 验证代码路径再换成真实协方差。在这份作业的流程里马氏距离不是单独用来分类的而是配合 GMM 使用GMM 拟合出每个高斯成分的均值 μ_k 和协方差 Σ_k 之后用样本到各成分均值的马氏距离来判断归属。这里有个容易翻车的点协方差矩阵在维度较高或样本较少时可能接近奇异直接求逆会报错或数值不稳定常见做法是加一个很小的对角正则项或者用伪逆。def mahalanobis(x, mu, cov, reg1e-6): # 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) # 以第一类样本为例计算类内均值与协方差 X0 X[y 0] mu0 X0.mean(axis0) cov0 np.cov(X0, rowvarFalse) d mahalanobis(X[0], mu0, cov0) print(样本0到第一类均值的马氏距离:, d)参数说明reg 取 1e-6 是经验值数据标准化后协方差量级在 1 附近时够用如果数据没标准化、量纲差异大这个值要相应放大。rowvarFalse 表示每列是一个变量这是 numpy 协方差函数的默认行为容易搞反的地方写错会得到一个 150×150 的矩阵而不是 4×4。马氏距离本身对量纲不敏感但前提是协方差估计准确所以样本量太小时它反而不如欧氏距离稳。3. PCA 与 LDA 降维两种投影思路的实现与差异3.1 主成分分析相关矩阵特征分解与降维程度选择PCA 的目标是把四个相关变量转换成一组线性不相关的变量按方差从大到小排列。作业里特别提到对于不同性质的数据要先用相关矩阵代替协方差矩阵以消除量纲影响。鸢尾花四个属性单位相同理论上用协方差矩阵也行但用相关矩阵更稳妥也符合作业报告里的推导。实现上就是对相关矩阵做特征分解取最大的 k 个特征值对应的特征向量组成投影矩阵。降维程度的选择是这份作业的一个分析重点。作业用留一法去比较不同 k 值下的表现而不是拍脑袋取 k2。常见做法是先看累计方差贡献率一般取到 85% 到 95%但贡献率只反映信息保留量不反映分类或聚类效果所以最终还是要靠刀切法这类外部评估来定。这一点在报告里体现得很清楚PCA 是无监督的它保留的是方差大的方向而方差大的方向未必是区分类别最有效的方向。class PCA: def __init__(self, n_components): self.n_components n_components def fit(self, X): # 用相关矩阵等价于先标准化再求协方差 self.mean_ X.mean(axis0) self.std_ X.std(axis0, ddof1) X_std (X - self.mean_) / self.std_ R np.corrcoef(X_std, rowvarFalse) eigvals, eigvecs np.linalg.eigh(R) # eigh 返回升序取最大的几个并保持降序 idx np.argsort(eigvals)[::-1][:self.n_components] self.components_ eigvecs[:, idx] self.explained_variance_ratio_ eigvals[idx] / eigvals.sum() return self def transform(self, X): X_std (X - self.mean_) / self.std_ return X_std self.components_ pca PCA(n_components2).fit(X) Z_pca pca.transform(X) print(前两个主成分方差贡献率:, pca.explained_variance_ratio_)这段实现里ddof1 表示样本标准差和 numpy 默认的总体标准差不同写错会让标准化后的数据和报告对不上。eigh 用于对称矩阵比 eig 更稳也更快返回的特征值升序所以要用 argsort 反转。components_ 的每一列是一个主成分方向transform 时用标准化后的数据右乘顺序不能反。如果只想要降维后的坐标做可视化到这一步就够了如果要接 GMM就把 Z_pca 直接喂进去。3.2 线性判别分析类间散度与类内散度的广义特征分解LDA 和 PCA 的根本区别在于它用了标签。作业里说得很直白LDA 在投影时同时放大类间离差、缩小类内离差所以有标签时效果通常比 PCA 好很多。实现上LDA 要求解的是广义特征值问题 S_b w λ S_w w其中 S_b 是类间散度矩阵S_w 是类内散度矩阵。投影方向就是 S_w^{-1} S_b 的最大特征值对应的特征向量。对于三分类问题LDA 最多只能降到 2 维因为类间散度矩阵的秩最多是类别数减一这一点在选降维程度时是硬约束不能像 PCA 那样随意取 k。def lda_fit(X, y, n_components2): classes np.unique(y) mu_all X.mean(axis0) S_w np.zeros((X.shape[1], X.shape[1])) S_b np.zeros((X.shape[1], X.shape[1])) for c in classes: Xc X[y c] mu_c Xc.mean(axis0) S_w (Xc - mu_c).T (Xc - mu_c) n_c Xc.shape[0] diff (mu_c - mu_all).reshape(-1, 1) S_b n_c * (diff diff.T) # 加正则项避免 S_w 奇异 S_w 1e-6 * np.eye(S_w.shape[0]) eigvals, eigvecs np.linalg.eig(np.linalg.inv(S_w) S_b) idx np.argsort(eigvals.real)[::-1][:n_components] return eigvecs[:, idx].real W lda_fit(X, y, n_components2) Z_lda X W print(LDA 投影后形状:, Z_lda.shape)这里 S_w 累加的是每个类内样本对类均值的离差平方和S_b 累加的是类均值对总均值的加权离差权重是类样本数。eig 返回的可能有复数取实部是因为理论上这个广义特征值问题的解是实的数值误差带来的虚部可以丢弃。正则项同样是为了防止 S_w 奇异尤其当某个类样本数少于特征维度时。投影矩阵 W 的列是判别方向Z_lda 就是降维后的坐标可以直接接 GMM 或分类器。3.3 两种降维方式的选型边界把 PCA 和 LDA 放在一起看选型逻辑其实很清楚没有标签、或者标签不可信时用 PCA它只依赖数据本身的方差结构有可靠标签、且目标是分类或聚类时优先 LDA因为它把类别信息编码进了投影方向。但 LDA 也有边界类别数少时降维上限低三分类最多两维如果原始特征维度很高而类别很少LDA 压得会比 PCA 狠可能损失过多信息。作业里用刀切法去比较两者在不同降维程度下的表现正是为了把这个边界量化出来而不是停留在“LDA 更好”这种定性判断上。4. 高斯混合模型与刀切法聚类拟合和泛化评估4.1 GMM 的极大似然估计与 EM 迭代高斯混合模型假设数据由若干个高斯成分按不同权重混合而成每个成分有自己的均值、协方差和混合系数。作业里用极大似然估计求参数但 GMM 的似然函数直接求导没有闭式解实际落地靠 EM 算法迭代E 步算每个样本属于各成分的后验概率M 步用这些后验概率加权更新均值、协方差和混合系数。和 k-均值相比GMM 允许各成分有不同方向的协方差能刻画椭球形的簇而不是只认球形簇这对鸢尾花第二、三类那种有相关性的分布更合适。from sklearn.mixture import GaussianMixture # 在 LDA 降维后的两维数据上拟合 GMM gmm GaussianMixture(n_components3, covariance_typefull, init_paramskmeans, max_iter200, random_state42) gmm.fit(Z_lda) labels_gmm gmm.predict(Z_lda) print(各成分权重:, gmm.weights_) print(收敛与否:, gmm.converged_)参数上n_components 取 3 是因为已知三个类别但实际聚类时这个数往往未知作业里是靠刀切法或 BIC 去辅助判断。covariance_typefull 表示每个成分有完整的协方差矩阵这是 GMM 相对 k-均值优势的来源但参数量也最大样本少时容易过拟合可以退到 diag 或 tied。init_paramskmeans 用 k-均值初始化均值比随机初始化稳random_state 固定后结果可复现。converged_ 为 False 时要调大 max_iter 或检查数据尺度。4.2 刀切法留一法的评估流程刀切法在这份作业里承担的是模型比较的角色。它的流程是每次从训练集里剔除一个样本用剩下的样本建立判别规则再对剔除的那个样本做预测重复到每个样本都被预测一次最后统计判错率。作业里特别指出刀切法比回报法更合理因为回报法用训练集自身评估会高估性能而刀切法每个样本都是在没见过的模型上被预测的。代价是运算量大n 个样本要训练 n 次但 150 条样本的量级在普通机器上完全跑得动。def loo_accuracy(X, y, fit_predict): n X.shape[0] correct 0 for i in range(n): mask np.ones(n, dtypebool) mask[i] False pred fit_predict(X[mask], y[mask], X[i]) if pred y[i]: correct 1 return correct / n def pca_gmm_predict(X_train, y_train, x_test, k2): pca PCA(n_componentsk).fit(X_train) Z_train pca.transform(X_train) Z_test pca.transform(x_test.reshape(1, -1)) gmm GaussianMixture(n_components3, covariance_typefull, random_state42).fit(Z_train) return gmm.predict(Z_test)[0] acc loo_accuracy(X, y, pca_gmm_predict) print(PCAGMM 留一法正确率:, acc)这段代码把降维和聚类包成一个 fit_predict 函数刀切法只负责切分和统计这样换降维方法时只改传入的函数即可比较逻辑不变。注意每次循环都要重新 fit PCA 和 GMM不能复用全量数据上拟合的模型否则就泄漏了测试样本的信息正确率会虚高。这也是刀切法运算量大的原因但换来的是更可信的评估。4.3 降维程度与降维方法的比对结果怎么读作业的核心分析落在 4.4 节的刀切法算法比对。读这类结果时重点不是看哪个数字最高而是看趋势PCA 随着 k 从 1 增到 4正确率怎么变LDA 在 k1 和 k2 时差多少同一 k 下两者谁更稳。常见现象是 LDA 在 k2 时明显优于 PCA 的 k2因为 LDA 用了标签但 PCA 的 k3 可能追上来因为保留了更多方差。这些差异要结合混淆矩阵看鸢尾花第二、三类之间的误判是主要误差来源第一类几乎总能分对。把这张比对表做出来作业的分析部分就立住了。5. 避坑与排查这份作业里最容易翻车的五个地方5.1 协方差矩阵奇异导致求逆报错现象跑马氏距离或 LDA 时抛 LinAlgError: Singular matrix。原因某个类样本数少于特征维度或特征之间完全共线协方差矩阵不满秩。解决加对角正则项正则系数从 1e-6 起试数据标准化后可适当减小或者改用 np.linalg.pinv 求伪逆但伪逆会改变距离的几何含义报告里要说明。5.2 标准化与否导致 PCA 结果对不上现象自己算的主成分方差贡献率和报告差很多。原因报告用相关矩阵等价于先标准化而代码里直接对原始数据求协方差量纲大的列主导了主成分。解决统一先做 (X-mean)/std 再求协方差ddof 取 1如果坚持用协方差矩阵报告里要注明且不能和用相关矩阵的结果混着比较。5.3 刀切法里模型复用造成信息泄漏现象留一法正确率高得离谱接近 100%。原因在循环外先在全量数据上 fit 了 PCA 或 GMM循环内只做 transform 和 predict测试样本参与了降维方向的估计。解决每次循环内重新 fit 整个流水线包括标准化、降维、聚类确保被剔除的样本对模型零贡献。5.4 GMM 不收敛或聚类标签与真实类别错位现象converged_ 为 False或 predict 出来的标签 0/1/2 和真实类别对不上。原因初始化不好、迭代次数不够或者聚类本身是无监督的成分编号和类别编号没有天然对应。解决用 kmeans 初始化、调大 max_iter标签错位不是 bug评估时要用匈牙利算法或混淆矩阵做最优匹配后再算正确率不能直接逐元素比。5.5 LDA 降维维度超过类别数减一现象LDA 投影后维度设成 3 或 4结果报错或得到无意义的方向。原因类间散度矩阵的秩最多是类别数减一三分类最多两个有效判别方向。解决n_components 不超过 min(特征维度, 类别数-1)超出部分特征值接近零对应的特征向量是数值噪声不要拿来用。6. 把评估做扎实一个可复用的比对脚本与我的固定习惯真正让这份作业从“能跑”变成“可信”的是把比对做成可复现的脚本而不是手工改参数跑几次记数字。我一般会写一个统一入口把降维方法、降维维度、聚类成分数作为参数内部固定随机种子输出正确率、混淆矩阵和运行时间三列这样换任何一组配置都能直接对比。下面这个脚本把 PCA 和 LDA 在 k1、2、3 下的留一法结果一次性打出来PCA 的 k 上限受特征维度约束LDA 的 k 上限受类别数约束脚本里做了截断保护。import time def build_predictor(method, k): def predict(X_train, y_train, x_test): if method pca: reducer PCA(n_componentsk).fit(X_train) else: W lda_fit(X_train, y_train, n_componentsk) reducer type(LDA, (), { transform: lambda self, X: X W })() Z_train reducer.transform(X_train) Z_test reducer.transform(x_test.reshape(1, -1)) gmm GaussianMixture(n_components3, covariance_typefull, random_state42).fit(Z_train) return gmm.predict(Z_test)[0] return predict results [] for method in [pca, lda]: max_k 4 if method pca else 2 for k in range(1, max_k 1): t0 time.time() acc loo_accuracy(X, y, build_predictor(method, k)) results.append((method, k, round(acc, 4), round(time.time() - t0, 2))) for r in results: print(f{r[0]:4} k{r[1]} acc{r[2]} time{r[3]}s)这个脚本的价值在于把“降维方法 × 降维维度”这个二维网格一次性铺开正确率和耗时并排看选型时不用来回改代码。参数上唯一要留意的是 LDA 的 max_k 写死为 2这是三分类的硬上限换成别的数据集要按类别数减一重算。跑完之后把结果整理成表格贴进报告再配一张 PCA 前两维的散点图按真实类别着色分析部分就有据可依了。血泪经验是这类统计作业的坑几乎都不在公式推导而在数据流的一致性——标准化用没用、ddof 取几、循环里模型有没有复用、标签映射对不对任何一处不一致都会让最终数字和预期差一截而且很难从结果反推是哪一步出的问题。从那以后我每次做这类留一法评估都强制先把全流程在一个样本上手工走一遍确认中间量的形状和数值合理再放开循环跑全量。希望帮到你。本文还有配套的精品资源点击获取