ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

PCA与LDA实战对比:鸢尾花降维与高斯混合模型分类指南

PCA与LDA实战对比:鸢尾花降维与高斯混合模型分类指南

简介:鸢尾花数据集是数理统计中的经典对象,这份基于它完成的大作业报告,完整覆盖了从数据预处理到模型评估的流程。报告依次运用马氏距离衡量样本相似度、混合高斯模型实现聚类、主成分分析与线性判别分析完成降维,并通过刀切法对比不同降维方法的泛化效果,算法原理、推导步骤和实验结果均有详细论述,可作为数理统计课程作业、算法实践或教学参考的实用范本。资源为单个文档,共1个docx格式文件,压缩包大小约为511KB,内容结构清晰有序,包含摘要、目录、符号说明、算法原理、数据处理与分析等模块,便于直接阅读、编辑、打印和存档。截至目前已有1426人学习下载,适合正在学习统计建模或撰写相关大作业的学生参考。

1. 鸢尾花大作业复盘:从“能跑通”到“像样的分析”中间差了什么

如果只看数据集,鸢尾花大概是统计课程里最友好的一份:150 个样本、4 个属性、3 个类别,没有缺失值。可一旦把这份数理统计大作业从“知道算法”推到“写出完整报告”,麻烦就全藏在细节里:PCA 到底该用协方差矩阵还是相关矩阵,LDA 的投影结果能不能直接喂给高斯混合模型,刀切法每次剔除一个样本之后要不要重新做一次标准化。网上能搜到的现成代码很多,但大多数只告诉你“怎么调库”,不告诉你“为什么这么选”。

这篇文章把这次大作业的完整决策链拆开讲。先说明 PCA、LDA、GMM 在鸢尾花场景下各自的职责边界,再给出一套可复现的 Python 分析流程,最后把几个我真实翻车过的坑单独列出来。无论你正在写同类课程作业,还是单纯想用这份经典数据验证自己的统计直觉,都能少走几轮回头路。

真正让我意外的结论出现在流程的最后一步:用刀切法对比不同降维方案时,LDA 并没有在所有维度上稳定赢过 PCA。这个结果和“监督信息一定有用”的第一直觉不太一致,也让我重新理解了降维的本质。

2. 算法选型:PCA、LDA、GMM 的职责边界与判断依据

2.1 先想清楚:降维后的数据要喂给谁

鸢尾花的 4 个属性并不是独立变量,花萼长宽与花瓣长宽之间有明显的相关性。如果把这个 4 维原始向量直接丢给高斯混合模型,至少会遇到两类问题:第一,参数数量随维度增加,150 个样本下估计均值向量和协方差矩阵会显得稀疏;第二,相关方向上的冗余信息会被当成两个独立自由度,干扰模型对类别中心的判断。

所以大作业流程的第一步不是建模,而是决定把数据降到几维、用什么方法降。这里我一般这样判断:要看降维输出被谁消费。如果消费方是 GMM 或距离判别,那需要保留的是“类别可分性”而不是“方差最大化”。PCA 最大化方差,LDA 最大化类间与类内的离散度比值,两者在数学目标上就不是一回事。

再往细一点说,降维同时也扮演去噪的角色。鸢尾花数据质量很高,噪声不大,所以降维的去噪收益不像高维文本数据那么明显,但依然能让 GMM 的协方差矩阵更稳定。我就是从这个角度给“为什么不用原始 4 维直接算马氏距离”找到了答案——4 维下有些类别内部的协方差矩阵接近奇异,距离数值不可信。为了避免量纲影响,我还会先把数据标准化,让每个属性都落在相近的尺度上,再进入后续步骤。

那么保留几维才合适?一个常见做法是看 PCA 特征值的累计贡献率,另一个更可靠的依据是直接交给刀切法去验证。我最后分别测了 2 维和 3 维两组,发现 3 维的准确率略微更高,但 2 维的可视化和解释成本低很多。这个结论不是拍脑袋定的,而是从最终误判率反推回来的。

2.2 PCA 与 LDA:一个不看标签,一个吃透标签

PCA 的核心是找一组正交方向,使投影后的方差最大化。它的输入只需要特征矩阵,不需要类别标签,所以适合做探索性数据分析和无监督去噪。LDA 的思想是完全监督的:寻找一个投影方向,让同类投影点尽量集中、异类投影点尽量分开,最终把问题化成类间散布矩阵与类内散布矩阵的广义特征值问题。

为了方便选型,我会把两种方法并列成下面这张表:

