简介这份数理统计课程大作业以鸢尾花数据集为分析对象面向正在学习多元统计、模式识别或准备课程报告的高校学生围绕花萼与花瓣四项属性完成分类、聚类与降维的完整实践。文档系统梳理马氏距离、高斯混合模型、主成分分析、线性判别分析及刀切法的原理并给出数据预处理、算法比对与结果讨论其中重点展示PCA与LDA降维后在极大似然估计下拟合混合高斯模型、再用留一法评估降维程度与模型差异的思路。资源包仅含1个docx文件大小约511KB结构包含摘要、引言、符号说明、算法原理、数据处理分析和结论等章节便于按目录查阅与引用。目前已有1429人学习适合作为课程作业参考、统计方法复习或机器学习入门的案例材料。1. 从一份看起来很水的数据集作业说起数理统计大作业挑鸢尾花数据集十个里有八个觉得这题被做烂了150 行、4 列、3 个类别随便调个库就能出结果。真动手才发现难点根本不在分类准确率而在于怎么把马氏距离、混合高斯模型、主成分分析、线性判别分析、刀切法这五个知识点串成一条自洽的链路——降维降几维、GMM 用几个成分、刀切法怎么比对两种降维方法的差异每一步都要有理由。这套方案就是围绕这条链路展开的从 4 维原始特征出发用 PCA 和 LDA 分别降到低维空间用高斯混合模型在极大似然意义下拟合分布用马氏距离判类最后用留一法给出可信的错误率对比。适合正在写数理统计大作业、统计学习课程设计的人直接照着改。2. 数据与符号150 个样本、4 个属性、3 个类别的落表方式鸢尾花数据集的结构已经被反复讲过但真到写报告那一章多数人卡在符号怎么定义、表怎么列、预处理做到哪一步算够。这一章把这个地基打牢后面三章的推导和代码都建立在这套符号上。2.1 数据集结构与符号约定数据集包含 150 个观测样本分属 3 个类别每类 50 个每个样本有 4 个数值型属性单位统一为厘米。原始数据在公开的机器学习数据仓库里可以拿到字段顺序常见的是花萼长度、花萼宽度、花瓣长度、花瓣宽度但不同来源的列顺序并不一致导入后第一件事是对着前几行核对列名别默认列序。符号体系建议在报告开头一次性定死后面全部沿用否则写到第四章会出现同一个量两个名字的尴尬。我一般用下面这张表符号含义说明n样本个数取 150x₁花萼长度连续型单位 cmx₂花萼宽度连续型单位 cmx₃花瓣长度连续型单位 cmx₄花瓣宽度连续型单位 cmy类别标签0 / 1 / 2仅作标识D样本矩阵n × 4μ均值向量4 维Σ协方差矩阵4 × 4R相关矩阵4 × 4标准化后等价于协方差阵u线性变换方向4 × 1 或 4 × kK高斯成分个数每类内独立设定标签编码有一个容易忽略的点0、1、2 只是类别代号不构成大小关系。后面所有涉及标签的运算只有同类/异类判断绝不能把标签当数值参与回归或者求均值否则会在 LDA 的类间散度计算里引入完全不存在的序关系。四个属性的取值范围差异不小。花瓣长度大致落在 1 到 7 厘米花瓣宽度只有 0.1 到 2.5 厘米花萼宽度的波动更窄。这种量级差异对欧氏距离和基于协方差矩阵的 PCA 都是干扰项后面会专门处理。2.2 数据清洗与标准化为什么先换成相关矩阵先做三件最基础的检查缺失值、重复行、类别计数。这组数据通常是干净的150 行无缺失但重复行值得看一眼——如果发现有两条四个属性完全相同的记录却标了不同类别那就是天然的噪声样本留一法里必然被判错一次报告里要能解释。类别计数确认 50/50/50后面按类分层抽样才有依据。真正需要设计的是标准化策略。PCA 的目标是找方差最大的方向如果直接对原始协方差矩阵做特征分解方差绝对值最大的花瓣长度会独占第一主成分另外三个属性几乎没有话语权降维结果就变成了只看花瓣长度。消除这种影响的常规做法是对每个属性做 z 分数标准化等价地用相关矩阵 R 代替协方差矩阵 Σ 做特征分解import numpy as np def corr_matrix(X): 按列标准化后求相关矩阵等价于对 X 做 z-score 再求协方差 mu X.mean(axis0) sd X.std(axis0, ddof1) # 样本标准差用 n-1 Z (X - mu) / sd # 每列均值 0、方差 1 R np.corrcoef(Z, rowvarFalse) # 4 x 4 相关矩阵 return mu, sd, Rddof1是无偏估计的常规选择样本量 150 时和ddof0差别很小但报告里写清楚更稳妥。rowvarFalse表示每一列是一个变量这一行写错会让相关矩阵变成 150 × 150后面特征分解直接炸掉是新手最常见的翻车点之一。标准化之后还有一个判断要做四个属性之间到底存不存在值得降维的相关性。把 R 打印出来看非对角元素如果花瓣长度和花瓣宽度的相关系数在 0.9 附近花萼长度和花瓣长度也在 0.8 以上说明信息确实重叠PCA 有存在价值如果全是 0.1 以下那降维就是纯粹丢信息这时候应该老实报告该数据在相关意义下近似正交而不是硬降。这组数据的实际相关结构属于前者花瓣两个属性高度相关这也是二维降维后仍能保留大部分信息的原因。2.3 训练/验证的划分方式与样本组织这份作业的评估方式不能用简单的训练集/测试集对半切。样本只有 150 条切一半出去每类只剩 25 条训练样本GMM 的协方差矩阵在 2 维空间还算稳4 维空间下每类 25 条估计 4 × 4 协方差阵就有点勉强误差估计的抖动会非常大。刀切法留一法在这里是更合适的选择每次留出 1 条做测试其余 149 条训练重复 150 次。好处是几乎用上了全部数据偏差小代价是要训练 150 次模型PCA、LDA、GMM 各跑一遍Python 里循环 150 次完全在可接受范围内单次实验通常几十秒到几分钟。如果打算做多组参数对比比如 PCA 降到 1、2、3 维LDA 降到 1、2 维总训练次数会乘上组数这时候建议把降维部分预先算好——PCA 和 LDA 的变换方向只依赖训练集但严格来说每次留一都要重新拟合偷懒复用会引入轻微的信息泄漏报告里要如实说明是哪种做法。样本组织上还有个细节留一法循环里要保证每类的样本都被均匀留出。直接按行号 0 到 149 顺序循环是没问题的因为原始数据通常是按类别排好序的前 50 条全是第一类。如果做过 shuffle记得固定随机种子让结果可复现。3. 算法原理与选型五个模块各自解决什么问题五个算法不是并列关系而是一条流水线上的五个工位两个降维工位PCA、LDA、一个建模工位GMM、一个判类工位马氏距离、一个质检工位刀切法。搞清楚每个工位的输入输出报告的逻辑链条才不会散。3.1 马氏距离把协方差结构带进距离里欧氏距离的问题在于它假设各维度独立且等方差四个属性量纲虽然都是厘米但方差差了一个量级直接用欧氏距离等于默认花瓣宽度的 1 厘米和花瓣长度的 1 厘米等价显然不合理。马氏距离的定义是d²(x, μ) (x − μ)ᵀ Σ⁻¹ (x − μ)其中 μ 是均值向量Σ 是协方差矩阵。三个性质值得在报告里点出来当 Σ 取单位阵时退化为欧氏距离它等价于先做标准化再做欧氏距离它是无量纲的不受量纲变化影响。第三点在跨数据集比较时特别有用。实际计算里 Σ 可能接近奇异尤其是在留一法某次训练中某类样本数偏少、或者降维后的维度接近样本数时。稳妥写法是用伪逆np.linalg.pinv或者给对角线加一个极小的正则项如 1e-6避免求逆时报错或者数值爆炸。这一步不做程序会在某次留一循环里突然抛LinAlgError排查起来很费时间。在这套方案里马氏距离承担的是判类职责GMM 拟合完之后每个高斯成分有自己的一组 (μₖ, Σₖ)新样本对每个成分算一次马氏距离取距离最小的成分所属类别作为预测结果。这个做法在高斯假设成立、各类先验相近时和最大后验判别是等价的。3.2 混合高斯模型用 EM 求极大似然单个高斯分布只能描述一个钟形的簇而一个类别内部的样本分布未必这么规整。高斯混合模型把若干个高斯成分按权重加起来p(x) Σₖ πₖ · N(x | μₖ, Σₖ)其中 πₖ 是混合权重满足求和为 1。对 n 个样本写对数似然直接对 π、μ、Σ 求导会得到一组互相耦合的方程没有闭式解所以走 EM 迭代E 步用当前参数算每个样本属于每个成分的后验概率责任度M 步用责任度加权更新 π、μ、Σ交替直到对数似然增量小于阈值。有几个设计决策要在报告里交代清楚。第一是成分数 K 怎么定每类内部用一个成分K1退化为单高斯是最简方案也可以用 BIC 在 K1 到 3 之间选BIC 会惩罚参数过多。第二是协方差类型full 最灵活但参数最多d 维下 d(d1)/2 个diag 只估方差不管相关性tied 让所有成分共享一个协方差阵。这组数据每类 50 个样本2 维下 full 只有 5 个参数完全撑得住4 维下 10 个参数也还行但留一法里每类只剩 49 条稳定性要盯一下。第三是初始化EM 对初值敏感常规做法是用 k-means 先粗分一次拿均值或者跑 10 次不同随机初值取对数似然最大的那组后者贵但省心。3.3 主成分分析无监督降维的上限在哪PCA 的思路是把 p 个相关变量线性变换成 p 个互不相关的新变量新变量的方差依次递减。方差大意味着信息量大所以保留前 k 个就完成了降维。数学上变换方向就是相关矩阵或协方差矩阵的前 k 个最大特征值对应的单位特征向量。维数怎么选看累积贡献率把特征值归一化后累加通常取到 85% 或 90% 就够。这组数据前两个主成分的累积贡献率一般在 95% 以上所以降到 2 维是很自然的选择降到 1 维会损失较多信息降到 3 维则收益很小、降维的意义被削弱。PCA 的局限必须点明它完全无视标签找的是数据整体最分散的方向而这个方向未必是最能把三个类别分开的方向。极端情况下方差最大的方向恰好是类内波动最大的方向降维之后三类样本反而糊在一起。这不是实现 bug是方法本身的属性也是后面引入 LDA 的直接动机。3.4 线性判别分析用上标签后的降维差异LDA 的目标和 PCA 相反不是让总体方差最大而是让投影之后同类样本尽量聚拢、异类样本尽量远离。定义类内散度矩阵 S_w 和类间散度矩阵 S_b投影方向 u 的目标函数是J(u) uᵀ S_b u / uᵀ S_w u最大化这个比值解归到广义特征值问题 S_w⁻¹ S_b u λu取最大若干特征值对应的特征向量组成变换矩阵。因为 S_b 的秩最多是 c − 1c 为类别数LDA 最多只能给出 c − 1 个有效判别方向3 类数据最多降到 2 维——这个硬上限是 PCA 没有的报告里要专门说明否则读者会问为什么 LDA 不降到 3 维。S_w 在样本数少于维度时容易奇异伪逆是常规兜底。另外 LDA 是监督方法留一法循环里每次只能用 149 条训练样本的标签重新拟合变换矩阵如果把全量数据的标签信息带进降维步骤错误率会被系统性地低估这是整个方案里最隐蔽的一处泄漏。3.5 刀切法留一法怎么给出可信的错误率刀切法的步骤可以拆成五步从某一类的 n₁ 个训练样本中剔除一个用剩余样本建立判别规则用这个规则判别被剔除的那个样本重复直到该类每个样本都被判别过一次记判错数 k₁对每个类别重复上述过程得到 kᵢ最后算总判错率。它比回报法用训练集自己测自己靠谱得多回报法在样本少的时候会给出虚高的准确率也比留出法稳定因为没有随机切分带来的方差。代价是训练次数等于样本数计算量线性增长。150 条样本、单次训练几十毫秒的量级总耗时完全可控。真正需要注意的是错误率估计本身的方差。150 个样本判错 3 个是 2%判错 4 个是 2.67%两个方法差 1 个样本就会让排名翻转。所以报告里不能只写LDA 比 PCA 好要写清楚差了几个样本、这个差异在 150 的样本量下是否站得住。保守的做法是除了错误率再报一个混淆矩阵看错分集中在哪两类之间——这组数据里第一类通常完全可分错分基本都发生在后两类的交界带。4. 代码实现PCA/LDA 类、GMM 的 EM 与刀切法主循环算法原理写清楚之后实现部分的核心是把每个模块封装成职责单一的类或函数让刀切法主循环只负责调度不掺和具体计算。这样调试时能单独验证每个环节报告里的 UML 图也有实际内容可画。4.1 PCA 类fit 与 transform 的职责边界PCA 用类实现fit只学参数均值、标准差、投影矩阵transform只做映射。这个边界划清楚留一法循环里才不会写错import numpy as np class PCA: def __init__(self, n_components2): self.n_components n_components def fit(self, X): # 1) 记录训练集均值与标准差供 transform 复用 self.mu_ X.mean(axis0) self.sd_ X.std(axis0, ddof1) Z (X - self.mu_) / self.sd_ # 标准化消除量纲 # 2) 相关矩阵做特征分解对称矩阵用 eigh结果按升序 R np.corrcoef(Z, rowvarFalse) eigval, eigvec np.linalg.eigh(R) order np.argsort(eigval)[::-1] # 改成降序 self.eigval_ eigval[order] self.eigvec_ eigvec[:, order] # 3) 取前 k 个方向记录累积贡献率 self.components_ self.eigvec_[:, :self.n_components] self.explained_ratio_ self.eigval_ / self.eigval_.sum() return self def transform(self, X): # 必须用训练集的 mu_ 和 sd_不能用测试集自己的统计量 Z (X - self.mu_) / self.sd_ return Z self.components_三个参数的含义n_components是保留的主成分个数做降维程度对比时依次取 1、2、3explained_ratio_是各主成分的方差贡献率报告里那张降维到 k 维保留多少信息的表直接从这取components_是 4 × k 的投影矩阵每一列是一个方向。transform里用训练集的mu_和sd_是关键。如果偷懒用X自己的均值方差训练集和测试集的标准化基准就不一致留一法每一折都引入一次系统性偏差最后错误率会莫名其妙地好看。这个坑我在第一次写的时候踩过当时怎么调参都解释不了为什么测试准确率比训练准确率还高。4.2 LDA 类广义特征值问题的数值稳定写法LDA 的实现在散度矩阵累加和广义特征值求解两步class LDA: def __init__(self, n_components2): self.n_components n_components def fit(self, X, y): d X.shape[1] classes np.unique(y) mu_all X.mean(axis0) S_w np.zeros((d, d)) S_b np.zeros((d, d)) for c in classes: Xc X[y c] mu_c Xc.mean(axis0) diff Xc - mu_c S_w diff.T diff # 类内散度 m (mu_c - mu_all).reshape(-1, 1) S_b len(Xc) * (m m.T) # 类间散度按类样本数加权 # 广义特征值问题pinv(S_w) S_b 的特征分解 A np.linalg.pinv(S_w) S_b eigval, eigvec np.linalg.eig(A) order np.argsort(eigval.real)[::-1] # 有效方向数上限为 类别数 - 1 k min(self.n_components, len(classes) - 1) self.components_ eigvec[:, order][:, :k].real return self def transform(self, X): return X self.components_np.linalg.eig处理的是非对称矩阵返回值可能是复数取.real是必要的否则后面算距离会报类型错误。pinv代替inv是为了应对 S_w 奇异3 类数据在 2 维空间下 S_w 通常是满秩的但留一法某次训练里如果某类样本数掉到 2 条以下秩就掉了稳一点没坏处。len(classes) - 1这个截断不能省。如果不截断直接取前n_components个特征向量当n_components3时第三个方向的对应特征值接近 0投影结果是纯噪声还会让后续 GMM 的协方差估计出现数值问题。4.3 GMME 步、M 步与参数初始化GMM 的核心是两个函数交替调用。为了代码短一些这里给 E 步和 M 步的骨架高斯密度用 SciPy 的实现from scipy.stats import multivariate_normal def e_step(X, weights, means, covs): n, k X.shape[0], len(weights) resp np.zeros((n, k)) for j in range(k): # 每个成分的加权密度 resp[:, j] weights[j] * multivariate_normal.pdf( X, meanmeans[j], covcovs[j], allow_singularTrue) resp / resp.sum(axis1, keepdimsTrue) # 归一化成责任度 return resp def m_step(X, resp): n, k resp.shape nk resp.sum(axis0) # 各成分的有效样本数 weights nk / n means (resp.T X) / nk[:, None] covs [] for j in range(k): diff X - means[j] covs.append((resp[:, j][:, None] * diff).T diff / nk[j]) return weights, means, covsallow_singularTrue是给留一法兜底的某次训练里成分的有效样本数nk[j]可能极小协方差阵接近奇异不加这个参数会直接抛异常。代价是密度值可能异常大所以初始化阶段要尽量避免产生空成分——用 k-means 初始化均值再把权重平均分比纯随机初始化稳得多。收敛判据用对数似然增量阈值设 1e-6 就够通常 20 到 50 轮内收敛。迭代次数上限设 200 轮防止数值问题导致死循环。每个类别单独拟合一个 GMM预测时对测试样本算它到每个成分的马氏距离取最小者所属类别。如果每类只用一个成分K1那就退化成最近高斯中心 协方差加权实现更简单报告里也能和 K1 的结果做个对比。4.4 刀切法主循环把降维与分类串起来主循环负责调度把降维、拟合、预测三步串起来并且统计每类的判错数def loo_evaluate(X, y, reducer_cls, reducer_kwargs, gmm_k1): n len(y) wrong 0 conf np.zeros((3, 3), dtypeint) # 混淆矩阵 for i in range(n): idx np.arange(n) ! i X_tr, y_tr X[idx], y[idx] X_te X[i:i1] # 1) 只用训练折拟合降维器避免标签泄漏 reducer reducer_cls(**reducer_kwargs).fit(X_tr, y_tr) \ if reducer_cls is LDA else reducer_cls(**reducer_kwargs).fit(X_tr) Z_tr, Z_te reducer.transform(X_tr), reducer.transform(X_te) # 2) 每类单独建 GMM返回该类的成分参数列表 models {c: fit_gmm(Z_tr[y_tr c], gmm_k) for c in np.unique(y_tr)} # 3) 马氏距离判类 pred predict_by_mahalanobis(Z_te[0], models) conf[y[i], pred] 1 wrong int(pred ! y[i]) return 1 - wrong / n, confreducer_cls is LDA那个三元表达式是为了兼容两种fit签名实际写的时候拆成两个函数更清楚。核心约束只有一条fit只看训练折transform才能碰测试样本。混淆矩阵conf的行是真实类别、列是预测类别答辩时这张表比一个准确率数字有说服力得多——它能直接告诉你错分发生在哪两类之间。跑对比实验时把reducer_cls换成 PCA 或 LDAreducer_kwargs里的n_components依次取 1、2、3LDA 会自动截断到 2gmm_k取 1 和 2就能得到一张完整的对比矩阵。每组配置跑 150 次留一一共十几组总耗时通常在几分钟量级。5. 避坑与排查五个最容易翻车的细节这套流程跑通不难跑出可信的结果难。下面五条是我在复现和改这份作业时踩过或者见别人踩过的按现象 → 原因 → 解决写对照着自查能省不少时间。5.1 降维环节的三个翻车点降维用了全量数据错误率虚低。现象是留一法跑出来的准确率比预期高出一大截比如 98% 以上换个随机种子结果也没变化。原因是 PCA 或 LDA 的变换矩阵在循环外一次性拟合了全部 150 条样本测试样本的标签信息已经提前渗进降维方向等于开卷考试。解决方式是把fit挪进循环体每折重新拟合如果为了提速确实要预计算报告里必须明确写降维方向在全体样本上估计存在轻微信息泄漏别装作没这回事。PCA 用了协方差矩阵而不是相关矩阵。现象是第一主成分的载荷几乎全部集中在花瓣长度上降维之后另外三个属性形同虚设。原因是四个属性方差差异大未标准化时方差绝对值主导了特征分解。解决办法是先 z-score 标准化再求协方差或者直接用相关矩阵——两者在数学上等价代码里选一种写清楚即可。LDA 降维维数超过类别数减一。现象是n_components3时程序不报错但第三个方向投影出来的数值极小GMM 拟合时协方差阵接近奇异偶尔抛LinAlgError。原因是 S_b 的秩最多为 c − 1 2。解决办法是在fit里显式截断k min(n_components, len(classes) - 1)并在报告里说明这个理论上限这本身就是个能加分的知识点。5.2 模型与评估环节的两个翻车点EM 初始化不当结果每次跑都不一样。现象是同样的数据、同样的 K两次运行的对数似然差一大截错误率上下浮动好几个百分点。原因是随机初始化落到了局部极值或者某个成分一开始就没有样本认领变成了空成分。解决办法是用 k-means 的聚类中心初始化均值、权重设为均分或者跑 10 次随机初值取对数似然最大的那组把随机种子固定下来写进代码。只报准确率不看混淆矩阵。现象是报告里写某方法准确率 96%答辩时被问错的那 6 个样本错在哪就答不上来。原因是单一指标掩盖了错误分布而这组数据的三类难度差异极大——第一类几乎线性可分后两类在边界上有重叠。解决办法是把混淆矩阵一并输出指出错分集中在哪两类之间再解释这是因为这两个类别的花瓣尺寸分布本身就存在交叠属于数据固有难度而非方法失效。这一句话往往比多跑十组实验更能体现你理解了这个模型。6. 验收与答辩用一张比对表证明方法差异是真的这套作业最后要交的不是代码是一个能自圆其说的结论。验收的核心动作只有一件把所有配置的结果压到一张表里让差异自己说话。我一般会做这样一张对比表行是降维方法加维数列是 GMM 成分数和两个指标降维方法维数成分数 K留一错误率主要错分区间不降维41待测后两类PCA11待测后两类明显增多PCA21待测后两类PCA22待测后两类PCA31待测后两类LDA11待测后两类LDA21待测后两类LDA22待测后两类表里的待测要你自己跑出来填我不给具体数字因为随机种子、初始化方式、是否分层都会影响小数点后那一位。填完之后按三条线读第一LDA 在同等维数下通常优于 PCA原因是它用上了标签信息这正是理论预期如果实测反了先查泄漏和实现 bug别急着写与理论不符第二维数从 1 涨到 2 收益明显从 2 涨到 3 收益很小说明二维已经抓住了主要判别信息第三K 从 1 涨到 2 未必带来提升样本量小时参数增多反而过拟合这个结论本身就可以写进报告。答辩被问到的概率最高的三个问题提前准备好答案。一是为什么 LDA 最多降到 2 维答 S_b 的秩上限是 c − 1。二是刀切法和交叉验证有什么区别答留一法是 k 折交叉验证在 k n 时的特例偏差最小但方差较大样本量小时更适用。三是错误率差一个样本算不算有差异这个问题最考验功底老实回答150 个样本下 2% 和 2.67% 的差距只有 1 个样本不足以支撑强结论只能说方向一致比硬吹方法优劣得分更高。有个习惯我从做这份作业之后一直保留任何涉及方法 A 比方法 B 好的结论都要先跑三遍不同随机种子看排名稳不稳再决定这句话能不能写进报告。数理统计这门课考的从来不是调包能力是你能不能在噪声里区分出信号。这套鸢尾花分析方案的完整推导、代码和对比表我整理成了一份文档需要照着改的话可以直接拿去对照。希望帮到你。本文还有配套的精品资源点击获取