简介：这份数理统计大作业文档面向高校统计学、机器学习课程的学习者，围绕经典鸢尾花数据集展开完整分析，帮助读者理解多重变量数据的建模流程。内容涵盖马氏距离度量样本相似度、混合高斯模型聚类、主成分分析降维去噪、线性判别分析兼顾类间与类内差异，以及刀切法评估降维方法与模型泛化能力，并附有数据预处理与结果讨论。资源包共1个docx文件，约511KB，结构完整、论述连贯，适合作为课程作业参考或复习数理统计知识点的范本。目前已有1429人学习下载，读者可从中获得从算法原理到数据处理分析的完整思路，对照梳理各方法的适用场景与实现要点，快速搭建自己的分析框架。

1. 从一份鸢尾花大作业说起：为什么降维和聚类要放在同一套流程里

数理统计大作业最容易翻车的地方，不是公式不会推，而是把 PCA、LDA、GMM 各写一个 demo，最后拼不出一条完整的分析链。鸢尾花数据集只有 150 个样本、4 个特征、3 个类别，看起来简单，但它恰好能把“降维—建模—评估”这条主线完整跑一遍：先用主成分分析（PCA）和线性判别分析（LDA）把四维数据压到低维，再在高斯混合模型（GMM）下做极大似然估计，最后用刀切法（留一法）比较不同降维程度和方法之间的差异。这套方案适合正在做数理统计课程设计、机器学习毕设开题，或者需要交一份可复现分析报告的人。下面按需求、设计、实现、避坑、验收的顺序拆开讲，代码和参数都能直接改。

2. 需求与数据底座：鸢尾花数据集到底要解决什么问题

2.1 四个特征、三个类别与一条分析主线

鸢尾花数据集由 Fisher 在 1936 年整理，1988 年公开，包含 150 个观测，每类 50 个，特征是花萼长度、花萼宽度、花瓣长度、花瓣宽度，单位都是厘米，标签是 Setosa、Versicolour、Virginica 三类。这个数据集的价值不在于“分类准确率刷到多高”，而在于它同时满足几个条件：特征都是数值型、量纲一致、近似服从高斯分布、特征之间存在明显线性相关。这几点决定了后面可以放心用马氏距离、GMM 和 PCA。

大作业要回答的问题可以拆成三层。第一层是降维：四个特征里花瓣长度和花瓣宽度高度相关，花萼长度和花萼宽度也有一定相关性，直接建模会引入冗余和噪声，需要 PCA 或 LDA 压到一维或二维。第二层是建模：降维后的数据用 GMM 拟合分布，在极大似然意义下估计均值、协方差和混合权重，再用样本到各高斯成分均值的马氏距离判定类别。第三层是评估：降维到几维、用哪种降维方法，不能拍脑袋，要用刀切法算判错率来比。

我一般会把整份作业的符号先固定下来，避免后面公式和代码对不上。样本矩阵记为 $D$，样本数 $n=150$，输入向量 $x=(x_1,x_2,x_3,x_4)^T$ 对应四个特征，类别 $y\in{0,1,2}$，均值向量 $\mu$，协方差矩阵 $\Sigma$，相关矩阵 $R$，线性变换方向 $u$。符号表看着是形式主义，但写到第四章比对结果时，你会发现没有它根本说不清“第几主成分对应哪个特征”。

2.2 数据读取与预处理的工程化写法

原始数据从公开机器学习仓库获取，格式是逗号分隔、无表头、最后一列是类别。读取时不要手写 150 行，直接用 pandas 加列名，顺手把类别映射成 0/1/2。预处理有两个关键动作：一是检查缺失值和重复值，鸢尾花数据一般干净，但作业里要写出检查过程；二是标准化，PCA 用相关矩阵代替协方差矩阵来消除量纲影响，代码里体现为对特征做 z-score。

import pandas as pd import numpy as np COLUMNS = ["sepal_len", "sepal_wid", "petal_len", "petal_wid", "label"] LABEL_MAP = {"Iris-setosa": 0, "Iris-versicolor": 1, "Iris-virginica": 2} def load_iris(path): df = pd.read_csv(path, header=None, names=COLUMNS) # 类别字符串映射为 0/1/2，方便后续按整数索引 df["label"] = df["label"].str.strip().map(LABEL_MAP) assert df["label"].notna().all(), "存在未识别的类别标签" assert df.shape[0] == 150, "样本数应为 150" return df def standardize(X): # 用总体标准差做 z-score，等价于用相关矩阵做 PCA mu = X.mean(axis=0) sigma = X.std(axis=0, ddof=0) return (X - mu) / sigma, mu, sigma

