简介本资源是高校数理统计课程的大作业实践报告面向统计学、数据科学及相关专业本科生聚焦经典鸢尾花数据集的多方法综合分析。报告系统应用马氏距离度量样本相似性、混合高斯模型实现无监督聚类、主成分分析与线性判别分析进行降维与可分性提升并采用刀切法留一法评估不同降维策略与模型的泛化性能完整覆盖课程核心知识点与建模流程。资源为1个511KB的Word文档.docx内容结构严谨含摘要、引言、六大算法原理详述3.1–3.5节、四阶段数据处理分析4.1–4.4节及关键词与目录公式推导、代码逻辑说明与结果比对清晰可循。目前已有1426人学习下载适合课程复习、大作业参考、统计建模入门与算法对比研究。1. 这份数理统计大作业不是跑通 sklearn 就算交差它要你亲手推导马氏距离、手写 GMM EM 步骤、用刀切法硬刚每一种降维组合的泛化能力很多同学拿到“孙海燕老师数理统计大作业”这个标题第一反应是不就是拿 iris 数据集调个PCA().fit_transform()GaussianMixture().fit()吗三小时出图五小时写完摘要答辩前一晚改格式——结果被导师一句“你这个马氏距离是直接调scipy.spatial.distance.mahalanobis吗协方差矩阵求逆前检查过条件数没”当场卡住。这不是编程课设这是数理统计的原理验证型大作业它要求你把 PCA 的特征向量投影、LDA 的类间/类内散度矩阵、GMM 的 E-M 迭代、刀切法Jackknife的留一循环全部拆解到矩阵运算层面再用 Python 一行行实现。它面向的是某高校数理统计课程的真实教学目标——不考你会不会用 API而考你是否真正理解每个公式的物理意义、数值稳定性边界和统计假设前提。适合两类人一是正在啃《多元统计分析》教材、卡在 LDA 推导第 3 行的同学二是想用一份可答辩、可复现、带完整推导链的方案替代网上零散博客拼凑的毕设/课设学生。它解决的不是“怎么画出分类图”而是“当你的 PCA 降维后 GMM 聚类效果反而变差时你能否定位是协方差矩阵病态、还是 LDA 投影方向未对齐类别可分性、抑或刀切法暴露了模型过拟合”。2. 总体架构设计四层结构与五个核心模块的职责划分这份大作业不是单点算法堆砌而是一个闭环验证系统。我把它拆成四层结构数据层 → 降维层 → 建模层 → 验证层每一层都承担明确的统计任务且层间接口必须显式暴露数学含义。比如降维层输出的不能是黑盒数组而必须是带维度标签的pd.DataFrame其列名需体现“PC1解释方差比 72.9%”或“LDA1最大化类间离差比”建模层输入必须接收该结构并在拟合时校验输入维度是否匹配 GMM 成分数验证层则强制所有模块支持predict_proba和score接口为刀切法统一调用铺路。整个流程围绕五个核心模块展开它们不是孤立函数而是有严格依赖关系的组件数据预处理模块负责加载 UCI Iris 原始 CSV执行中心化非标准化因马氏距离对尺度敏感但需保证各变量同量纲、标签编码Setosa→0, Versicolour→1, Virginica→2并生成训练/验证分离的索引映射PCA 降维模块不调用sklearn.decomposition.PCA而是手动计算相关矩阵R (X^T X)/(n-1)调用np.linalg.eig(R)获取特征值/向量按累计解释方差比 ≥85% 自动截断主成分数通常为 2并返回正交投影矩阵W_pca ∈ ℝ^{4×2}LDA 降维模块手动构建类内散度矩阵S_w Σ_i Σ_{x∈C_i} (x - μ_i)(x - μ_i)^T与类间散度矩阵S_b Σ_i n_i (μ_i - μ)(μ_i - μ)^T求解广义特征值问题S_w^{-1} S_b v λv取前k个最大特征值对应特征向量组成W_lda ∈ ℝ^{4×2}GMM 建模模块基于降维后数据实现完整的 EM 算法E 步计算后验概率γ(z_{nk}) π_k N(x_n|μ_k,Σ_k) / Σ_j π_j N(x_n|μ_j,Σ_j)M 步更新π_k Σ_n γ(z_{nk})/N,μ_k Σ_n γ(z_{nk}) x_n / Σ_n γ(z_{nk}),Σ_k Σ_n γ(z_{nk}) (x_n - μ_k)(x_n - μ_k)^T / Σ_n γ(z_{nk})收敛判据为对数似然增量 1e-6刀切法验证模块对每个样本i临时剔除它用剩余n-1个样本重新训练 PCA/LDAGMM 全流程再用该模型预测x_i所属类别记录错误数最终计算总错判率err_rate k_total / 150。提示所有模块必须返回dict结构包含model模型参数字典、transform_matrix若适用、log_likelihoodGMM、error_count刀切法。这能避免答辩时被问“你 PCA 选了几个主成分依据是什么”时只答“我试了2个效果好”——你得拿出explained_variance_ratio_数组和累计曲线。2.1 数据预处理为什么中心化是必须的而标准化反而是陷阱Iris 数据四个属性单位不同厘米但马氏距离的核心前提是变量已处于同一量纲下且协方差矩阵可逆。直接标准化z-score会破坏原始数据的协方差结构——因为cov(z_i,z_j) cov((x_i-μ_i)/σ_i, (x_j-μ_j)/σ_j) cov(x_i,x_j)/(σ_i σ_j)这导致后续 PCA 特征向量不再反映原始变量间的相关性强度。正确做法是仅做中心化zero-mean保留原始尺度import numpy as np import pandas as pd # 加载原始 UCI Iris 数据150×5最后一列为字符串标签 df pd.read_csv(iris.data, headerNone, names[sepal_length,sepal_width,petal_length,petal_width,class]) X_raw df.iloc[:, :4].values # shape: (150, 4) y_raw df[class].map({Iris-setosa:0, Iris-versicolor:1, Iris-virginica:2}).values # 仅中心化减去均值不除标准差 X_centered X_raw - np.mean(X_raw, axis0) # shape: (150, 4) # 验证各列均值应接近0 print(Centered data column means:, np.mean(X_centered, axis0)) # [~0, ~0, ~0, ~0]这段代码的关键在于np.mean(X_raw, axis0)计算的是每列即每个属性的均值然后从整列中减去该值。它不改变变量间的协方差比例却为后续计算马氏距离d_M(x,y) sqrt((x-y)^T Σ^{-1} (x-y))提供了必要前提Σ必须是X_centered的协方差矩阵。如果误用标准化Σ会变成相关矩阵R此时马氏距离退化为欧氏距离因RI时d_M d_E完全失去度量多变量相关性的意义——这正是答辩时最容易被揪住的“概念性翻车点”。2.2 PCA 降维模块从相关矩阵到投影矩阵的完整推导链PCA 的本质是找一个正交基{u_1,...,u_k}使得数据在该基上的投影方差最大。教科书常写“求协方差矩阵C (X^T X)/(n-1)的特征向量”但对 Iris 这种小样本n150, p4直接算X^T X可能引入数值误差。更稳健的做法是先中心化X_c再计算相关矩阵R (X_c^T X_c)/(n-1)因各变量量纲一致R与C成比例特征向量相同。以下是手动实现 PCA 的核心步骤def manual_pca(X_centered, n_components2): n, p X_centered.shape # 1. 计算相关矩阵 R (p x p) R np.dot(X_centered.T, X_centered) / (n - 1) # 无偏估计 # 2. 求特征值分解R V diag(λ) V.T eigenvals, eigenvecs np.linalg.eig(R) # 3. 按特征值降序排列最大方差优先 idx np.argsort(eigenvals)[::-1] eigenvals eigenvals[idx] eigenvecs eigenvecs[:, idx] # 列向量为特征向量 # 4. 截取前 n_components 个主成分投影方向 W_pca eigenvecs[:, :n_components] # shape: (p, n_components) # 5. 投影数据X_pca X_centered W_pca X_pca np.dot(X_centered, W_pca) # 6. 计算解释方差比 explained_ratio eigenvals[:n_components] / np.sum(eigenvals) return { W_pca: W_pca, X_pca: X_pca, eigenvals: eigenvals, explained_ratio: explained_ratio, cumsum_ratio: np.cumsum(explained_ratio) } # 调用示例 pca_result manual_pca(X_centered, n_components2) print(Explained variance ratio for PC1 PC2:, pca_result[explained_ratio]) print(Cumulative sum:, pca_result[cumsum_ratio][-1]) # 通常 0.95这段代码的深层逻辑是R的第i个特征值λ_i表示数据在第i个主成分方向u_i上的投影方差u_i是R的第i个特征向量。因此W_pca的列就是我们要的正交基。注意np.linalg.eig返回的特征向量是列向量所以eigenvecs[:, idx]是正确的索引方式。若误写成eigenvecs[idx, :]得到的将是错误的行向量矩阵导致投影X_centered W_pca维度不匹配150×4 4×2 150×2这是初学者高频报错点。此外explained_ratio必须用eigenvals[:n_components] / np.sum(eigenvals)计算而非eigenvals[:n_components] / n_components——后者是常见误解混淆了“单个主成分方差”与“总方差”的关系。2.3 LDA 降维模块为什么 LDA 必须用带标签的类内/类间散度矩阵PCA 是无监督的它只看数据整体分布而 LDA 是有监督的它的目标是让同类样本投影后尽可能紧凑类内散度小异类样本投影后尽可能分离类间散度大。这决定了 LDA 的投影矩阵W_lda必须由标签y驱动。核心公式是最大化 Rayleigh 商J(W) tr(W^T S_b W) / tr(W^T S_w W)其中S_b是类间散度矩阵S_w是类内散度矩阵。手动实现的关键在于正确计算这两个矩阵def manual_lda(X_centered, y, n_components2): n, p X_centered.shape classes np.unique(y) n_classes len(classes) # 1. 计算全局均值 μ mu_global np.mean(X_centered, axis0) # shape: (p,) # 2. 初始化类内散度矩阵 S_w (p x p) 和类间散度矩阵 S_b (p x p) S_w np.zeros((p, p)) S_b np.zeros((p, p)) # 3. 对每个类别计算 S_w 和 S_b for cls in classes: X_cls X_centered[y cls] # 取出该类所有样本 n_cls X_cls.shape[0] mu_cls np.mean(X_cls, axis0) # 该类均值 # 类内散度S_w Σ (x - μ_cls)(x - μ_cls)^T diff_cls X_cls - mu_cls S_w np.dot(diff_cls.T, diff_cls) # 类间散度S_b n_cls * (μ_cls - μ_global)(μ_cls - μ_global)^T diff_global mu_cls - mu_global S_b n_cls * np.outer(diff_global, diff_global) # 4. 求解广义特征值问题 S_w^{-1} S_b v λv # 注意S_w 可能奇异需用 np.linalg.pinv 求伪逆 S_w_inv np.linalg.pinv(S_w) eigenvals, eigenvecs np.linalg.eig(np.dot(S_w_inv, S_b)) # 5. 按特征值降序排列越大表示类间分离越强 idx np.argsort(eigenvals)[::-1] eigenvals eigenvals[idx] eigenvecs eigenvecs[:, idx] # 6. 取前 n_components 个判别向量 W_lda eigenvecs[:, :n_components] # shape: (p, n_components) # 7. 投影数据 X_lda np.dot(X_centered, W_lda) return { W_lda: W_lda, X_lda: X_lda, eigenvals: eigenvals, S_w: S_w, S_b: S_b } # 调用示例 lda_result manual_lda(X_centered, y_raw, n_components2) print(LDA eigenvalues (largest first):, lda_result[eigenvals][:2])这里最易出错的是S_w和S_b的计算逻辑。S_w是所有类的类内离差平方和之和它衡量的是类内紧凑性S_b是各类均值相对于全局均值的加权离差平方和它衡量的是类间分离性。若漏掉n_cls权重如S_b np.outer(diff_global, diff_global)会导致S_b被低估LDA 方向偏向样本数少的类别。另一个坑是S_w可能奇异秩亏直接求逆np.linalg.inv(S_w)会报错必须用伪逆np.linalg.pinv。这也是为什么 LDA 在小样本、高维场景下不如 PCA 稳健——它对S_w的条件数极度敏感。3. 核心算法实现GMM 的 EM 迭代与马氏距离分类器的耦合设计GMM 不是黑箱它是对数据分布的显式建模。在 Iris 场景下我们假设每个类别c ∈ {0,1,2}对应一个高斯成分GMM 的目标是估计这三个成分的参数(π_c, μ_c, Σ_c)。但直接用sklearn.mixture.GaussianMixture会掩盖关键细节EM 算法的收敛性、协方差矩阵的约束全协方差对角球形、以及如何用马氏距离替代默认的predict_proba进行硬分类。本节将 GMM 与马氏距离深度耦合形成一条从概率建模到确定性判别的完整链路。3.1 GMM 的 EM 算法手写迭代拒绝黑盒EM 算法分为 EExpectation步和 MMaximization步。E 步计算每个样本属于每个成分的后验概率γ(z_{nk})M 步用这些后验概率加权更新参数。关键在于M 步更新的协方差矩阵Σ_k必须是正定的否则后续马氏距离计算会失败。因此我们在 M 步后添加协方差正则化加一个小的单位阵εIdef gmm_em(X, n_components3, max_iter100, tol1e-6, reg_covar1e-6): n, d X.shape # 初始化参数K-means 启发式 from sklearn.cluster import KMeans kmeans KMeans(n_clustersn_components, n_init10, random_state42) labels kmeans.fit_predict(X) pi_k np.array([np.sum(labelsk)/n for k in range(n_components)]) mu_k np.array([np.mean(X[labelsk], axis0) for k in range(n_components)]) # 初始协方差每个类的样本协方差 正则项 sigma_k np.array([ np.cov(X[labelsk].T) reg_covar * np.eye(d) for k in range(n_components) ]) log_likelihoods [] for iteration in range(max_iter): # E-step: 计算后验概率 γ(z_{nk}) gamma np.zeros((n, n_components)) for k in range(n_components): # 计算高斯概率密度 N(x_n | μ_k, Σ_k) diff X - mu_k[k] # (n, d) # Σ_k 是 dxd需求逆 try: sigma_inv np.linalg.inv(sigma_k[k]) except np.linalg.LinAlgError: # 若奇异加更大正则项 sigma_k[k] 1e-3 * np.eye(d) sigma_inv np.linalg.inv(sigma_k[k]) # 马氏距离平方: (x-μ)^T Σ^{-1} (x-μ) maha_sq np.sum(diff sigma_inv * diff, axis1) # (n,) # 高斯密度忽略归一化常数 norm_const np.sqrt(np.linalg.det(2*np.pi*sigma_k[k])) gauss_prob np.exp(-0.5 * maha_sq) / norm_const gamma[:, k] pi_k[k] * gauss_prob # 归一化 gamma使其每行和为1 gamma_sum np.sum(gamma, axis1, keepdimsTrue) gamma gamma / gamma_sum # M-step: 更新参数 N_k np.sum(gamma, axis0) # (n_components,) # 更新 π_k pi_k N_k / n # 更新 μ_k mu_k np.dot(gamma.T, X) / N_k[:, None] # (n_components, d) # 更新 Σ_k for k in range(n_components): diff X - mu_k[k] # (n, d) sigma_k[k] np.dot(gamma[:, k] * diff.T, diff) / N_k[k] # 添加正则化确保正定 sigma_k[k] reg_covar * np.eye(d) # 计算对数似然 log_likelihood np.sum(np.log(gamma_sum 1e-30)) log_likelihoods.append(log_likelihood) # 检查收敛 if iteration 0 and abs(log_likelihoods[-1] - log_likelihoods[-2]) tol: break return { pi_k: pi_k, mu_k: mu_k, sigma_k: sigma_k, gamma: gamma, log_likelihoods: log_likelihoods, n_iter: iteration 1 } # 示例在 PCA 降维后的数据上训练 GMM X_pca pca_result[X_pca] # (150, 2) gmm_result gmm_em(X_pca, n_components3, max_iter50) print(GMM converged in, gmm_result[n_iter], iterations) print(Final log-likelihood:, gmm_result[log_likelihoods][-1])这段代码的精妙之处在于它把 GMM 的核心计算E 步的后验概率、M 步的参数更新全部显式写出而非调用封装函数。特别是sigma_k的更新np.dot(gamma[:, k] * diff.T, diff) / N_k[k]是加权协方差的标准公式其中gamma[:, k]是每个样本属于第k类的置信度。reg_covar1e-6的正则项至关重要——它防止sigma_k奇异为后续马氏距离计算提供稳定基础。若跳过此步np.linalg.inv(sigma_k[k])在sigma_k[k]接近奇异时会崩溃这是 GMM 实现中最常见的“玄学翻车点”。3.2 马氏距离分类器用统计距离替代软概率直击判别本质GMM 默认的predict方法基于最大后验概率MAPargmax_k γ(z_{nk})。但这忽略了距离的几何意义。马氏距离d_M(x, μ_k) sqrt((x-μ_k)^T Σ_k^{-1} (x-μ_k))才是衡量样本x到第k类中心μ_k的“统计距离”——它考虑了该类的协方差结构Σ_k比欧氏距离更能反映真实分布形态。因此我们设计一个马氏距离分类器它不依赖γ(z_{nk})而是直接计算x到每个μ_k的马氏距离选择距离最小的类别def mahalanobis_classifier(X, mu_k, sigma_k): X: (n, d) 待分类样本 mu_k: (k, d) 各类均值 sigma_k: (k, d, d) 各类协方差矩阵 返回: (n,) 预测类别标签 n, d X.shape k mu_k.shape[0] distances np.zeros((n, k)) for i in range(k): try: sigma_inv np.linalg.inv(sigma_k[i]) except np.linalg.LinAlgError: # 协方差矩阵奇异时用伪逆 sigma_inv np.linalg.pinv(sigma_k[i]) diff X - mu_k[i] # (n, d) # 计算马氏距离平方(x-μ)^T Σ^{-1} (x-μ) # 使用 einsum 避免显式循环 maha_sq np.einsum(ij,jk,ik-i, diff, sigma_inv, diff) distances[:, i] np.sqrt(np.maximum(maha_sq, 0)) # 防止负数开方 # 返回距离最小的类别索引 return np.argmin(distances, axis1) # 在 PCA 降维数据上测试 y_pred_maha mahalanobis_classifier(X_pca, gmm_result[mu_k], gmm_result[sigma_k]) # 计算准确率 accuracy_maha np.mean(y_pred_maha y_raw) print(Mahalanobis classifier accuracy on PCA data:, accuracy_maha)这个分类器的设计哲学是GMM 建模的是概率但 Iris 分类的本质是判别。马氏距离直接量化了“样本离哪个类中心更近”且其计算过程复用了 GMM 已估计的μ_k和Σ_k无需额外训练。np.einsum(ij,jk,ik-i, diff, sigma_inv, diff)是高效计算马氏距离平方的向量化方法比双重 for 循环快一个数量级。np.maximum(maha_sq, 0)是安全防护——浮点误差可能导致maha_sq为极小负数开方会报nan。这种细节能让你的代码在答辩演示时稳如磐石。3.3 刀切法Jackknife验证逐样本剔除硬核检验泛化能力刀切法是比交叉验证更严格的模型评估方法尤其适合小样本n150。它的逻辑是每次只剔除一个样本用剩余 149 个样本重新训练整个 pipelinePCA/LDA → GMM → 马氏分类再用该模型预测被剔除的样本。150 次循环后统计总错判数。这种方法能暴露模型对单个样本的敏感性是检验“过拟合”的黄金标准。实现难点在于如何高效地“剔除一个样本”并重训我们采用索引掩码而非物理删除避免重复加载数据def jackknife_validation(X_centered, y,降维_func, gmm_func, classifier_func, n_components2): X_centered: (150, 4) 中心化数据 y: (150,) 标签 降维_func: 如 manual_pca 或 manual_lda接受 X 和 n_components gmm_func: 如 gmm_em接受降维后 X classifier_func: 如 mahalanobis_classifier接受降维后 X 和 GMM 参数 n len(y) errors 0 # 预分配存储避免循环中重复创建 X_pca_all np.zeros((n, n_components)) y_pred_all np.zeros(n, dtypeint) for i in range(n): # 创建掩码剔除第 i 个样本 mask np.ones(n, dtypebool) mask[i] False X_train X_centered[mask] # (149, 4) y_train y[mask] # (149,) # 1. 降维用训练集学习投影矩阵 dim_result 降维_func(X_train, n_componentsn_components) X_train_dim dim_result[X_pca] if X_pca in dim_result else dim_result[X_lda] # 2. GMM 训练 gmm_result gmm_func(X_train_dim, n_components3) # 3. 对被剔除的样本 x_i 进行预测 x_test X_centered[i:i1] # (1, 4) x_test_dim np.dot(x_test, dim_result[W_pca]) if W_pca in dim_result else np.dot(x_test, dim_result[W_lda]) y_pred classifier_func(x_test_dim, gmm_result[mu_k], gmm_result[sigma_k]) if y_pred[0] ! y[i]: errors 1 error_rate errors / n return { error_count: errors, error_rate: error_rate, accuracy: 1 - error_rate } # 执行 PCA GMM 马氏距离的刀切法 pca_jackknife jackknife_validation( X_centered, y_raw, manual_pca, gmm_em, mahalanobis_classifier, n_components2 ) print(PCAGMMMahalanobis Jackknife error rate:, pca_jackknife[error_rate]) # 执行 LDA GMM 马氏距离的刀切法 lda_jackknife jackknife_validation( X_centered, y_raw, manual_lda, gmm_em, mahalanobis_classifier, n_components2 ) print(LDAGMMMahalanobis Jackknife error rate:, lda_jackknife[error_rate])这段代码的工程价值在于它把刀切法封装成一个可插拔的验证器。只需传入不同的降维_funcmanual_pca或manual_lda就能对比两种降维策略的泛化能力。注意mask的使用——它通过布尔索引高效提取子集比np.delete(X_centered, i, axis0)快得多。x_test_dim的计算必须复用dim_result中的W_pca或W_lda确保测试样本与训练样本在同一个投影空间。若误用X_centered[i]直接降维会得到错误结果。刀切法耗时较长150×训练但这是值得的——它给出的error_rate是模型真实鲁棒性的铁证。4. 关键避坑指南5 个让答辩老师眼前一亮的“血泪经验”做这份大作业踩坑是常态。但有些坑踩一次就足以让答辩陷入被动。以下是我在指导某高校数理统计课程设计时学生反馈最高频、后果最严重的 5 个问题按“现象 → 原因 → 解决”结构整理每一条都来自真实翻车现场。4.1 现象PCA 降维后 GMM 聚类效果反而比原始数据差原因PCA 投影后数据维度降低但 GMM 的协方差矩阵Σ_k仍按降维后维度如 2D估计其行列式det(Σ_k)极小导致高斯密度计算溢出exp(-large_number)→ 0后验概率γ(z_{nk})全为 0EM 算法崩溃。解决在 GMM 的 M 步中对Σ_k添加更强正则项reg_covar如1e-3或改用“对角协方差”假设sigma_k[k] np.diag(np.diag(sigma_k[k]))牺牲部分建模能力换取数值稳定。验证时打印np.linalg.cond(sigma_k[k])条件数若 1e12 则必须正则化。4.2 现象LDA 投影后三类样本在二维空间严重重叠分类准确率低于 70%原因LDA 最大化tr(W^T S_b W)/tr(W^T S_w W)但当S_w接近奇异时如某类样本过少S_w^{-1}放大噪声投影方向被污染。Iris 中 Setosa 类在花瓣长度/宽度上天然分离但 Versicolour 和 Virginica 高度重叠LDA 可能过度优化后两者分离而牺牲 Setosa。解决计算S_w的奇异值分解U, s, Vh np.linalg.svd(S_w)若最小奇异值s[-1] 1e-8则用S_w_reg S_w 1e-5 * np.eye(p)替代或改用“分步 LDA”先用 PCA 降到 3D再在此子空间上运行 LDA规避原始 4D 的病态性。4.3 现象马氏距离分类器在某个类别上总是判错但 GMM 的后验概率显示该类置信度很高原因马氏距离d_M(x, μ_k)依赖Σ_k^{-1}而Σ_k是 GMM 估计的协方差。若某类如 Virginica样本在某个方向如花瓣宽度方差极小Σ_k对应元素接近 0其逆极大导致该方向微小扰动被放大距离失真。解决在计算马氏距离前对Σ_k进行“协方差缩放”sigma_k_scaled sigma_k[k] * np.max(np.diag(sigma_k[k]))即用最大方差归一化所有方差消除尺度差异。或改用“加权欧氏距离”d_weighted sqrt(Σ_j (x_j - μ_{kj})^2 / σ_{kjj}^2)其中σ_{kjj}是Σ_k对角线元素各维度方差。4.4 现象刀切法运行到第 87 次循环时抛出LinAlgError: Singular matrix原因第 87 个样本剔除后剩余 149 个样本中某类如 Versicolour只剩 49 个但S_w计算需至少p15个样本才能满秩p449 个样本虽够但若其分布高度线性相关如所有 Versicolour 样本的花萼长宽比恒为 1.5S_w仍可能秩亏。解决在manual_lda函数开头添加秩检查rank_sw np.linalg.matrix_rank(S_w)若rank_sw p则对S_w添加1e-4 * np.eye(p)并重算或在刀切法循环中捕获异常对该次迭代跳过并计为“无法评估”最后报告valid_iterations/150。4.5 现象不同降维方法PCA vs LDA的刀切法错误率相差不到 0.5%无法说明优劣原因只比较总错误率忽略了错误模式。PCA 可能在 Setosa 类全对但在另两类各错 5 个LDA 可能在 Setosa 错 2 个另两类各错 3 个——总错数相同但统计意义不同。解决输出混淆矩阵Confusion Matrix和 per-class 错误率。例如from sklearn.metrics import confusion_matrix; cm confusion本文还有配套的精品资源点击获取 
长河君 · 2026/10/19 6:34:15 · 阅读 1.2k 手写PCA/LDA/GMM与马氏距离：数理统计大作业原理实现指南
简介本资源是高校数理统计课程的大作业实践报告面向统计学、数据科学及相关专业本科生聚焦经典鸢尾花数据集的多方法综合分析。报告系统应用马氏距离度量样本相似性、混合高斯模型实现无监督聚类、主成分分析与线性判别分析进行降维与可分性提升并采用刀切法留一法评估不同降维策略与模型的泛化性能完整覆盖课程核心知识点与建模流程。资源为1个511KB的Word文档.docx内容结构严谨含摘要、引言、六大算法原理详述3.1–3.5节、四阶段数据处理分析4.1–4.4节及关键词与目录公式推导、代码逻辑说明与结果比对清晰可循。目前已有1426人学习下载适合课程复习、大作业参考、统计建模入门与算法对比研究。1. 这份数理统计大作业不是跑通 sklearn 就算交差它要你亲手推导马氏距离、手写 GMM EM 步骤、用刀切法硬刚每一种降维组合的泛化能力很多同学拿到“孙海燕老师数理统计大作业”这个标题第一反应是不就是拿 iris 数据集调个PCA().fit_transform()GaussianMixture().fit()吗三小时出图五小时写完摘要答辩前一晚改格式——结果被导师一句“你这个马氏距离是直接调scipy.spatial.distance.mahalanobis吗协方差矩阵求逆前检查过条件数没”当场卡住。这不是编程课设这是数理统计的原理验证型大作业它要求你把 PCA 的特征向量投影、LDA 的类间/类内散度矩阵、GMM 的 E-M 迭代、刀切法Jackknife的留一循环全部拆解到矩阵运算层面再用 Python 一行行实现。它面向的是某高校数理统计课程的真实教学目标——不考你会不会用 API而考你是否真正理解每个公式的物理意义、数值稳定性边界和统计假设前提。适合两类人一是正在啃《多元统计分析》教材、卡在 LDA 推导第 3 行的同学二是想用一份可答辩、可复现、带完整推导链的方案替代网上零散博客拼凑的毕设/课设学生。它解决的不是“怎么画出分类图”而是“当你的 PCA 降维后 GMM 聚类效果反而变差时你能否定位是协方差矩阵病态、还是 LDA 投影方向未对齐类别可分性、抑或刀切法暴露了模型过拟合”。2. 总体架构设计四层结构与五个核心模块的职责划分这份大作业不是单点算法堆砌而是一个闭环验证系统。我把它拆成四层结构数据层 → 降维层 → 建模层 → 验证层每一层都承担明确的统计任务且层间接口必须显式暴露数学含义。比如降维层输出的不能是黑盒数组而必须是带维度标签的pd.DataFrame其列名需体现“PC1解释方差比 72.9%”或“LDA1最大化类间离差比”建模层输入必须接收该结构并在拟合时校验输入维度是否匹配 GMM 成分数验证层则强制所有模块支持predict_proba和score接口为刀切法统一调用铺路。整个流程围绕五个核心模块展开它们不是孤立函数而是有严格依赖关系的组件数据预处理模块负责加载 UCI Iris 原始 CSV执行中心化非标准化因马氏距离对尺度敏感但需保证各变量同量纲、标签编码Setosa→0, Versicolour→1, Virginica→2并生成训练/验证分离的索引映射PCA 降维模块不调用sklearn.decomposition.PCA而是手动计算相关矩阵R (X^T X)/(n-1)调用np.linalg.eig(R)获取特征值/向量按累计解释方差比 ≥85% 自动截断主成分数通常为 2并返回正交投影矩阵W_pca ∈ ℝ^{4×2}LDA 降维模块手动构建类内散度矩阵S_w Σ_i Σ_{x∈C_i} (x - μ_i)(x - μ_i)^T与类间散度矩阵S_b Σ_i n_i (μ_i - μ)(μ_i - μ)^T求解广义特征值问题S_w^{-1} S_b v λv取前k个最大特征值对应特征向量组成W_lda ∈ ℝ^{4×2}GMM 建模模块基于降维后数据实现完整的 EM 算法E 步计算后验概率γ(z_{nk}) π_k N(x_n|μ_k,Σ_k) / Σ_j π_j N(x_n|μ_j,Σ_j)M 步更新π_k Σ_n γ(z_{nk})/N,μ_k Σ_n γ(z_{nk}) x_n / Σ_n γ(z_{nk}),Σ_k Σ_n γ(z_{nk}) (x_n - μ_k)(x_n - μ_k)^T / Σ_n γ(z_{nk})收敛判据为对数似然增量 1e-6刀切法验证模块对每个样本i临时剔除它用剩余n-1个样本重新训练 PCA/LDAGMM 全流程再用该模型预测x_i所属类别记录错误数最终计算总错判率err_rate k_total / 150。提示所有模块必须返回dict结构包含model模型参数字典、transform_matrix若适用、log_likelihoodGMM、error_count刀切法。这能避免答辩时被问“你 PCA 选了几个主成分依据是什么”时只答“我试了2个效果好”——你得拿出explained_variance_ratio_数组和累计曲线。2.1 数据预处理为什么中心化是必须的而标准化反而是陷阱Iris 数据四个属性单位不同厘米但马氏距离的核心前提是变量已处于同一量纲下且协方差矩阵可逆。直接标准化z-score会破坏原始数据的协方差结构——因为cov(z_i,z_j) cov((x_i-μ_i)/σ_i, (x_j-μ_j)/σ_j) cov(x_i,x_j)/(σ_i σ_j)这导致后续 PCA 特征向量不再反映原始变量间的相关性强度。正确做法是仅做中心化zero-mean保留原始尺度import numpy as np import pandas as pd # 加载原始 UCI Iris 数据150×5最后一列为字符串标签 df pd.read_csv(iris.data, headerNone, names[sepal_length,sepal_width,petal_length,petal_width,class]) X_raw df.iloc[:, :4].values # shape: (150, 4) y_raw df[class].map({Iris-setosa:0, Iris-versicolor:1, Iris-virginica:2}).values # 仅中心化减去均值不除标准差 X_centered X_raw - np.mean(X_raw, axis0) # shape: (150, 4) # 验证各列均值应接近0 print(Centered data column means:, np.mean(X_centered, axis0)) # [~0, ~0, ~0, ~0]这段代码的关键在于np.mean(X_raw, axis0)计算的是每列即每个属性的均值然后从整列中减去该值。它不改变变量间的协方差比例却为后续计算马氏距离d_M(x,y) sqrt((x-y)^T Σ^{-1} (x-y))提供了必要前提Σ必须是X_centered的协方差矩阵。如果误用标准化Σ会变成相关矩阵R此时马氏距离退化为欧氏距离因RI时d_M d_E完全失去度量多变量相关性的意义——这正是答辩时最容易被揪住的“概念性翻车点”。2.2 PCA 降维模块从相关矩阵到投影矩阵的完整推导链PCA 的本质是找一个正交基{u_1,...,u_k}使得数据在该基上的投影方差最大。教科书常写“求协方差矩阵C (X^T X)/(n-1)的特征向量”但对 Iris 这种小样本n150, p4直接算X^T X可能引入数值误差。更稳健的做法是先中心化X_c再计算相关矩阵R (X_c^T X_c)/(n-1)因各变量量纲一致R与C成比例特征向量相同。以下是手动实现 PCA 的核心步骤def manual_pca(X_centered, n_components2): n, p X_centered.shape # 1. 计算相关矩阵 R (p x p) R np.dot(X_centered.T, X_centered) / (n - 1) # 无偏估计 # 2. 求特征值分解R V diag(λ) V.T eigenvals, eigenvecs np.linalg.eig(R) # 3. 按特征值降序排列最大方差优先 idx np.argsort(eigenvals)[::-1] eigenvals eigenvals[idx] eigenvecs eigenvecs[:, idx] # 列向量为特征向量 # 4. 截取前 n_components 个主成分投影方向 W_pca eigenvecs[:, :n_components] # shape: (p, n_components) # 5. 投影数据X_pca X_centered W_pca X_pca np.dot(X_centered, W_pca) # 6. 计算解释方差比 explained_ratio eigenvals[:n_components] / np.sum(eigenvals) return { W_pca: W_pca, X_pca: X_pca, eigenvals: eigenvals, explained_ratio: explained_ratio, cumsum_ratio: np.cumsum(explained_ratio) } # 调用示例 pca_result manual_pca(X_centered, n_components2) print(Explained variance ratio for PC1 PC2:, pca_result[explained_ratio]) print(Cumulative sum:, pca_result[cumsum_ratio][-1]) # 通常 0.95这段代码的深层逻辑是R的第i个特征值λ_i表示数据在第i个主成分方向u_i上的投影方差u_i是R的第i个特征向量。因此W_pca的列就是我们要的正交基。注意np.linalg.eig返回的特征向量是列向量所以eigenvecs[:, idx]是正确的索引方式。若误写成eigenvecs[idx, :]得到的将是错误的行向量矩阵导致投影X_centered W_pca维度不匹配150×4 4×2 150×2这是初学者高频报错点。此外explained_ratio必须用eigenvals[:n_components] / np.sum(eigenvals)计算而非eigenvals[:n_components] / n_components——后者是常见误解混淆了“单个主成分方差”与“总方差”的关系。2.3 LDA 降维模块为什么 LDA 必须用带标签的类内/类间散度矩阵PCA 是无监督的它只看数据整体分布而 LDA 是有监督的它的目标是让同类样本投影后尽可能紧凑类内散度小异类样本投影后尽可能分离类间散度大。这决定了 LDA 的投影矩阵W_lda必须由标签y驱动。核心公式是最大化 Rayleigh 商J(W) tr(W^T S_b W) / tr(W^T S_w W)其中S_b是类间散度矩阵S_w是类内散度矩阵。手动实现的关键在于正确计算这两个矩阵def manual_lda(X_centered, y, n_components2): n, p X_centered.shape classes np.unique(y) n_classes len(classes) # 1. 计算全局均值 μ mu_global np.mean(X_centered, axis0) # shape: (p,) # 2. 初始化类内散度矩阵 S_w (p x p) 和类间散度矩阵 S_b (p x p) S_w np.zeros((p, p)) S_b np.zeros((p, p)) # 3. 对每个类别计算 S_w 和 S_b for cls in classes: X_cls X_centered[y cls] # 取出该类所有样本 n_cls X_cls.shape[0] mu_cls np.mean(X_cls, axis0) # 该类均值 # 类内散度S_w Σ (x - μ_cls)(x - μ_cls)^T diff_cls X_cls - mu_cls S_w np.dot(diff_cls.T, diff_cls) # 类间散度S_b n_cls * (μ_cls - μ_global)(μ_cls - μ_global)^T diff_global mu_cls - mu_global S_b n_cls * np.outer(diff_global, diff_global) # 4. 求解广义特征值问题 S_w^{-1} S_b v λv # 注意S_w 可能奇异需用 np.linalg.pinv 求伪逆 S_w_inv np.linalg.pinv(S_w) eigenvals, eigenvecs np.linalg.eig(np.dot(S_w_inv, S_b)) # 5. 按特征值降序排列越大表示类间分离越强 idx np.argsort(eigenvals)[::-1] eigenvals eigenvals[idx] eigenvecs eigenvecs[:, idx] # 6. 取前 n_components 个判别向量 W_lda eigenvecs[:, :n_components] # shape: (p, n_components) # 7. 投影数据 X_lda np.dot(X_centered, W_lda) return { W_lda: W_lda, X_lda: X_lda, eigenvals: eigenvals, S_w: S_w, S_b: S_b } # 调用示例 lda_result manual_lda(X_centered, y_raw, n_components2) print(LDA eigenvalues (largest first):, lda_result[eigenvals][:2])这里最易出错的是S_w和S_b的计算逻辑。S_w是所有类的类内离差平方和之和它衡量的是类内紧凑性S_b是各类均值相对于全局均值的加权离差平方和它衡量的是类间分离性。若漏掉n_cls权重如S_b np.outer(diff_global, diff_global)会导致S_b被低估LDA 方向偏向样本数少的类别。另一个坑是S_w可能奇异秩亏直接求逆np.linalg.inv(S_w)会报错必须用伪逆np.linalg.pinv。这也是为什么 LDA 在小样本、高维场景下不如 PCA 稳健——它对S_w的条件数极度敏感。3. 核心算法实现GMM 的 EM 迭代与马氏距离分类器的耦合设计GMM 不是黑箱它是对数据分布的显式建模。在 Iris 场景下我们假设每个类别c ∈ {0,1,2}对应一个高斯成分GMM 的目标是估计这三个成分的参数(π_c, μ_c, Σ_c)。但直接用sklearn.mixture.GaussianMixture会掩盖关键细节EM 算法的收敛性、协方差矩阵的约束全协方差对角球形、以及如何用马氏距离替代默认的predict_proba进行硬分类。本节将 GMM 与马氏距离深度耦合形成一条从概率建模到确定性判别的完整链路。3.1 GMM 的 EM 算法手写迭代拒绝黑盒EM 算法分为 EExpectation步和 MMaximization步。E 步计算每个样本属于每个成分的后验概率γ(z_{nk})M 步用这些后验概率加权更新参数。关键在于M 步更新的协方差矩阵Σ_k必须是正定的否则后续马氏距离计算会失败。因此我们在 M 步后添加协方差正则化加一个小的单位阵εIdef gmm_em(X, n_components3, max_iter100, tol1e-6, reg_covar1e-6): n, d X.shape # 初始化参数K-means 启发式 from sklearn.cluster import KMeans kmeans KMeans(n_clustersn_components, n_init10, random_state42) labels kmeans.fit_predict(X) pi_k np.array([np.sum(labelsk)/n for k in range(n_components)]) mu_k np.array([np.mean(X[labelsk], axis0) for k in range(n_components)]) # 初始协方差每个类的样本协方差 正则项 sigma_k np.array([ np.cov(X[labelsk].T) reg_covar * np.eye(d) for k in range(n_components) ]) log_likelihoods [] for iteration in range(max_iter): # E-step: 计算后验概率 γ(z_{nk}) gamma np.zeros((n, n_components)) for k in range(n_components): # 计算高斯概率密度 N(x_n | μ_k, Σ_k) diff X - mu_k[k] # (n, d) # Σ_k 是 dxd需求逆 try: sigma_inv np.linalg.inv(sigma_k[k]) except np.linalg.LinAlgError: # 若奇异加更大正则项 sigma_k[k] 1e-3 * np.eye(d) sigma_inv np.linalg.inv(sigma_k[k]) # 马氏距离平方: (x-μ)^T Σ^{-1} (x-μ) maha_sq np.sum(diff sigma_inv * diff, axis1) # (n,) # 高斯密度忽略归一化常数 norm_const np.sqrt(np.linalg.det(2*np.pi*sigma_k[k])) gauss_prob np.exp(-0.5 * maha_sq) / norm_const gamma[:, k] pi_k[k] * gauss_prob # 归一化 gamma使其每行和为1 gamma_sum np.sum(gamma, axis1, keepdimsTrue) gamma gamma / gamma_sum # M-step: 更新参数 N_k np.sum(gamma, axis0) # (n_components,) # 更新 π_k pi_k N_k / n # 更新 μ_k mu_k np.dot(gamma.T, X) / N_k[:, None] # (n_components, d) # 更新 Σ_k for k in range(n_components): diff X - mu_k[k] # (n, d) sigma_k[k] np.dot(gamma[:, k] * diff.T, diff) / N_k[k] # 添加正则化确保正定 sigma_k[k] reg_covar * np.eye(d) # 计算对数似然 log_likelihood np.sum(np.log(gamma_sum 1e-30)) log_likelihoods.append(log_likelihood) # 检查收敛 if iteration 0 and abs(log_likelihoods[-1] - log_likelihoods[-2]) tol: break return { pi_k: pi_k, mu_k: mu_k, sigma_k: sigma_k, gamma: gamma, log_likelihoods: log_likelihoods, n_iter: iteration 1 } # 示例在 PCA 降维后的数据上训练 GMM X_pca pca_result[X_pca] # (150, 2) gmm_result gmm_em(X_pca, n_components3, max_iter50) print(GMM converged in, gmm_result[n_iter], iterations) print(Final log-likelihood:, gmm_result[log_likelihoods][-1])这段代码的精妙之处在于它把 GMM 的核心计算E 步的后验概率、M 步的参数更新全部显式写出而非调用封装函数。特别是sigma_k的更新np.dot(gamma[:, k] * diff.T, diff) / N_k[k]是加权协方差的标准公式其中gamma[:, k]是每个样本属于第k类的置信度。reg_covar1e-6的正则项至关重要——它防止sigma_k奇异为后续马氏距离计算提供稳定基础。若跳过此步np.linalg.inv(sigma_k[k])在sigma_k[k]接近奇异时会崩溃这是 GMM 实现中最常见的“玄学翻车点”。3.2 马氏距离分类器用统计距离替代软概率直击判别本质GMM 默认的predict方法基于最大后验概率MAPargmax_k γ(z_{nk})。但这忽略了距离的几何意义。马氏距离d_M(x, μ_k) sqrt((x-μ_k)^T Σ_k^{-1} (x-μ_k))才是衡量样本x到第k类中心μ_k的“统计距离”——它考虑了该类的协方差结构Σ_k比欧氏距离更能反映真实分布形态。因此我们设计一个马氏距离分类器它不依赖γ(z_{nk})而是直接计算x到每个μ_k的马氏距离选择距离最小的类别def mahalanobis_classifier(X, mu_k, sigma_k): X: (n, d) 待分类样本 mu_k: (k, d) 各类均值 sigma_k: (k, d, d) 各类协方差矩阵 返回: (n,) 预测类别标签 n, d X.shape k mu_k.shape[0] distances np.zeros((n, k)) for i in range(k): try: sigma_inv np.linalg.inv(sigma_k[i]) except np.linalg.LinAlgError: # 协方差矩阵奇异时用伪逆 sigma_inv np.linalg.pinv(sigma_k[i]) diff X - mu_k[i] # (n, d) # 计算马氏距离平方(x-μ)^T Σ^{-1} (x-μ) # 使用 einsum 避免显式循环 maha_sq np.einsum(ij,jk,ik-i, diff, sigma_inv, diff) distances[:, i] np.sqrt(np.maximum(maha_sq, 0)) # 防止负数开方 # 返回距离最小的类别索引 return np.argmin(distances, axis1) # 在 PCA 降维数据上测试 y_pred_maha mahalanobis_classifier(X_pca, gmm_result[mu_k], gmm_result[sigma_k]) # 计算准确率 accuracy_maha np.mean(y_pred_maha y_raw) print(Mahalanobis classifier accuracy on PCA data:, accuracy_maha)这个分类器的设计哲学是GMM 建模的是概率但 Iris 分类的本质是判别。马氏距离直接量化了“样本离哪个类中心更近”且其计算过程复用了 GMM 已估计的μ_k和Σ_k无需额外训练。np.einsum(ij,jk,ik-i, diff, sigma_inv, diff)是高效计算马氏距离平方的向量化方法比双重 for 循环快一个数量级。np.maximum(maha_sq, 0)是安全防护——浮点误差可能导致maha_sq为极小负数开方会报nan。这种细节能让你的代码在答辩演示时稳如磐石。3.3 刀切法Jackknife验证逐样本剔除硬核检验泛化能力刀切法是比交叉验证更严格的模型评估方法尤其适合小样本n150。它的逻辑是每次只剔除一个样本用剩余 149 个样本重新训练整个 pipelinePCA/LDA → GMM → 马氏分类再用该模型预测被剔除的样本。150 次循环后统计总错判数。这种方法能暴露模型对单个样本的敏感性是检验“过拟合”的黄金标准。实现难点在于如何高效地“剔除一个样本”并重训我们采用索引掩码而非物理删除避免重复加载数据def jackknife_validation(X_centered, y,降维_func, gmm_func, classifier_func, n_components2): X_centered: (150, 4) 中心化数据 y: (150,) 标签 降维_func: 如 manual_pca 或 manual_lda接受 X 和 n_components gmm_func: 如 gmm_em接受降维后 X classifier_func: 如 mahalanobis_classifier接受降维后 X 和 GMM 参数 n len(y) errors 0 # 预分配存储避免循环中重复创建 X_pca_all np.zeros((n, n_components)) y_pred_all np.zeros(n, dtypeint) for i in range(n): # 创建掩码剔除第 i 个样本 mask np.ones(n, dtypebool) mask[i] False X_train X_centered[mask] # (149, 4) y_train y[mask] # (149,) # 1. 降维用训练集学习投影矩阵 dim_result 降维_func(X_train, n_componentsn_components) X_train_dim dim_result[X_pca] if X_pca in dim_result else dim_result[X_lda] # 2. GMM 训练 gmm_result gmm_func(X_train_dim, n_components3) # 3. 对被剔除的样本 x_i 进行预测 x_test X_centered[i:i1] # (1, 4) x_test_dim np.dot(x_test, dim_result[W_pca]) if W_pca in dim_result else np.dot(x_test, dim_result[W_lda]) y_pred classifier_func(x_test_dim, gmm_result[mu_k], gmm_result[sigma_k]) if y_pred[0] ! y[i]: errors 1 error_rate errors / n return { error_count: errors, error_rate: error_rate, accuracy: 1 - error_rate } # 执行 PCA GMM 马氏距离的刀切法 pca_jackknife jackknife_validation( X_centered, y_raw, manual_pca, gmm_em, mahalanobis_classifier, n_components2 ) print(PCAGMMMahalanobis Jackknife error rate:, pca_jackknife[error_rate]) # 执行 LDA GMM 马氏距离的刀切法 lda_jackknife jackknife_validation( X_centered, y_raw, manual_lda, gmm_em, mahalanobis_classifier, n_components2 ) print(LDAGMMMahalanobis Jackknife error rate:, lda_jackknife[error_rate])这段代码的工程价值在于它把刀切法封装成一个可插拔的验证器。只需传入不同的降维_funcmanual_pca或manual_lda就能对比两种降维策略的泛化能力。注意mask的使用——它通过布尔索引高效提取子集比np.delete(X_centered, i, axis0)快得多。x_test_dim的计算必须复用dim_result中的W_pca或W_lda确保测试样本与训练样本在同一个投影空间。若误用X_centered[i]直接降维会得到错误结果。刀切法耗时较长150×训练但这是值得的——它给出的error_rate是模型真实鲁棒性的铁证。4. 关键避坑指南5 个让答辩老师眼前一亮的“血泪经验”做这份大作业踩坑是常态。但有些坑踩一次就足以让答辩陷入被动。以下是我在指导某高校数理统计课程设计时学生反馈最高频、后果最严重的 5 个问题按“现象 → 原因 → 解决”结构整理每一条都来自真实翻车现场。4.1 现象PCA 降维后 GMM 聚类效果反而比原始数据差原因PCA 投影后数据维度降低但 GMM 的协方差矩阵Σ_k仍按降维后维度如 2D估计其行列式det(Σ_k)极小导致高斯密度计算溢出exp(-large_number)→ 0后验概率γ(z_{nk})全为 0EM 算法崩溃。解决在 GMM 的 M 步中对Σ_k添加更强正则项reg_covar如1e-3或改用“对角协方差”假设sigma_k[k] np.diag(np.diag(sigma_k[k]))牺牲部分建模能力换取数值稳定。验证时打印np.linalg.cond(sigma_k[k])条件数若 1e12 则必须正则化。4.2 现象LDA 投影后三类样本在二维空间严重重叠分类准确率低于 70%原因LDA 最大化tr(W^T S_b W)/tr(W^T S_w W)但当S_w接近奇异时如某类样本过少S_w^{-1}放大噪声投影方向被污染。Iris 中 Setosa 类在花瓣长度/宽度上天然分离但 Versicolour 和 Virginica 高度重叠LDA 可能过度优化后两者分离而牺牲 Setosa。解决计算S_w的奇异值分解U, s, Vh np.linalg.svd(S_w)若最小奇异值s[-1] 1e-8则用S_w_reg S_w 1e-5 * np.eye(p)替代或改用“分步 LDA”先用 PCA 降到 3D再在此子空间上运行 LDA规避原始 4D 的病态性。4.3 现象马氏距离分类器在某个类别上总是判错但 GMM 的后验概率显示该类置信度很高原因马氏距离d_M(x, μ_k)依赖Σ_k^{-1}而Σ_k是 GMM 估计的协方差。若某类如 Virginica样本在某个方向如花瓣宽度方差极小Σ_k对应元素接近 0其逆极大导致该方向微小扰动被放大距离失真。解决在计算马氏距离前对Σ_k进行“协方差缩放”sigma_k_scaled sigma_k[k] * np.max(np.diag(sigma_k[k]))即用最大方差归一化所有方差消除尺度差异。或改用“加权欧氏距离”d_weighted sqrt(Σ_j (x_j - μ_{kj})^2 / σ_{kjj}^2)其中σ_{kjj}是Σ_k对角线元素各维度方差。4.4 现象刀切法运行到第 87 次循环时抛出LinAlgError: Singular matrix原因第 87 个样本剔除后剩余 149 个样本中某类如 Versicolour只剩 49 个但S_w计算需至少p15个样本才能满秩p449 个样本虽够但若其分布高度线性相关如所有 Versicolour 样本的花萼长宽比恒为 1.5S_w仍可能秩亏。解决在manual_lda函数开头添加秩检查rank_sw np.linalg.matrix_rank(S_w)若rank_sw p则对S_w添加1e-4 * np.eye(p)并重算或在刀切法循环中捕获异常对该次迭代跳过并计为“无法评估”最后报告valid_iterations/150。4.5 现象不同降维方法PCA vs LDA的刀切法错误率相差不到 0.5%无法说明优劣原因只比较总错误率忽略了错误模式。PCA 可能在 Setosa 类全对但在另两类各错 5 个LDA 可能在 Setosa 错 2 个另两类各错 3 个——总错数相同但统计意义不同。解决输出混淆矩阵Confusion Matrix和 per-class 错误率。例如from sklearn.metrics import confusion_matrix; cm confusion本文还有配套的精品资源点击获取