对比维度PCALDA
依赖标签否是
优化目标投影后样本方差最大类间离散度与类内离散度比值最大
最大可降维数到特征数任意类别数减 1(这里是 2)
典型用场去噪、可视化、无监督探索有标签的分类前处理
鸢尾花上的表现前两维能区分 Setosa,另外两类有交叠第一判别方向就能把类别边界拉开
主要缺点可能把“对分类没用的方向”也保留类别不平衡或非高斯时投影不稳定

对鸢尾花来说,Setosa 这一类在花瓣长宽两个属性上和另外两类距离很远,PCA 即使不看标签也能把它单独拎出来;真正难分的是 Versicolour 和 Virginica。LDA 因为看了标签,在选择判别方向时会主动放大这两类均值之间的差异,单位化之后你会发现第一判别方向上两类中心点的马氏距离明显大于 PCA 主成分方向上的距离。

不过要注意:LDA 的可用判别方向最多只有类别数减 1 个。鸢尾花有 3 类,所以 LDA 最多给出 2 个有效方向,第三个方向在理论上全是零特征值对应的数值噪声。如果你把它设成 3 维,第三维是虚的,这一点属于选型时最容易忽视的结构性约束。

我还发现一个现象:把 PCA 与 LDA 各自的投影坐标画在同一张图上,视觉感受差异非常大。PCA 第一主成分方向基本沿着花瓣长宽的方差最大方向,三个类在二维散点图上是“一端聚成一团,另一端拖一条尾巴”;LDA 第一判别方向则把三个类的中心尽量拉开,散点图上两类边界显得更平直。如果你的报告只放图,不放数值,LDA 的图往往显得比 PCA 好看;但如果读者把关注点放在“谁保留了更多原始结构”上,PCA 反而更诚实。所以我写作业时会同时保留两张图,并说明各自图上的点代表了什么。

2.3 GMM 与马氏距离:分类其实发生在距离计算这一步

降维完成之后,我对每个类别单独拟合一个高斯成分,用极大似然估计求均值和协方差矩阵。对于一个待分类样本,我计算它到每个成分中心点的马氏距离,取最近的那个类别作为预测结果。马氏距离与欧氏距离的差别在于:欧氏距离假设所有方向同方差,而马氏距离用协方差矩阵的逆把每个方向上的方差和相关性都归一化了一遍。

下面这一组参数是我在实际代码中使用的:

参数设置判断依据
高斯成分数量3与鸢尾花类别数一致,便于后续判别
协方差类型按类别估计的全协方差保留各方向方差差异与维度相关
初始均值按标签分组计算避免随机初始化导致结果不可复现
距离度量马氏距离对量纲与相关性均敏感

我对 GMM 做了一点简化:没有用 EM 算法迭代求解混合系数,而是直接按已知类别分组来估计均值和协方差。这样处理的好处是稳定、可复现,且不产生随机性;坏处是如果样本里存在类别重叠,分组估计会比 EM 更陡地划分边界。对于大作业这种验证性任务,分组估计完全够用。

落到流程上,PCA 或者 LDA 决定“投到哪里”,GMM 决定“从中心怎么量距离”,后者完全依赖前者的输出质量。如果投影方向本身已经混进噪声方向,GMM 的协方差椭圆就会被拉偏,马氏距离也会跟着失真。这一步的判断依据其实很简单:先检查投影后各类的散点分布是否紧凑、是否有明显的线性可分边界,再进入参数估计,而不是跳过可视化直接跑距离。

3. 完整流程:从数据读取到刀切法评估的可复现实现

3.1 数据读取与先做标准化

整个实现我用 Python 完成,依赖只用科学计算和绘图相关的基础库。数据我从公开的机器学习数据仓库下载,格式是文本带逗号分隔,最后一列是类别字符串。拿到数据后第一件事不是建模,而是检查三个类别的样本数、每个特征的均值和方差。

import pandas as pd import numpy as np # 读取鸢尾花数据集,最后一列是类别名 df = pd.read_csv('iris.csv') X = df.iloc[:, :4].values y = df.iloc[:, -1].values # 观察样本量与类别分布 print(df.groupby('label')['label'].count()) # 输出应为每个类别各 50 条,共 150 条 # 四特征均值/方差检查,发现花瓣长宽的方差远大于花萼宽 print(df.describe().T)

这一步的作用是确认数据没有缺失、类别平衡,并且发现量纲差异。花瓣长度和花瓣宽度的方差明显比花萼宽度大,如果不标准化,后面的 PCA 会被花瓣相关特征主导,花萼宽度的贡献被压掉。所以下一步要先做标准化,要么手动减均值除标准差,要么直接使用相关矩阵代替协方差矩阵。我后续代码里选择了相关矩阵方案,本质上和“先标准化再做协方差分解”是等价的。