这段代码里ddof=0是有意为之：PCA 用相关矩阵时，标准化要除以总体标准差，和协方差矩阵的特征分解保持口径一致。如果这里用ddof=1，特征值会有细微偏差，虽然对鸢尾花这种量级影响不大，但答辩时被问到“你的相关矩阵怎么算的”会答不圆。assert不是摆设，数据行数或标签映射出错时能第一时间拦住，比后面 GMM 跑出诡异结果再回头查要省事得多。

3. 算法链路设计：PCA、LDA、GMM 与刀切法怎么串

3.1 马氏距离：把协方差结构带进分类判据

马氏距离是整条链路的粘合剂。欧氏距离假设各方向方差相同且特征独立，鸢尾花显然不满足：花瓣长度的方差远大于花萼宽度，特征之间还相关。马氏距离的定义是 $d(x,\mu)=\sqrt{(x-\mu)^T\Sigma^{-1}(x-\mu)}$，当 $\Sigma=I$ 时退化为欧氏距离，本质上是把数据白化后再算欧氏距离，所以它无量纲、不受量纲变化影响。

在方案里，马氏距离承担两个职责。一是 GMM 内部计算样本到各高斯成分均值的“归属强度”，虽然 GMM 的 E 步用的是带协方差的行列式和指数项，但判别阶段用马氏距离更直观；二是最终分类规则，样本归到马氏距离最小的那个成分对应的类别。这里有个容易忽略的点：协方差矩阵要加正则项，否则某类样本在降维后维度接近样本数时，$\Sigma$ 可能奇异，求逆直接报错。通常的做法是 $\Sigma+\lambda I$，$\lambda$ 取 $10^{-6}$ 到 $10^{-3}$ 量级。

3.2 高斯混合模型：极大似然估计与 EM 迭代

GMM 把总体分布写成 $p(x)=\sum_{k=1}^{K}\pi_k \mathcal{N}(x|\mu_k,\Sigma_k)$，$K$ 是成分数，$\pi_k$ 是混合权重。对数似然 $\sum_{i}\log p(x_i)$ 没有闭式解，用 EM 迭代：E 步算后验责任度 $\gamma_{ik}$，M 步更新 $\pi_k$、$\mu_k$、$\Sigma_k$。鸢尾花三类，$K$ 取 3 是自然选择，但作业里值得加一组 $K=2$ 和 $K=4$ 的对照，说明 $K$ 不是拍出来的。

from scipy.stats import multivariate_normal def gmm_em(X, K=3, max_iter=200, tol=1e-6, reg=1e-6): n, d = X.shape # 用分位数初始化均值，避免所有成分挤在一起 idx = np.linspace(0, n - 1, K).astype(int) mu = X[np.argsort(X[:, 0])][idx].copy() sigma = [np.cov(X.T) + reg * np.eye(d) for _ in range(K)] pi = np.ones(K) / K prev_ll = -np.inf for _ in range(max_iter): # E 步：责任度矩阵，行和为 1 resp = np.zeros((n, K)) for k in range(K): resp[:, k] = pi[k] * multivariate_normal.pdf(X, mu[k], sigma[k]) ll = np.log(resp.sum(axis=1) + 1e-12).sum() resp /= resp.sum(axis=1, keepdims=True) # M 步：加权更新权重、均值、协方差 nk = resp.sum(axis=0) pi = nk / n mu = (resp.T @ X) / nk[:, None] for k in range(K): diff = X - mu[k] sigma[k] = (resp[:, k][:, None] * diff).T @ diff / nk[k] + reg * np.eye(d) if abs(ll - prev_ll) < tol: break prev_ll = ll return pi, mu, sigma, resp

参数说明：reg是协方差正则项，防止奇异；tol控制对数似然收敛阈值，鸢尾花这种规模 200 次迭代足够；初始化用排序后等间隔取点，比纯随机稳，能减少 EM 陷入局部最优的概率。这段代码没有调库的GaussianMixture，是为了让作业里“极大似然估计”这句话有对应的实现，答辩时能指着 E 步 M 步讲。

3.3 PCA 与 LDA：无监督降维和有监督降维的对照

PCA 的目标是找线性变换 $u$，使投影方差最大，解是协方差矩阵（或相关矩阵）最大特征值对应的特征向量。鸢尾花四维降到二维，通常前两个主成分能解释 95% 以上的方差。LDA 则利用标签，目标是最大化类间散度与类内散度的比值，解是 $S_w^{-1}S_b$ 的最大特征值对应特征向量。两者放在一起对比，正好说明“有标签时 LDA 通常优于 PCA”这个结论，而不是空口断言。

