简介这份数理统计大作业资源面向高校学生与数据分析初学者围绕经典鸢尾花数据集展开完整分析帮助读者理解多重变量分析在实际问题中的落地方式。内容涵盖马氏距离、混合高斯模型、主成分分析、线性判别分析及刀切法等核心知识点并给出数据预处理、降维去噪、聚类分类与模型评估的完整流程可作为课程作业参考或自学范例。资源包共1个docx文档约511KB内含摘要、算法原理、数据处理与结果比对等章节结构完整、条理清晰。目前已有1424人学习下载适合需要完成数理统计课程大作业、或希望系统梳理降维与聚类方法的读者参考借鉴。1. 从鸢尾花到刀切法一份数理统计大作业的完整复现路径鸢尾花数据集大概是每个学统计或机器学习的人绕不开的第一课但多数人只停留在调个sklearn的load_iris()然后跑个分类器看准确率。这份孙海燕老师布置的数理统计大作业要求的东西比调包深得多用马氏距离做判别、用高斯混合模型做聚类、用主成分分析和线性判别分析做降维最后用刀切法留一法去量化不同降维方式对分类效果的影响。整套流程覆盖了从数据清洗、降维、概率建模到模型评估的完整链路适合正在做统计课程设计、需要一份可复现参考实现的人。如果你手头有鸢尾花数据但不知道怎么把 PCA、LDA、GMM 串成一条分析线这份作业的代码和思路可以直接拿来拆解。2. 数据准备与马氏距离判别为什么不能用欧氏距离直接算2.1 鸢尾花数据的加载与标准化处理鸢尾花数据集共 150 个样本3 个类别各 50 个每个样本有 4 个特征花萼长度、花萼宽度、花瓣长度、花瓣宽度。原始数据从 UCI 仓库下载后是.data格式没有表头需要手动指定列名。常见做法是用pandas读入后直接做标准化因为后面 PCA 和 LDA 对量纲敏感。import pandas as pd import numpy as np from sklearn.preprocessing import StandardScaler # 列名按 UCI 原始顺序花萼长度、花萼宽度、花瓣长度、花瓣宽度、类别 col_names [sepal_length, sepal_width, petal_length, petal_width, label] df pd.read_csv(iris.data, headerNone, namescol_names) # 特征与标签分离 X_raw df.iloc[:, :4].values y df[label].astype(category).cat.codes.values # 0/1/2 编码 # 标准化均值为 0方差为 1 scaler StandardScaler() X scaler.fit_transform(X_raw) print(f样本数: {X.shape[0]}, 特征数: {X.shape[1]}) print(f各类别样本数: {np.bincount(y)})这里用StandardScaler而不是手动减均值除标准差是因为后面计算协方差矩阵和相关矩阵时标准化后的数据能直接消除量纲影响。注意cat.codes会把类别按字母顺序编码Setosa 对应 0Versicolour 对应 1Virginica 对应 2和作业里的符号说明一致。2.2 马氏距离的计算与判别逻辑马氏距离和欧氏距离的核心区别在于它考虑了特征之间的相关性。欧氏距离假设各维度独立且等方差但鸢尾花的四个特征明显相关——花瓣长的花花萼往往也长。马氏距离的公式是 ( D_M(x, \mu) \sqrt{(x-\mu)^T \Sigma^{-1} (x-\mu)} )其中 (\Sigma) 是协方差矩阵。def mahalanobis_distance(x, mu, cov): 计算样本 x 到均值 mu 的马氏距离 diff x - mu inv_cov np.linalg.inv(cov) return np.sqrt(diff.T inv_cov diff) # 按类别计算均值和协方差 classes np.unique(y) means {} covs {} for c in classes: X_c X[y c] means[c] np.mean(X_c, axis0) covs[c] np.cov(X_c, rowvarFalse) # 对第一个样本做判别 sample X[0] distances {c: mahalanobis_distance(sample, means[c], covs[c]) for c in classes} predicted min(distances, keydistances.get) print(f样本真实类别: {y[0]}, 马氏距离判别结果: {predicted}) print(f到各类别的马氏距离: {distances})这段代码的关键在于np.cov默认按行作为变量所以rowvarFalse表示每列是一个特征。实际跑的时候会发现用马氏距离直接做最近邻判别在鸢尾花上的准确率大概在 95% 左右错的主要是 Versicolour 和 Virginica 之间的边界样本。这也解释了为什么后面要引入 GMM——马氏距离只用了均值和协方差而 GMM 能建模更复杂的分布形状。注意计算协方差矩阵前一定要确保每个类别的样本数大于特征数否则协方差矩阵奇异求逆会报错。鸢尾花每类 50 个样本、4 个特征不存在这个问题但换成高维小样本数据就得先降维或用伪逆。3. PCA 与 LDA 降维实战从 4 维到 2 维的信息保留率3.1 PCA 的实现与主成分个数选择PCA 的目标是找到一组正交基使得数据投影后的方差最大。具体做法是对协方差矩阵做特征分解取最大的几个特征值对应的特征向量。作业里要求自己实现 PCA 类而不是直接调sklearn.decomposition.PCA目的是理解特征值排序和投影的每一步。class PCA: def __init__(self, n_components): self.n_components n_components self.components None self.mean None self.explained_variance_ratio None def fit(self, X): self.mean np.mean(X, axis0) X_centered X - self.mean cov np.cov(X_centered, rowvarFalse) eigenvalues, eigenvectors np.linalg.eigh(cov) # eigh 返回升序需要反转 idx np.argsort(eigenvalues)[::-1] eigenvalues eigenvalues[idx] eigenvectors eigenvectors[:, idx] self.components eigenvectors[:, :self.n_components] self.explained_variance_ratio eigenvalues[:self.n_components] / np.sum(eigenvalues) return self def transform(self, X): X_centered X - self.mean return X_centered self.components pca PCA(n_components2) X_pca pca.fit(X).transform(X) print(f前两个主成分解释方差比例: {pca.explained_variance_ratio}) print(f累计解释方差: {np.sum(pca.explained_variance_ratio):.4f})跑完会发现前两个主成分累计解释了约 95% 的方差这意味着从 4 维降到 2 维只损失了 5% 的信息。np.linalg.eigh用于对称矩阵比通用的eig更稳定且返回实特征值。特征值排序后取前n_components个对应的特征向量就是投影方向。3.2 LDA 的类间散度与类内散度矩阵LDA 和 PCA 的本质区别在于 LDA 用了标签信息。它要最大化类间散度与类内散度的比值也就是让不同类别的投影点尽可能分开同一类别的投影点尽可能聚集。对于多分类问题LDA 最多能降到类别数 - 1维鸢尾花有 3 类所以最多降到 2 维。class LDA: def __init__(self, n_components): self.n_components n_components self.components None def fit(self, X, y): n_features X.shape[1] classes np.unique(y) mean_overall np.mean(X, axis0) S_W np.zeros((n_features, n_features)) S_B np.zeros((n_features, n_features)) for c in classes: X_c X[y c] mean_c np.mean(X_c, axis0) S_W (X_c - mean_c).T (X_c - mean_c) n_c X_c.shape[0] mean_diff (mean_c - mean_overall).reshape(-1, 1) S_B n_c * (mean_diff mean_diff.T) # 求解 S_W^{-1} S_B 的特征值 eig_vals, eig_vecs np.linalg.eig(np.linalg.inv(S_W) S_B) eig_vals np.real(eig_vals) eig_vecs np.real(eig_vecs) idx np.argsort(eig_vals)[::-1] eig_vecs eig_vecs[:, idx] self.components eig_vecs[:, :self.n_components] return self def transform(self, X): return X self.components lda LDA(n_components2) X_lda lda.fit(X, y).transform(X) print(fLDA 投影后前两维的类别均值:) for c in np.unique(y): print(f 类别 {c}: {np.mean(X_lda[y c], axis0)})LDA 的核心是求解 ( S_W^{-1} S_B ) 的特征向量。np.linalg.eig返回的特征值可能是复数因为数值计算误差所以用np.real取实部。投影后可以看到三个类别的均值在 LDA 空间里分得很开尤其是 Setosa 和其他两类几乎完全分离。这也是为什么 LDA 在鸢尾花上做分类通常比 PCA 效果好——它利用了标签信息来指导降维方向。提示如果S_W不可逆常见做法是加一个小的正则项比如S_W 1e-6 * np.eye(n_features)或者先做 PCA 降维再用 LDA。4. 高斯混合模型聚类从 EM 算法到马氏距离判别4.1 GMM 的 EM 算法实现高斯混合模型假设数据由若干个高斯分布混合生成每个高斯成分有自己的均值、协方差和权重。EM 算法通过交替执行 E 步计算后验概率和 M 步更新参数来最大化似然函数。作业里要求用极大似然估计求参数下面是一个简化版的 GMM 实现。class GMM: def __init__(self, n_components, max_iter100, tol1e-4): self.n_components n_components self.max_iter max_iter self.tol tol def _gaussian_pdf(self, X, mean, cov): d X.shape[1] diff X - mean inv_cov np.linalg.inv(cov) det_cov np.linalg.det(cov) norm 1.0 / np.sqrt((2 * np.pi) ** d * det_cov) exp_term np.exp(-0.5 * np.sum(diff inv_cov * diff, axis1)) return norm * exp_term def fit(self, X): n, d X.shape # 初始化随机选 n_components 个样本作为均值 idx np.random.choice(n, self.n_components, replaceFalse) self.means X[idx] self.covs [np.eye(d) for _ in range(self.n_components)] self.weights np.ones(self.n_components) / self.n_components log_likelihood_old 0 for iteration in range(self.max_iter): # E 步计算后验概率 responsibilities np.zeros((n, self.n_components)) for k in range(self.n_components): responsibilities[:, k] self.weights[k] * self._gaussian_pdf(X, self.means[k], self.covs[k]) responsibilities / responsibilities.sum(axis1, keepdimsTrue) # M 步更新参数 N_k responsibilities.sum(axis0) for k in range(self.n_components): self.means[k] (responsibilities[:, k] X) / N_k[k] diff X - self.means[k] self.covs[k] (responsibilities[:, k] * diff.T) diff / N_k[k] self.weights[k] N_k[k] / n # 计算对数似然 log_likelihood np.sum(np.log(np.sum([ self.weights[k] * self._gaussian_pdf(X, self.means[k], self.covs[k]) for k in range(self.n_components) ], axis0))) if np.abs(log_likelihood - log_likelihood_old) self.tol: print(fEM 在第 {iteration} 次迭代收敛) break log_likelihood_old log_likelihood return self def predict(self, X): responsibilities np.zeros((X.shape[0], self.n_components)) for k in range(self.n_components): responsibilities[:, k] self.weights[k] * self._gaussian_pdf(X, self.means[k], self.covs[k]) return np.argmax(responsibilities, axis1)E 步计算每个样本属于每个高斯成分的后验概率M 步用这些概率加权更新均值、协方差和权重。_gaussian_pdf里用np.sum(diff inv_cov * diff, axis1)计算马氏距离的平方这是多元高斯密度函数的核心。收敛条件用对数似然的变化量判断tol1e-4是常见取值。4.2 降维后 GMM 聚类效果对比把 PCA 和 LDA 降维后的数据分别喂给 GMM观察聚类结果和真实标签的匹配程度。由于 GMM 是无监督的聚类编号和真实类别编号不一定对应需要用匈牙利算法做标签匹配。from scipy.optimize import linear_sum_assignment def cluster_accuracy(y_true, y_pred): 用匈牙利算法匹配聚类标签和真实标签 n_classes max(y_true.max(), y_pred.max()) 1 confusion np.zeros((n_classes, n_classes), dtypeint) for t, p in zip(y_true, y_pred): confusion[t, p] 1 row_ind, col_ind linear_sum_assignment(-confusion) return confusion[row_ind, col_ind].sum() / len(y_true) # PCA 降维后 GMM 聚类 gmm_pca GMM(n_components3) gmm_pca.fit(X_pca) y_pred_pca gmm_pca.predict(X_pca) acc_pca cluster_accuracy(y, y_pred_pca) # LDA 降维后 GMM 聚类 gmm_lda GMM(n_components3) gmm_lda.fit(X_lda) y_pred_lda gmm_lda.predict(X_lda) acc_lda cluster_accuracy(y, y_pred_lda) print(fPCA GMM 聚类准确率: {acc_pca:.4f}) print(fLDA GMM 聚类准确率: {acc_lda:.4f})实际跑下来LDA GMM 的准确率通常比 PCA GMM 高几个百分点因为 LDA 降维时已经利用了标签信息投影后的数据类别可分性更强。但 GMM 本身是无监督的所以即使 LDA 用了标签GMM 聚类时并没有用到标签这个对比仍然有意义。linear_sum_assignment来自 scipy用于求解二分图最大匹配这里把混淆矩阵取负后求最小代价匹配等价于最大化正确匹配数。注意GMM 对初始化敏感不同的随机种子可能得到不同的聚类结果。建议多跑几次取平均或者用 k-means 的聚类中心作为 GMM 的初始均值。5. 刀切法评估与避坑指南留一法到底怎么用才不翻车5.1 刀切法的实现与降维方法对比刀切法留一法的核心思想是每次留一个样本作为测试集其余样本训练模型重复 n 次后统计判错率。虽然计算量大但对小数据集来说能给出更稳定的评估结果。下面用刀切法对比 PCA 和 LDA 在不同降维维度下的分类效果。def jackknife_evaluation(X, y, reducer_class, n_components, classifiermahalanobis): 刀切法评估每次留一个样本做测试 n X.shape[0] errors 0 for i in range(n): X_train np.delete(X, i, axis0) y_train np.delete(y, i) X_test X[i:i1] y_test y[i] # 降维 reducer reducer_class(n_componentsn_components) if reducer_class LDA: reducer.fit(X_train, y_train) else: reducer.fit(X_train) X_train_reduced reducer.transform(X_train) X_test_reduced reducer.transform(X_test) # 用马氏距离做判别 classes np.unique(y_train) means {c: np.mean(X_train_reduced[y_train c], axis0) for c in classes} covs {c: np.cov(X_train_reduced[y_train c], rowvarFalse) for c in classes} distances {c: mahalanobis_distance(X_test_reduced[0], means[c], covs[c]) for c in classes} pred min(distances, keydistances.get) if pred ! y_test: errors 1 return errors / n # 对比不同降维维度 for n_comp in [1, 2]: err_pca jackknife_evaluation(X, y, PCA, n_comp) err_lda jackknife_evaluation(X, y, LDA, n_comp) print(f降维到 {n_comp} 维: PCA 判错率{err_pca:.4f}, LDA 判错率{err_lda:.4f})这段代码里np.delete每次剔除一个样本然后重新训练降维器和判别器。注意 LDA 的fit需要标签所以单独判断了reducer_class LDA。跑完会发现 LDA 降到 2 维时判错率最低PCA 降到 2 维时判错率略高但差距不大。如果降到 1 维两者判错率都会上升因为信息损失太多。5.2 常见踩坑记录现象一协方差矩阵奇异导致np.linalg.inv报错。原因是个别类别的样本数少于特征数或者特征之间存在完全共线性。解决办法是先用 PCA 降维或者在协方差矩阵上加一个小的对角正则项cov 1e-6 * np.eye(d)。现象二GMM 的 EM 算法不收敛或收敛到局部最优。原因是初始化均值太接近或者协方差矩阵初始值设得太大。解决办法是用 k-means 先跑一遍得到初始均值协方差初始值设为各类别的样本协方差。现象三LDA 降维后维度超过类别数 - 1。比如鸢尾花 3 类LDA 最多降到 2 维如果设n_components3np.linalg.eig会返回复数特征值投影结果不可用。解决办法是限制n_components n_classes - 1。现象四刀切法跑得太慢。150 个样本跑 150 次训练每次都要重新计算协方差矩阵和特征分解在 Python 里可能要几十秒。解决办法是把降维步骤提前算好或者用joblib并行化。现象五标准化后再做 PCA 和直接做 PCA 结果差异很大。原因是鸢尾花四个特征的量纲虽然都是厘米但花瓣长度的方差明显大于花萼宽度不标准化的话 PCA 会被大方差特征主导。解决办法是统一用StandardScaler处理后再降维。6. 从判错率到模型选择一个容易被忽略的验证技巧刀切法跑完得到判错率之后很多人直接取最小值对应的模型就结束了。但这里有个细节判错率本身是个估计量也有方差。150 个样本的留一法判错率的标准误大约是 ( \sqrt{p(1-p)/n} )当 ( p0.05 ) 时标准误约 0.018。这意味着两个模型的判错率差 1 个百分点可能只是随机波动不能说明谁真的更好。我一般会补一个 McNemar 检验专门比较两个模型在同一批样本上的判错情况是否显著不同。具体做法是统计两个模型都判错、一个判错一个判对、都判对的样本数构造 2x2 列联表然后用卡方检验。from statsmodels.stats.contingency_tables import mcnemar def mcnemar_test(y_true, pred_a, pred_b): McNemar 检验比较两个模型的判错差异 both_wrong np.sum((pred_a ! y_true) (pred_b ! y_true)) a_wrong_b_right np.sum((pred_a ! y_true) (pred_b y_true)) a_right_b_wrong np.sum((pred_a y_true) (pred_b ! y_true)) both_right np.sum((pred_a y_true) (pred_b y_true)) table [[both_right, a_right_b_wrong], [a_wrong_b_right, both_wrong]] result mcnemar(table, exactTrue) return result.pvalue # 假设已经得到 PCA 和 LDA 的留一法预测结果 # p_value mcnemar_test(y, y_pred_pca_loo, y_pred_lda_loo) # print(fMcNemar 检验 p 值: {p_value:.4f})如果 p 值大于 0.05说明两个模型的判错率差异不显著选哪个都行优先选计算量小的。如果 p 值小于 0.05才能说 LDA 显著优于 PCA。这个检验在作业里没要求但实际做模型对比时很有用能避免被随机波动带偏。从那以后我每次做留一法对比都会顺手跑一遍 McNemar 检验确认差异不是玄学。希望帮到你。本文还有配套的精品资源点击获取