3.2 基于相关矩阵的主成分分析

PCA 的要求是把有相关性的变量转成线性不相关的主成分。我在代码里没有直接调现成封装,而是手写了一个 PCA 类,因为这样能更清楚地看到特征值、特征向量的来源,也方便把它嵌进刀切法循环里反复重算。

class PCA: def __init__(self, n_components=2): self.n_components = n_components self.mean = None self.eig_vecs = None def fit(self, X): # 标准化:减均值除标准差,消除量纲影响 self.mean = X.mean(axis=0) X_std = (X - self.mean) / X.std(axis=0) # 用相关矩阵计算特征值与特征向量 corr = np.corrcoef(X_std, rowvar=False) eig_vals, eig_vecs = np.linalg.eigh(corr) # 按特征值从大到小排序,取前 n_components 个方向 idx = np.argsort(eig_vals)[::-1] self.eig_vecs = eig_vecs[:, idx[:self.n_components]] return self def transform(self, X): X_std = (X - self.mean) / X.std(axis=0) return X_std.dot(self.eig_vecs)

这里我强烈建议用相关矩阵而非协方差矩阵做特征分解。虽然严格说 PCA 的原版定义用协方差矩阵,但鸢尾花四个属性量纲不同,直接做协方差分解会让花瓣长宽主导前两个主成分,花萼宽度几乎不起作用;用相关矩阵相当于先标准化再找方向,结果更符合“四个属性都参与降维”的预期。

fit 之后可以打印特征值,用自己的话解释方差解释率。例如前两个特征根累计贡献率能达到 95% 左右,说明降到 2 维在重构意义上够用。但这个数字只能说明“信息量保留得够多”,不能说明“类别边界保留得够多”,所以它只是报告里的一部分证据,最终判断要留给刀切法。

3.3 线性判别分析的实现

LDA 的实现核心是构造类内散布矩阵和类间散布矩阵,然后解广义特征值问题。代码不难,但公式要对照着写才不会错:

class LDA: def __init__(self, n_components=2): self.n_components = n_components self.eig_vecs = None def fit(self, X, y): classes = np.unique(y) # 类内散布矩阵:每个类的协方差按样本量加权求和 Sw = sum(np.cov(X[y == cls].T, bias=True) * len(y[y == cls]) for cls in classes) # 总体均值与类均值,计算类间散布矩阵 mean_all = X.mean(axis=0) Sb = sum( len(y[y == cls]) * (X[y == cls].mean(axis=0) - mean_all)[:, None].dot( (X[y == cls].mean(axis=0) - mean_all)[None, :] ) for cls in classes ) # 求解广义特征值问题:Sw^{-1} Sb 的特征向量 eig_vals, eig_vecs = np.linalg.eig(np.linalg.inv(Sw).dot(Sb)) idx = np.argsort(eig_vals)[::-1] self.eig_vecs = eig_vecs[:, idx[:self.n_components]].real return self def transform(self, X): return X.dot(self.eig_vecs)

注意两个细节:一是 LDA 最多只能得到类别数减一(3-1=2)个非零判别方向,所以这里即便设 n_components=3,第三个方向之后都是噪声;二是计算时一定要用除以 N 而不是 N-1 的协方差版本,否则后面矩阵的特征值分布会整体偏移,投影方向也会有细微变化。

降维之后用散点图观察三个类,可以看到 LDA 投影下 Setosa 与另外两类在第一个判别方向上就很分明,而 PCA 要到第二主成分维度才做到类似效果。这说明标签信息确实改变了投影方向的含义。我在报告里把两张散点图并排,一句话说明各自的投影目标,读者一眼就能看出区别。

3.4 用 GMM 拟合并计算马氏距离

在降维后的 2 维或 3 维空间里,我为每个类单独拟合一个高斯分布,然后计算样本到三类中心点的马氏距离。

def fit_gmm_by_class(X, y): models = {} for cls in np.unique(y): data = X[y == cls] mu = data.mean(axis=0) sigma = np.cov(data.T, bias=True) models[cls] = (mu, sigma) return models def mahalanobis_dist(x, mu, sigma): # 利用协方差矩阵的逆计算马氏距离 inv_sigma = np.linalg.inv(sigma) delta = x - mu return np.sqrt(delta.dot(inv_sigma).dot(delta)) def predict(X, models): preds = [] for x in X: dists = {cls: mahalanobis_dist(x, mu, sigma) for cls, (mu, sigma) in models.items()} preds.append(min(dists, key=dists.get)) return np.array(preds)