def pca_fit(X, n_components=2): R = np.corrcoef(X, rowvar=False) # 相关矩阵消除量纲 eigvals, eigvecs = np.linalg.eigh(R) order = np.argsort(eigvals)[::-1] eigvals, eigvecs = eigvals[order], eigvecs[:, order] return eigvecs[:, :n_components], eigvals def lda_fit(X, y, n_components=2): classes = np.unique(y) mu_all = X.mean(axis=0) Sw = np.zeros((X.shape[1], X.shape[1])) Sb = np.zeros_like(Sw) for c in classes: Xc = X[y == c] muc = Xc.mean(axis=0) Sw += (Xc - muc).T @ (Xc - muc) Sb += len(Xc) * np.outer(muc - mu_all, muc - mu_all) eigvals, eigvecs = np.linalg.eigh(np.linalg.inv(Sw + 1e-6 * np.eye(Sw.shape[0])) @ Sb) order = np.argsort(eigvals)[::-1] return eigvecs[:, order[:n_components]], eigvals[order]

pca_fit用np.corrcoef而不是np.cov，对应前面说的“用相关矩阵代替协方差矩阵消除量纲”。lda_fit里Sw加正则项同理，鸢尾花每类 50 个样本、四维特征，Sw本来可逆，但降维到一维再算时数值稳定性会变差，加正则是个保险。LDA 最多只能降到 $C-1=2$ 维，这一点在作业里要写清楚，否则有人试图降到三维会得到无意义的第三个方向。

3.4 刀切法：用留一法给降维方案打分

刀切法（留一法）的流程是：每次从训练集剔除一个样本，用剩余样本建立判别规则，再预测被剔除的样本，重复 150 次，统计判错个数，判错率 $= \sum k_i / n$。它比简单回报法更客观，因为每个样本都当过测试集；比 k 折交叉验证更彻底，因为训练集几乎用满。代价是运算量大，150 次 GMM 拟合在普通笔记本上要跑几分钟，但完全可接受。

方案里刀切法要比较的对象是“降维方法 × 降维维度”的组合：PCA 降到 1 维、2 维，LDA 降到 1 维、2 维，再加一组不降维的原始四维作为基线。每组跑一遍留一法，记录判错率，最后画一张对比表。这张表就是整份作业的结论支撑，比任何文字描述都有说服力。

4. 实现与结果分析：从类设计到对比表

4.1 模块划分与类接口

按作业里的 UML 思路，PCA 和 LDA 各设计成一个类，fit接收训练样本求变换矩阵，transform做投影，fit_transform串起来。GMM 单独一个类，fit跑 EM，predict用马氏距离判类。刀切法写成函数，接收“降维器 + 分类器”的组合，返回判错率。这样拆的好处是刀切法里换降维方法只改一个参数，不用复制粘贴五份代码。

class PCAReducer: def __init__(self, n_components=2): self.n_components = n_components def fit(self, X, y=None): self.W_, self.eigvals_ = pca_fit(X, self.n_components) return self def transform(self, X): return X @ self.W_ class LDAReducer: def __init__(self, n_components=2): self.n_components = n_components def fit(self, X, y): self.W_, self.eigvals_ = lda_fit(X, y, self.n_components) return self def transform(self, X): return X @ self.W_

接口统一后，刀切法的主循环可以写成通用形式：对每个样本，用剩余样本fit降维器和 GMM，再transform被留出的样本、算马氏距离、取最小距离对应的类别。注意 LDA 的fit需要标签，PCA 不需要，所以通用循环里要把y传进去，PCA 类用y=None兼容。

4.2 降维维度与判错率的对比结果

跑完五组配置，典型结果会呈现这样的规律：原始四维基线判错率在 2% 到 4% 之间；PCA 降到二维判错率接近基线，降到一维明显变差，因为第一主成分主要承载花瓣长度方向的信息，Setosa 和另外两类能分开，但 Versicolour 和 Virginica 有重叠；LDA 降到二维判错率最低，通常能到 2% 以下，降到一维也比 PCA 一维好，因为它用了标签，投影方向直接冲着类间可分性去。

这张表要配合解释：PCA 二维接近基线说明四维里确实有冗余，两个主成分保住了大部分判别信息；PCA 一维掉得厉害说明“方差最大”不等于“判别性最强”；LDA 二维最优印证了有监督降维在分类任务上的优势。数字不必和上表完全一致，但趋势要能自圆其说，答辩时被追问“为什么你的 PCA 一维这么差”要能答到“第一主成分被花瓣长度主导，类间重叠没被拉开”。

4.3 结果解读里容易漏掉的两个细节

第一个细节是标准化对 PCA 的影响。如果不标准化直接用协方差矩阵，花瓣长度的方差远大于其他特征，第一主成分几乎就是花瓣长度本身，降维效果和标准化后完全不同。作业里要明确写“采用相关矩阵”，并在附录放一张标准化前后的特征值对比，说明为什么这么选。

