高斯判别分析GDA：从数学推导到Python实战

1. 项目概述GDA是做什么的为什么值得学1.1 从分类问题说起生成模型和判别模型的本质区别各位做机器学习的朋友如果你已经学过逻辑回归那我问你一个问题逻辑回归到底在做什么表面上看它是在拟合一个决策边界把正负样本分开。但更本质地说它是在直接建模 P(y|x)也就是给定特征 x直接估计类别 y 的概率。这类方法统称为判别模型它们不关心数据是怎么生成的只管找一条边界。而高斯判别分析Gaussian Discriminant Analysis简称 GDA走的是完全不同的路子。它是一个典型的生成模型核心思想是先假设每一类数据分别服从一个高斯分布然后反过来用贝叶斯公式算出 P(y|x)。打个比方判别模型就像你要判断一个人是南方人还是北方人你直接看他的口音特征来分类生成模型则是你先建立南方人的身高、体重、口音分布和北方人的身高、体重、口音分布两个完整档案来了新人对照两份档案看哪边更匹配。GDA 这个名字听起来很数学但它其实就是假设数据服从高斯分布 用最大似然估计学参数 贝叶斯决策这三个标准步骤的组合拳。它在小样本场景下尤其有优势因为高斯分布的假设本身就带有很强的结构性先验数据量少的时候也能学得比较稳。而且它是理解 LDA线性判别分析、QDA二次判别分析、朴素贝叶斯等一系列模型的金钥匙。1.2 这篇文章适合谁能解决什么问题这篇文章适合以下几类读者正在学吴恩达机器学习课程被 GDA 推导卡住的学生做分类任务时只知道调 sklearn 包但对模型原理一知半解的工程师以及准备面试、想深入理解生成模型与判别模型差别的求职者。我会从数学原理讲到代码实现再补充实战中的坑和调参经验。读完这篇文章你至少能回答以下问题GDA 的参数是怎么估计出来的为什么共享协方差矩阵时决策边界是线性的GDA 和逻辑回归到底该选哪个协方差矩阵奇异了怎么办这些都是在实际项目中真正会碰到的困惑我尽量用大白话和代码给你讲透。2. 数学原理与模型推导知其然更知其所以然2.1 多元高斯分布理解协方差矩阵比想象中重要要理解 GDA首先得熟悉多元高斯分布。它的概率密度函数长这样p(x; μ, Σ) 1 / ((2π)^(n/2) |Σ|^(1/2)) * exp(-1/2 * (x-μ)^T Σ^(-1) (x-μ))这个公式看起来吓人但拆开理解很简单。μ 是均值向量决定了分布的中心位置Σ 是协方差矩阵决定了分布的形状和方向。你可以把二维高斯分布想象成一个椭圆Σ 的特征向量决定了椭圆的长轴短轴方向特征值决定了椭圆在各方向上的伸展程度。在 GDA 里我们通常会做一个关键假设所有类别共享同一个协方差矩阵 Σ。这个假设是 LDA 和 QDA 的分水岭。如果共享 Σ决策边界是线性的如果每个类别有自己的 Σ_k决策边界是二次的QDA。为什么后面我会给出几何直觉和数学推导。实操中很多人忽略了一件事检查特征是否近似高斯分布。GDA 的G不是白叫的它默认数据本身是从高斯分布里采样出来的。如果特征严重偏态或者存在明显长尾直接套 GDA 效果会很难看。常见的补救办法是取对数、开根号、Box-Cox 变换等先把特征掰成近似高斯再说。2.2 模型假设与参数估计最大似然估计的完整推导GDA 的完整模型假设如下y ~ Bernoulli(φ) # 类别标签服从伯努利分布 x | y0 ~ N(μ0, Σ) # 负类样本服从高斯分布 x | y1 ~ N(μ1, Σ) # 正类样本服从高斯分布注意两类共享同一个 Σ但均值不同μ0 和 μ1。我们的目标是根据训练数据估计四个参数φ、μ0、μ1、Σ。方法是最经典的最大似然估计MLE。对数似然函数写出来L(φ, μ0, μ1, Σ) Σ log p(x^i, y^i; φ, μ0, μ1, Σ) Σ log p(x^i | y^i; μ0, μ1, Σ) Σ log p(y^i; φ)接下来对每个参数求偏导并令其为 0。这里我不做冗长的矩阵微分推导直接给出直观理解φ 1/m * Σ 1{y^i 1} # 正类样本占比 μ0 (Σ 1{y^i 0} x^i) / (Σ 1{y^i 0}) # 负类样本的均值 μ1 (Σ 1{y^i 1} x^i) / (Σ 1{y^i 1}) # 正类样本的均值 Σ 1/m * Σ (x^i - μ_y^i)(x^i - μ_y^i)^T # 每个样本减去它所属类别的均值这几个公式的直觉非常清晰φ 就是训练集中正类的比例μ0 是负类样本的质心μ1 是正类样本的质心Σ 是所有样本围绕各自类别质心的平均散开程度因为共享 Σ所以要把两类数据放一起算但每个样本要减去它自己类别的均值。共享协方差矩阵的几何含义是假设两个类别的高斯分布形状和朝向完全相同只是中心位置不同。这在二维平面上就是两个大小、方向一样的椭圆只是圆心挪了个位置。如果现实里两个类别的数据散布明显不同比如正类数据很集中、负类数据很分散那共享 Σ 的假设就不合适应该改用 QDA。2.3 判别边界的生成为什么共享协方差时边界是线性的有了参数估计接下来是决策。根据贝叶斯公式p(y1 | x) p(x | y1) * p(y1) / p(x)我们选概率更大的类别作为预测结果。比较 p(y1|x) 和 p(y0|x)等价于比较p(x | y1) * p(y1) vs p(x | y0) * p(y0)两边取对数代入高斯分布的密度函数一堆项会抵消。当 Σ 共享时二次项 -1/2 x^T Σ^(-1) x 在两边是相同的互相抵消剩下的只有一次项(μ1 - μ0)^T Σ^(-1) x (1/2)(μ0^T Σ^(-1) μ0 - μ1^T Σ^(-1) μ1) log(φ / (1-φ)) 0这个表达式是关于 x 的线性函数大于 0 决策为正类小于 0 决策为负类等于 0 就是决策边界。线性函数等于 0 在二维平面上是一条直线在高维空间是一个超平面。这就是为什么 LDA 的决策边界是线性的。如果每个类别有自己的协方差矩阵 Σ_k二次项无法抵消决策边界就成了二次曲面这就是 QDA 的来历。理解这一点很重要并不是 GDA 本身决定了边界形状而是你对协方差矩阵的假设决定了边界形状。我在实际项目中还发现一个有趣的现象当数据量比较小时LDA 的线性假设虽然在偏差上稍大但方差小整体泛化能力往往优于 QDA。而数据量足够大时QDA 的灵活性就开始体现优势。这背后是机器学习里最经典的偏差-方差权衡。3. 核心实操从零实现高斯判别分析3.1 环境准备与数据构造纸上谈兵不如动手跑代码。我用 Python 从零实现一个 GDA不依赖 sklearn 的高级封装只用 NumPy 手动算参数这样每一步都透明可控。先构造一个二维数据集方便可视化。我会生成两类数据每类各 100 个样本第一类以 (1, 2) 为中心第二类以 (4, 5) 为中心协方差矩阵共享import numpy as np import matplotlib.pyplot as plt np.random.seed(42) # 生成两类数据共享协方差矩阵 cov np.array([[1.2, 0.5], [0.5, 1.0]], dtypenp.float64) # 第一类100个样本均值为 (1, 2) X0 np.random.multivariate_normal([1.0, 2.0], cov, size100) # 第二类100个样本均值为 (4, 5) X1 np.random.multivariate_normal([4.0, 5.0], cov, size100) X np.vstack([X0, X1]) y np.array([0] * 100 [1] * 100) # 打乱顺序 idx np.random.permutation(len(X)) X, y X[idx], y[idx] # 可视化 plt.figure(figsize(8, 6)) plt.scatter(X[y 0][:, 0], X[y 0][:, 1], alpha0.6, labelclass 0) plt.scatter(X[y 1][:, 0], X[y 1][:, 1], alpha0.6, labelclass 1) plt.xlabel(x1) plt.ylabel(x2) plt.legend() plt.title(Training Data) plt.grid(True, alpha0.3) plt.show()这里我故意选择了有正相关性的协方差矩阵0.5 的正相关这样两个类别的高斯分布都是斜着的椭圆决策边界的走向会比较有意思。3.2 手写GDA参数估计与预测的完整代码接下来实现 GDA 的核心逻辑每一步都对应前面推导的公式class GDA: def __init__(self): self.phi None self.mu0 None self.mu1 None self.Sigma None def fit(self, X, y): m, n X.shape # 正类样本的索引和数量 pos_idx np.where(y 1)[0] neg_idx np.where(y 0)[0] m_pos len(pos_idx) m_neg len(neg_idx) # 1. phi正类样本占比 self.phi m_pos / m # 2. mu0 和 mu1各类别均值 self.mu0 np.mean(X[neg_idx], axis0) self.mu1 np.mean(X[pos_idx], axis0) # 3. Sigma共享协方差矩阵 # 对所有样本减去其所属类别的均值再外积累加 Sigma np.zeros((n, n)) for i in range(m): if y[i] 1: diff X[i] - self.mu1 else: diff X[i] - self.mu0 Sigma np.outer(diff, diff) self.Sigma Sigma / m return self def predict_proba(self, X): # 计算 p(y1|x) 和 p(y0|x)返回正类概率 # 用 log 形式避免数值下溢 m, n X.shape # 多元高斯密度的log def log_gaussian(x, mu, Sigma_inv, logdet): diff x - mu return -0.5 * (n * np.log(2 * np.pi) logdet diff Sigma_inv diff) Sigma_inv np.linalg.inv(self.Sigma) sign, logdet np.linalg.slogdet(self.Sigma) log_p0 log_gaussian(X, self.mu0, Sigma_inv, logdet) np.log(1 - self.phi) log_p1 log_gaussian(X, self.mu1, Sigma_inv, logdet) np.log(self.phi) # log-sum-exp 归一化得到概率 log_p np.column_stack([log_p0, log_p1]) # 减去每行最大值保证 exp 不溢出 max_log np.max(log_p, axis1, keepdimsTrue) p np.exp(log_p - max_log) p p / np.sum(p, axis1, keepdimsTrue) return p[:, 1] def predict(self, X, threshold0.5): proba self.predict_proba(X) return (proba threshold).astype(int) # 训练 model GDA() model.fit(X, y) # 查看学到的参数 print(phi , model.phi) print(mu0 , model.mu0) print(mu1 , model.mu1) print(Sigma ) print(model.Sigma)predict_proba 里有个细节值得多说一句直接对高斯密度函数做数值计算如果维度高或者数据分布比较极端exp(-1/2 * diff^T Σ^(-1) diff) 可能会下溢成 0导致概率全是 0/1看起来好像模型很自信其实是数值问题。所以我用 log 空间计算再用 log-sum-exp 技巧归一化这样既稳又准。顺便看一下拟合出的参数是否接近真实的生成参数真实 φ0.5真实 μ0(1,2)真实 μ1(4,5)真实 Σ [[1.2, 0.5], [0.5, 1.0]]。因为模拟数据本身是从这些分布里采样出来的所以估计值应该非常接近。这就是 GDA 的另一个好处参数有明确的可解释性你可以直接检查学到的分布是否合理。3.3 决策边界可视化与效果评估参数学完了接下来做完整的预测和可视化。决策边界的绘制方法在二维平面上均匀撒点用模型预测每个点的类别用等高线画出 p(y1|x) 0.5 的位置这其实就是决策边界。# 在平面上生成网格点 xx, yy np.meshgrid(np.linspace(-3, 8, 300), np.linspace(-2, 9, 300)) grid np.c_[xx.ravel(), yy.ravel()] proba model.predict_proba(grid) proba proba.reshape(xx.shape) plt.figure(figsize(10, 8)) # 绘制概率热力图 plt.contourf(xx, yy, proba, levelsnp.linspace(0, 1, 20), cmapcoolwarm, alpha0.6) # 绘制决策边界 plt.contour(xx, yy, proba, levels[0.5], colorsblack, linewidths2) # 绘制数据点 plt.scatter(X[y 0][:, 0], X[y 0][:, 1], alpha0.7, labelclass 0, edgecolorsk) plt.scatter(X[y 1][:, 0], X[y 1][:, 1], alpha0.7, labelclass 1, edgecolorsk) # 绘制均值点 plt.scatter(model.mu0[0], model.mu0[1], cblue, marker*, s200, labelmu0) plt.scatter(model.mu1[0], model.mu1[1], cred, marker*, s200, labelmu1) plt.xlabel(x1) plt.ylabel(x2) plt.legend() plt.title(GDA Decision Boundary and Probability Heatmap) plt.grid(True, alpha0.3) plt.show()从热力图可以清楚看到靠近 class 0 中心的地方概率趋近 0靠近 class 1 中心的地方概率趋近 1中间有一条平滑过渡带。黑色直线就是决策边界方向大致与两个类别质心的连线垂直——这个直觉很重要GDA 的线性决策边界本质上是在两个高斯山峰之间找一个等概率分界线。训练集上的准确率评估y_pred model.predict(X) accuracy np.mean(y_pred y) print(f训练集准确率: {accuracy:.4f})理论上这个模型在测试集上的表现会不错因为数据确实是按高斯分布生成的。但在真实场景中数据生成过程往往没那么理想这时候 GDA 的表现取决于它的假设与真实情况的偏离程度。这就是为什么实际工程中我们通常会把 GDA 和逻辑回归都跑一遍对比结果再选模型。4. 常见问题与实战经验4.1 协方差矩阵奇异数据维度太高怎么办实战中第一个拦路虎就是特征维度 n 大于样本量 m。比如你做基因表达数据分类样本只有 50 个特征却有两万维。这时候估计出来的 Σ 是奇异矩阵根本求不了逆。即便 m 略大于 nΣ 也可能病态严重逆矩阵中的数值会极大。我踩过这个坑有一次做文本分类TF-IDF 特征抽取完有两千多维样本只有三百条直接跑 GDAnp.linalg.inv(Σ) 报错 LinAlgError: Singular matrix。后来学乖了常见的解决方案有这么几种降维先用 PCA 把特征降到远小于样本量的维度再套 GDA。这是最推荐的做法因为 PCA 能去相关而 GDA 的高斯假设在去相关后的数据上更合理。正则化给 Σ 的对角线加一个小的正数 λI即用 Σ λI 代替 Σ。这就是 LDA 里常见的 shrinkage 技巧。sklearn 的 LinearDiscriminantAnalysis 里 solverlsqr 和 shrinkage 参数就是干这个的。改用朴素贝叶斯假设各特征独立协方差矩阵退化为对角矩阵需要估计的参数骤减奇异问题自然消失。代价是忽略了特征间的相关性。我建议的优先级是PCA 降维 正则化 换模型。PCA 既解决了奇异性又顺手去噪一举两得。4.2 GDA与逻辑回归到底该选哪个这是面试高频题也是实际选型常纠结的问题。结论是GDA 和逻辑回归各有适用场景没有绝对的好坏取决于你对数据生成过程的信心。GDA 的优势在于当数据确实近似高斯分布时它需要更少的样本就能达到不错的效果因为它利用了很强的分布假设它还能输出完整的 P(x|y)可以做异常检测、数据生成等逻辑回归做不到的事参数可解释性强能直接给出每个类别的均值、协方差对业务分析很有价值。逻辑回归的优势在于它不对特征分布做假设鲁棒性更强。数据不是高斯分布、特征间有复杂相关性、有离群点逻辑回归照样能拟合出不错的决策边界。它直接建模决策边界在分类问题上是更直接的解法。我个人的经验法则如果样本量很小几百级别且特征经过变换后大致满足高斯分布优先试 GDA如果样本量大、特征类型复杂离散、混合类型、稀疏或者不确定分布形态直接用逻辑回归。还有一种情况我强烈推荐 GDA类不平衡严重但你想解释每个类到底长什么样的时候GDA 的均值向量就是天然的类别画像。4.3 实战中的几个坑坑一特征尺度不统一。高斯分布对尺度敏感比如一个特征取值范围是 (0, 1)另一个是 (0, 10000)协方差矩阵会严重偏向大尺度的特征。务必先做标准化或归一化让所有特征大致在同一个量级。这个操作对很多模型都有用但对 GDA 尤其重要因为高斯分布的形状完全由均值和协方差定义。坑二离群点影响大。因为高斯分布的尾部衰减很快一个远离中心的数据点会对均值估计产生极大影响。我在真实项目里发现训练集里混入三五个异常样本μ 和 Σ 的估计值就明显偏移决策边界跟着歪。处理办法在拟合前做离群点检测比如用 IQR、MAD或者直接可视化把肉眼可见的异常点剔除。坑三类别先验 φ 的使用。当类别不平衡时φ 会影响决策边界的位置让分类器偏向样本多的类别。这对某些业务可能是好事比如你希望倾向预测多数类以降低误报率但有时你并不想要这种倾向。你可以显式地把 φ 设置为 0.5忽略先验的偏差只看似然项的贡献。坑四过拟合风险。GDA 要估计的参数包括 μ0、μ1 和 Σ协方差矩阵有 n(n1)/2 个参数。特征维度稍高一点参数数量就蹭蹭涨很容易过拟合。所以高维场景记得配 PCA 和交叉验证。5. 扩展与应用GDA 的工业级用法5.1 多分类与 QDA从二分类走向更复杂场景GDA 天然支持多分类只需要把二分类的伯努利分布换成类别分布Categorical每个类别有自己的均值 μ_k共享协方差 Σ对应 LDA或者每个类别有自己的协方差 Σ_k对应 QDA。决策规则就是算每个类别的后验概率取最大者。sklearn 里的 LinearDiscriminantAnalysis 和 QuadraticDiscriminantAnalysis 就是工程化之后的 GDA。前者默认求解器是 SVD数值稳定性好不需要显式求逆后者直接按类别估计各自的协方差矩阵。我建议你在项目里优先用这些现成的实现它们已经处理了无数边界情况比自己手写的版本可靠得多。但有一点要注意sklearn 的 LDA 实现里有个 n_components 参数它可以做监督降维——把高维数据投影到一个低维子空间使得类间散度与类内散度的比值最大。这个用法的本质和 GDA 的分类是一样的数学结构但在降维场景下非常好用。5.2 异常检测与半监督学习中的应用GDA 的生成特性让它不仅能做分类还能做异常检测既然我们学到了每个类别的 P(x|y)那对一个新的样本我们可以计算它在各类别下的似然值 P(x|μ_k, Σ_k)。如果这个值低于某个阈值说明该样本与所有已知类别都不匹配很可能是一个异常或新类别。我在工业项目里就做过类似的事设备传感器数据里用正常工况的数据拟合一个高斯分布然后监控新样本的似然值一旦似然值显著下降就触发告警。这本质上就是多元高斯异常检测也是 GDA 最自然的延伸。它的最大好处是训练只需要正常数据不需要标注异常样本在异常数据极度稀缺的场景下非常实用。半监督学习的方向也值得一提假设你有少量标注数据还有大量无标注数据。可以用少量标注数据初始化 GDA 参数然后用 EM 算法迭代(E步) 对无标注数据计算它属于各类别的概率(M步) 用这些概率加权重新估计各类别的均值和协方差。这就是高斯混合模型GMM的思路而 Gaussian Mixture Model 就是 GDA 的无监督堂兄弟。所以理解了 GDA你就等于掌握了通向 GMM、EM 算法的一级台阶。5.3 特征相关性的影响一个容易被忽视的细节最后再聊一个实操中特别容易忽视的细节GDA 对特征相关性的处理方式。如果特征之间有强相关性协方差矩阵的非对角元素会很大。这时候 GDA 的线性判别边界仍然可以正常计算但参数的含义会变得微妙——某个特征的单独贡献会被其他相关特征稀释。举个例子假设 x1 和 x2 高度相关相关系数 0.9并且 x1 对区分两类有重要作用那么你单独看 μ1 和 μ0 在 x1 维度上的差异会觉得分类应该很容易但因为 x2 和 x1 几乎同步变化协方差矩阵的某些方向方差很大导致整体判别信号被削弱。这引出了一个实用建议当你用 GDA 做特征解释时不要只看均值差异要把协方差矩阵考虑进来。更稳妥的做法是先用 PCA 去相关再跑 GDA这样学出来的均值和协方差更干净决策边界也更稳定。我自己的经验是在大多数表格数据场景下标准化 PCA 去相关 GDA的组合效果往往优于直接用原始特征跑 GDA而且收敛更快、数值更稳。这一步做完你会发现自己对数据分布的理解也更深了——因为 GDA 逼着你去审视每个特征的分布形态、特征之间的相关结构这种透过模型看数据的视角是纯调包选手永远体会不到的。