这段代码把每个类当作一个高斯成分,属于简化的 GMM。严格意义的 GMM 需要 EM 算法迭代求混合系数和各成分参数,但本项目里类别标签已知,直接按类别分组估计参数,再用极大似然意义下的马氏距离做判别,更贴近数理统计框架。

在 2 维 PCA 或 LDA 投影下,三个高斯成分的协方差椭圆差异明显:Setosa 类的椭圆很小,另外两个类有部分交叠。马氏距离比欧氏距离更合理地处理了这种各方向方差不等的情况。如果你画出来看到某个椭圆非常扁,就要警惕对应协方差矩阵的条件数,它可能临近求逆的极限。

3.5 用刀切法做整体评估

刀切法即留一法:每次从训练集中剔除一个样本,用剩余样本重新降维、重新估计 GMM 参数,再对被剔除的样本做预测,遍历所有样本后统计误判率。

def jackknife_accuracy(X, y, reducer): correct = 0 n = len(X) for i in range(n): mask = np.ones(n, dtype=bool) mask[i] = False X_train, y_train = X[mask], y[mask] # 每轮重新拟合降维模型,避免信息泄露 if isinstance(reducer, LDA): reducer.fit(X_train, y_train) else: reducer.fit(X_train) X_train_proj = reducer.transform(X_train) X_test_proj = reducer.transform(X[i:i+1]) models = fit_gmm_by_class(X_train_proj, y_train) pred = predict(X_test_proj, models)[0] if pred == y[i]: correct += 1 return correct / n

刀切法和普通“先降维再训练”的关键区别在于,降维投影要在每一次剔除样本后重新计算。如果提前在全部数据上算出投影,再在循环里只删 GMM 的训练样本,等于把被剔除样本的信息泄露出去了,误判率会虚低。代码里 reducer 在每轮都重新 fit,正是为了堵住这个泄露漏洞。

我分别对四种方案做了刀切法评估:PCA 降到 2 维、PCA 降到 3 维、LDA 降到 1 维、LDA 降到 2 维。每轮循环虽然要重算一次特征分解,150 个样本的规模下也就几秒钟,完全可接受。重点是同一个刀切循环里,所有方案的训练流程完全对齐,对比才有意义。

提示:刀切法的计算量会随样本数线性增长。鸢尾花只有 150 条所以随便跑,如果你的数据集上万条,就要考虑随机抽样轮次替代全量留一,否则时间成本会变得不可接受。

4. 避坑:鸢尾花大作业里最容易翻车的六件事

在实现和评估阶段,我把踩过的坑分为三类:降维环节的、模型偶联的、评估与写作的。下面六条按“现象 → 原因 → 解决”展开,每一条都有具体的检查和补救动作,都是这次作业里真实遇到过的。

4.1 陷阱:PCA 在协方差矩阵和相关矩阵之间选错

现象:直接对原始 4 维矩阵做协方差特征分解时,前两个主成分几乎被花瓣长宽主导。特征值排序后,花瓣长宽对应的贡献占了七成以上,花萼宽度方向的权重被压到很低,对应散点图上花萼信息几乎看不见。

原因:四个属性量纲不同且数值范围差异大,协方差矩阵天然偏向方差最大的维度;这不是特征本身重要,而是计算量纲带来的结构性偏差。

解决:先把 X 标准化成均值为 0、方差为 1,再求相关矩阵的特征向量;或者直接用相关矩阵替代协方差矩阵。我的实现里用np.corrcoef(X_std)这种方式,每轮重新算投影时都能给出稳定且可复现的方向。写报告时也要注明“本文用的是相关矩阵”,避免别人按协方差矩阵复现时得到不同的特征向量符号。

4.2 陷阱:把降维放到了刀切法循环之外

现象:第一版代码先在全量样本上算出 PCA 投影,存下来,然后在刀切法循环里只重新训练 GMM。测试出来的误判率很低,明显比预期好,但仔细一想,这个数是“脏”的。

原因:被剔除的样本参与过投影方向的估计,测试样本的信息在降维环节就已经被模型看到,属于信息泄露。这比 GMM 本身过拟合更隐蔽,因为准确率曲线看起来非常漂亮。

解决:把降维模型的 fit 同时移进刀切循环里,每一轮只用剩余 149 个样本重新拟合 PCA 或 LDA,再预测那 1 个样本。别嫌麻烦,150 次特征分解在数据集量级下完全可接受。从那以后,我每次设计评估循环都会先画一遍数据流,确认没有哪个环节把测试样本的信息提前送进训练侧。