第二个细节是 GMM 的类别分配。GMM 是无监督的，跑出来的成分编号和真实标签没有对应关系，判类前要先做一次匹配，比如用混淆矩阵找最优排列，或者直接用马氏距离到各类均值的最近邻规则绕开编号问题。我一般选后者，逻辑简单，答辩时也好解释：GMM 负责估计分布参数，判别用马氏距离，两步分开说。

5. 避坑与排查：这类大作业最容易翻车的五个点

5.1 协方差矩阵奇异导致求逆报错

现象是 GMM 迭代到一半抛LinAlgError: Singular matrix，或者马氏距离算出inf。原因通常是某类样本在降维后维度接近或超过样本数，或者某个高斯成分只分到极少数样本，协方差估计退化。解决办法是给协方差加正则项 $\lambda I$，$\lambda$ 从 $10^{-6}$ 起调；同时在 EM 里加一个最小样本数保护，某成分责任度之和小于阈值就重新初始化，而不是硬算。

5.2 标准化口径不一致导致 PCA 和 LDA 结果对不上

现象是同一份数据，PCA 用相关矩阵、LDA 用原始协方差，两者降维后的尺度差一个量级，对比表看起来像两组不同数据。原因是 PCA 前做了 z-score，LDA 前没做。解决方法是整条链路统一：先标准化一次，PCA 和 LDA 都在标准化后的数据上跑，LDA 的Sw、Sb自然也在同一尺度下。这一点在代码里体现为standardize只调用一次，结果传给两个降维器。

5.3 刀切法循环里重复拟合造成结果不可复现

现象是同一份代码跑两次判错数不一样。原因是 GMM 初始化带随机性，而刀切法每次循环都重新fit，随机种子没固定。解决办法是在刀切法入口设np.random.seed(42)，并且 GMM 初始化改成确定性的分位数取点（如 3.2 节代码所示）。如果坚持随机初始化，就在每次fit前重置种子，保证 150 次循环的随机序列可复现。

5.4 LDA 降维维度超过 C-1 得到无意义方向

现象是试图把鸢尾花 LDA 降到三维，第三个特征值接近零，投影后数据几乎塌缩。原因是 LDA 的类间散度矩阵秩最多为 $C-1$，三类数据最多两个有效判别方向。解决办法是在lda_fit里加断言n_components <= len(classes) - 1，作业正文里也写明这个约束，避免读者照抄时踩同一个坑。

5.5 判错率统计把训练样本和测试样本混在一起

现象是判错率异常低，接近 0。原因是判别规则用了被预测样本本身的信息，比如先对所有样本拟合 GMM 再留一预测，等于测试集泄漏进训练。解决办法是严格按刀切法流程：每次剔除一个样本，用剩余 149 个样本重新fit降维器和 GMM，再预测被剔除的那个。运算量是大，但这是留一法的定义，不能偷懒。

6. 验收与答辩：怎么证明这套分析链站得住

验收这份大作业，我一般看四样东西，按重要性排序。第一是流程完整性：数据读取、标准化、PCA、LDA、GMM、马氏距离判别、刀切法比对，七个环节缺一不可，缺了哪个，答辩时被问到就是硬伤。第二是结果可复现：代码里固定随机种子，附录给出运行环境和依赖版本，别人拿到能跑出同样的判错率表。第三是对比有解释：PCA 一维为什么差、LDA 二维为什么好，要落到特征相关性和标签信息上，不能只说“实验结果符合预期”。第四是符号和公式自洽：第二章的符号表、第三章的公式、第四章的代码变量名要对得上，这是数理统计作业区别于纯编程作业的地方。

答辩时最常被追问的三个问题，提前准备好答案。一是“为什么用相关矩阵而不是协方差矩阵做 PCA”，答量纲不一致，花瓣长度方差远大于花萼宽度，不标准化会让第一主成分被单一特征主导。二是“GMM 是无监督的，你怎么保证成分和类别对应”，答判别阶段用马氏距离到各类均值的最近邻规则，GMM 只负责估计分布参数，两步解耦。三是“刀切法运算量大，有没有考虑 k 折”，答留一法是刀切法的标准形式，150 个样本规模可接受，且训练集利用率最高，适合小样本场景。

一个具体技巧：把五组配置的判错率表放在结论页最显眼的位置，旁边配一句“LDA 二维在留一法下判错率最低，验证了有监督降维在分类任务上的优势”。评委翻到这一页，三秒内能抓住你的核心结论，后面的追问也会围绕这张表展开，你准备起来更有针对性。

从那以后我每次带这类分析型作业，都强制要求先跑通“原始基线 + 两种降维 × 两个维度”这五组留一法，再动笔写正文。因为只有这张表出来了，你才知道哪些结论能写、哪些话是拍脑袋。希望帮到你。

本文还有配套的精品资源，点击获取