简介：这份数理统计大作业资源面向高校学生与数据分析初学者，围绕经典鸢尾花数据集展开完整分析，帮助读者理解多重变量分析在实际问题中的落地方式。内容涵盖马氏距离、混合高斯模型、主成分分析、线性判别分析及刀切法等核心知识点，并给出数据预处理、降维去噪、聚类分类与模型评估的完整流程，可作为课程作业参考或自学范例。资源包共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', header=None, names=col_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 对应 0，Versicolour 对应 1，Virginica 对应 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, axis=0) covs[c] = np.cov(X_c, rowvar=False) # 对第一个样本做判别 sample = X[0] distances = {c: mahalanobis_distance(sample, means[c], covs[c]) for c in classes} predicted = min(distances, key=distances.get) print(f"样本真实类别: {y[0]}, 马氏距离判别结果: {predicted}") print(f"到各类别的马氏距离: {distances}")这段代码的关键在于np.cov默认按行作为变量，所以rowvar=False表示每列是一个特征。实际跑的时候会发现，用马氏距离直接做最近邻判别，在鸢尾花上的准确率大概在 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, axis=0) X_centered = X - self.mean cov = np.cov(X_centered, rowvar=False) 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_components=2) 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, axis=0) 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, axis=0) 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_components=2) X_lda = lda.fit(X, y).transform(X) print(f"LDA 投影后前两维的类别均值:") for c in np.unique(y): print(f" 类别 {c}: {np.mean(X_lda[y == c], axis=0)}")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_iter=100, tol=1e-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, axis=1)) return norm * exp_term def fit(self, X): n, d = X.shape # 初始化：随机选 n_components 个样本作为均值 idx = np.random.choice(n, self.n_components, replace=False) 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(axis=1, keepdims=True) # M 步：更新参数 N_k = responsibilities.sum(axis=0) 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) ], axis=0))) if np.abs(log_likelihood - log_likelihood_old) < self.tol: print(f"EM 在第 {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, axis=1)E 步计算每个样本属于每个高斯成分的后验概率，M 步用这些概率加权更新均值、协方差和权重。_gaussian_pdf里用np.sum(diff @ inv_cov * diff, axis=1)计算马氏距离的平方，这是多元高斯密度函数的核心。收敛条件用对数似然的变化量判断，tol=1e-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), dtype=int) 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_components=3) 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_components=3) gmm_lda.fit(X_lda) y_pred_lda = gmm_lda.predict(X_lda) acc_lda = cluster_accuracy(y, y_pred_lda) print(f"PCA + GMM 聚类准确率: {acc_pca:.4f}") print(f"LDA + 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, classifier='mahalanobis'): """刀切法评估：每次留一个样本做测试""" n = X.shape[0] errors = 0 for i in range(n): X_train = np.delete(X, i, axis=0) y_train = np.delete(y, i) X_test = X[i:i+1] y_test = y[i] # 降维 reducer = reducer_class(n_components=n_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], axis=0) for c in classes} covs = {c: np.cov(X_train_reduced[y_train == c], rowvar=False) for c in classes} distances = {c: mahalanobis_distance(X_test_reduced[0], means[c], covs[c]) for c in classes} pred = min(distances, key=distances.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_components=3，np.linalg.eig会返回复数特征值，投影结果不可用。解决办法是限制n_components <= n_classes - 1。
现象四：刀切法跑得太慢。150 个样本跑 150 次训练，每次都要重新计算协方差矩阵和特征分解，在 Python 里可能要几十秒。解决办法是把降维步骤提前算好，或者用joblib并行化。
现象五：标准化后再做 PCA 和直接做 PCA 结果差异很大。原因是鸢尾花四个特征的量纲虽然都是厘米，但花瓣长度的方差明显大于花萼宽度，不标准化的话 PCA 会被大方差特征主导。解决办法是统一用StandardScaler处理后再降维。
6. 从判错率到模型选择：一个容易被忽略的验证技巧
刀切法跑完得到判错率之后，很多人直接取最小值对应的模型就结束了。但这里有个细节：判错率本身是个估计量，也有方差。150 个样本的留一法，判错率的标准误大约是 ( \sqrt{p(1-p)/n} )，当 ( p=0.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, exact=True) return result.pvalue # 假设已经得到 PCA 和 LDA 的留一法预测结果 # p_value = mcnemar_test(y, y_pred_pca_loo, y_pred_lda_loo) # print(f"McNemar 检验 p 值: {p_value:.4f}")如果 p 值大于 0.05，说明两个模型的判错率差异不显著，选哪个都行，优先选计算量小的。如果 p 值小于 0.05，才能说 LDA 显著优于 PCA。这个检验在作业里没要求，但实际做模型对比时很有用，能避免被随机波动带偏。
从那以后我每次做留一法对比，都会顺手跑一遍 McNemar 检验，确认差异不是玄学。希望帮到你。
本文还有配套的精品资源，点击获取