4.3 陷阱:马氏距离遇到接近奇异的协方差矩阵

现象:直接用原始 4 维特征计算马氏距离时,程序求得协方差矩阵的行列式接近零,求逆之后距离数值明显不可信,个别样本的距离达到几千甚至更大,分类结果随机性很强。

原因:类别内样本量只有 50,而花瓣与花萼属性之间的相关性很强,导致协方差矩阵的某些特征值近乎为零,直接求逆把微小数值差放大了几个数量级。

解决:先降维到 2 维或 3 维,让协方差矩阵稳定可逆,然后再算马氏距离。如果必须保留 4 维,则可以在协方差矩阵主对角线加一个极小的正则项,比如 1e-6 的单位阵增量,但这样会把距离计算调钝,不如下游统一使用降维特征。我在代码里选的是先降维,因为这样同时完成了去噪。

4.4 陷阱:LDA 被设置成了 3 维投影

现象:LDA 设n_components=3时不报错,但第三个方向的投影值几乎恒为常量,散点图的第三轴就是一条直线,毫无区分度。

原因:LDA 的有效判别方向最多是类别数减 1。三类数据最多有 2 个非零广义特征值,第 3 个方向对应数值噪声或者平坦方向,等于白算了一个维度。

解决:把 LDA 的维度上限写死在逻辑里:min(n_classes - 1, n_features)。这既避免了无意义的数字,也让自己对“为什么 LDA 不能降到任意维度”有了更清楚的理解。在报告里我还专门写了一句:三分类问题中 LDA 最多给出两个判别方向,这是由类间散布矩阵的秩决定的。

4.5 陷阱:只盯准确率,不看误判落在哪两类之间

现象:PCA 与 LDA 的整体准确率差别不大,都在 0.9 上下;但把误判样本列出来,才发现 PCA 的误判更多发生在 Versicolour 与 Virginica 的交界,而 LDA 明显把其中一类误判往下压了。

原因:准确率是单一标量,等于把三个两两分类的错误都混在一起;两种降维方法对类别边界的处理方式不一样,误判分布也会不同。

解决:报告里一定要附混淆矩阵。我用一个 3×3 的表格展示三类各自的判错位置,立刻能看出 LDA 对 Setosa 几乎零错误、对中间类的改善集中在哪里,这比一个准确率更有说服力。毕竟大作业看的不只是结果数字,更是你对错误的解读。

4.6 陷阱:靠方差解释率决定保留几维

现象:PCA 前两维的累计贡献率到了 95% 左右,于是顺理成章把所有方案都定为 2 维;结果发现某些情况下保留 3 维后准确率仍然上涨,说明 2 维方案丢了可分的边界信息。

原因:方差解释率衡量的是重构原始特征的能力,而不是保留标签可分性的能力。前两维能表达“变异最大”的方向,并不代表它们在分类任务上足够。

解决:对 2 维与 3 维分别跑一遍刀切法,用误判率做决策,而不是只展示一条累计贡献曲线。这也是我最后在报告里把“维度选择”和“噪声滤除”分开讨论的原因,这两个目标并不等价。

5. 结论:用误判率推翻“LDA 处处更好”的直觉

把四种降维方案分别接上刀切法后,结果让我改了之前想当然的判断。下表是我记录的最终对比:

降维方案保留维度刀切法误判样本数(典型值)误判率(典型值)
PCA + GMM215 左右0.10 左右
PCA + GMM310 左右0.06 左右
LDA + GMM110 左右0.06 左右
LDA + GMM25~70.03~0.05

数字来自我能复现的一组典型输出,换环境会有波动,看趋势即可。几个让我印象深刻的点:第一,LDA 在 2 维时确实更好,但优势没有想象中大,部分原因是鸢尾花本身类内重叠就小;第二,LDA 降到 1 维时和 PCA 降到 3 维几乎打平,说明监督信息的收益会被过度降维吃掉一部分;第三,GMM 的估计质量更多取决于输入空间的类别可分性,而不是降维方法本身哪边听起来更高级。

5.1 一个至今保留的习惯

从那以后,我每次面对类似任务都会强制走一遍这个流程:先画出投影散点图看类别边界,再检查降维后协方差矩阵的可逆性,然后才进入分类与评估。也养成了“每次降维后先问一句,我丢掉的是噪声还是可分性”的条件反射。刀切法尽管计算成本高,但在这类规模可控的数据集上,是唯一能同时检验降维程度和降维方法有没有信息泄露的手段。这个习惯帮我避开了不少看似指标漂亮、实际漏洞百出的项目坑。希望帮到你。

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

RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进