机器学习算法系列（十一）-线性判别分析算法（二）（Linear Discriminant Analysis Algorithm）

阅读本文需要的背景知识点：正态分布、线性判别分析、一丢丢编程知识

一、引言

  前面一节介绍了基本的线性判别分析算法，最后留下了一个问题，我们在使用 sklearn 得到的结果与上一节中自己实现的结果不同，这一节就来看看 sklearn 中通过概率分布的角度是如何实现线性判别分析的。

二、模型介绍

一元正态分布

  在介绍模型之前先来回顾一下一个在统计学中特别常见的分布——正态分布¹（Normal distribution），同时也被称为高斯分布（Gaussian distribution），该分布是下面介绍模型的基础。
$$f(x;\mu ;\sigma )=\frac{1}{\sigma \sqrt{2\pi }}\exp \left(-\frac{(x-\mu {)}^2}{2{\sigma }^2}\right)$$

  上面为一元正态分布的概率密度函数，其中 $x$ 为一元随机变量，$\mu$ 为位置参数，$\sigma$ 为尺度参数。

不同尺度参数的正态分布图示

多元正态分布

  在机器学习中我们的输入特征一般是多维的，多元正态分布的概率密度函数可以从一元正态分布的概率密度函数得到的。为了简单起见，假设随机变量之间是独立的，即相互之间不存在线性相关性。
（1）随机变量之间相互独立，则联合概率密度函数为每一个随机变量的概率密度函数的乘积，其中 x 为 $p$ 维
（2）分别带入各个一元概率密度函数
（3）将 e 的指数用 z 的平方来表示
（4）z 的平方可以写成向量矩阵的形式
（5）观察中间的对角矩阵，用 $\Sigma$ 来表示该对角矩阵的逆矩阵
（6）则 z 的平方可以简写成该式，x 为特征向量，$\mu$ 为均值向量。该表达式被称为马哈拉诺比斯距离²（Mahalanobis distance），用于表示数据的协方差距离
（7）观察 $\Sigma$，$\Sigma$ 的行列式为对角线上的元素相乘，开根号后为各 $\sigma$ 的乘积，$|\Sigma |$ 表示 $\Sigma$ 的行列式
（8）带入回多元正态分布的概率密度函数中，得到最终的形式
$$\begin{array}{rlr}f(x)& =f\left(x_1\right)f\left(x_2\right)\dots f\left(x_p\right)& (1)\\ & =\frac{1}{{\sigma }_1{\sigma }_2\dots {\sigma }_p(\sqrt{2\pi }{)}^p}\exp \left(-\frac{{\left(x_1-{\mu }_1\right)}^2}{2{\sigma }_1^2}-\frac{{\left(x_2-{\mu }_2\right)}^2}{2{\sigma }_2^2}-\cdots -\frac{{\left(x_p-{\mu }_p\right)}^2}{2{\sigma }_p^2}\right)& (2)\\ z^2& =\frac{{\left(x_1-{\mu }_1\right)}^2}{{\sigma }_1^2}+\frac{{\left(x_2-{\mu }_2\right)}^2}{{\sigma }_2^2}+\cdots +\frac{{\left(x_p-{\mu }_p\right)}^2}{{\sigma }_p^2}& (3)\\ & ={\left[\begin{array}{c}{x}_1-{\mu }_1\\ x_2-{\mu }_2\\ \cdots \\ x_p-{\mu }_p\end{array}\right]}^T\left[\begin{array}{cccc}\frac{1}{{\sigma }_1^2}& 0& \cdots & 0\\ 0& \frac{1}{{\sigma }_2^2}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & \frac{1}{{\sigma }_p^2}\end{array}\right]\left[\begin{array}{c}{x}_1-{\mu }_1\\ x_2-{\mu }_2\\ \cdots \\ x_p-{\mu }_p\end{array}\right]& (4)\\ \Sigma & =\left[\begin{array}{cccc}{\sigma }_1^2& 0& \cdots & 0\\ 0& {\sigma }_2^2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\sigma }_p^2\end{array}\right]& (5)\\ z^2& =(x-\mu {)}^T{\Sigma }^{-1}(x-\mu )& (6)\\ |\Sigma {|}^{\frac{1}{2}}& ={\sigma }_1{\sigma }_2\dots {\sigma }_p& (7)\\ f(x)& =\frac{1}{|\Sigma {|}^{\frac{1}{2}}(2\pi {)}^{\frac{p}{2}}}\exp \left(-\frac{(x-\mu {)}^T{\Sigma }^{-1}(x-\mu )}{2}\right)& (8)\end{array}$$

  以上为相互独立的随机变量的多元正态分布的概率密度函数，对于任意的随机变量，该式亦成立，只是Σ表示随机变量的协方差矩阵。

贝叶斯定理

  再来回顾一下概率论中的一个重要定理——贝叶斯定理³，在 B 事件发生时，A 事件发生的概率等于在 A 事件发生时，B 事件发生的概率乘以 A 事件发生的概率除以 B 事件发生的概率，用公式表示如下：
$$P(A\mid B)=\frac{P(B\mid A)P(A)}{P(B)}$$

线性判别分析

  前面铺垫了这么多前置知识，下面来介绍以概率的角度来实现线性判别分析的方法。首先我们需要假设样本点是服从正态分布，并且每一类样本的协方差矩阵相同，满足上面两点我们就可以开始推导线性判别分析的实现过程。
（1）我们先来看下在输入为 x 的情况下分类为k的概率，根据前面介绍的贝叶斯定理可得到下面的（1）式
（2）样本点服从正态分布，在分类为 k 的情况下，输入为 x 的概率可以用多元正态分布的概率密度函数表示
$$\begin{array}{rlr}P(k\mid x)& =\frac{P(x\mid k)P(k)}{P(x)}& (1)\\ & =\frac{f_k(x)P(k)}{P(x)}& (2)\end{array}$$

（1）我们的目的就是求在输入为 x 的情况下分类为 k 的概率最大的分类，所以我们可以写出假设函数如下图（1）式，该方法被称为最大后验概率估计⁴（maximum a posteriori probability (MAP) estimate）
（2）对其概率取对数，不影响函数的最后结果
（3）带入上面的 $P(k|x)$ 的表达式，由于 $P(x)$ 对最后结果也没有影响，也可以直接去掉
（4）带入多元正态分布的概率密度函数表达式
（5）将（4）式中的对数化简得到
（6）观察（5）式第二项，由于每一类样本的协方差矩阵相同，不影响最后的结果，可以直接去掉
$$\begin{array}{rlr}h(x)& =\underset{k}{\mathrm{argmax}}P(k\mid x)& (1)\\ & =\underset{k}{\mathrm{argmax}}\ln P(k\mid x)& (2)\\ & =\underset{k}{\mathrm{argmax}}\ln f_k(x)+\ln P(k)& (3)\\ & =\underset{k}{\mathrm{argmax}}\ln \left(\frac{e^{-\frac{{\left(x-{\mu }_k\right)}^T{\Sigma }_k^{-1}\left(x-{\mu }_k\right)}{2}}}{{\left|{\Sigma }_k\right|}^{\frac{1}{2}}(2\pi {)}^{\frac{p}{2}}}\right)+\ln P(k)& (4)\\ & =\underset{k}{\mathrm{argmax}}-\frac{1}{2}{\left(x-{\mu }_k\right)}^T{\Sigma }_k^{-1}\left(x-{\mu }_k\right)-\ln \left({\left|{\Sigma }_k\right|}^{\frac{1}{2}}(2\pi {)}^{\frac{p}{2}}\right)+\ln P(k)& (5)\\ & =\underset{k}{\mathrm{argmax}}-\frac{1}{2}{\left(x-{\mu }_k\right)}^T{\Sigma }_k^{-1}\left(x-{\mu }_k\right)+\ln P(k)& (6)\end{array}$$

（1）我们再来看下上面（6）式中的第一项，将其展开
（2）观察第三项，可以合并成整体的转置的形式，注意协方差矩阵 $\Sigma$ 为对称矩阵，其转置等于其本身
（3）上面第二三项互为转置，结果又是一个实数，所以其结果相等，可以直接合并为一项
$$\begin{array}{rlr}{\left(x-{\mu }_k\right)}^T{\Sigma }^{-1}\left(x-{\mu }_k\right)& =x^T{\Sigma }^{-1}x-x^T{\Sigma }^{-1}{\mu }_k-{\mu }_k^T{\Sigma }^{-1}x+{\mu }_k^T{\Sigma }^{-1}{\mu }_k& (1)\\ & =x^T{\Sigma }^{-1}x-x^T{\Sigma }^{-1}{\mu }_k-{\left(x^T{\Sigma }^{-1}{\mu }_k\right)}^T+{\mu }_k^T{\Sigma }^{-1}{\mu }_k& (2)\\ & =x^T{\Sigma }^{-1}x-2x^T{\Sigma }^{-1}{\mu }_k+{\mu }_k^T{\Sigma }^{-1}{\mu }_k& (3)\end{array}$$

（1）前面所得的假设函数
（2）用上面展开的结果替换第一项
（3）在（2）式中第一项的结果不影响最后的结果，也可以直接去掉，得到最后的假设函数
$$\begin{array}{rlr}h(x)& =\underset{k}{\mathrm{argmax}}-\frac{1}{2}{\left(x-{\mu }_k\right)}^T{\Sigma }^{-1}\left(x-{\mu }_k\right)+\ln P(k)& (1)\\ & =\underset{k}{\mathrm{argmax}}-\frac{1}{2}{x}^T{\Sigma }^{-1}x+x^T{\Sigma }^{-1}{\mu }_k-\frac{1}{2}{\mu }_k^T{\Sigma }^{-1}{\mu }_k+\ln P(k)& (2)\\ & =\underset{k}{\mathrm{argmax}}{x}^T{\Sigma }^{-1}{\mu }_k-\frac{1}{2}{\mu }_k^T{\Sigma }^{-1}{\mu }_k+\ln P(k)& (3)\end{array}$$

（1）将假设函数的函数部分写成新的一个函数 L(x)
（2）将其中两类做差，当结果等于零时，即此时既可以将其分为 i 类，也可以分为 j 类。其本质为分类 i 与分类 j 的决策边界
（3）可以看到该决策边界是一个关于 x 的线性函数
$$\begin{array}{rlr}{L}_k(x)& =x^T{\Sigma }^{-1}{\mu }_k-\frac{1}{2}{\mu }_k^T{\Sigma }^{-1}{\mu }_k+\ln P(k)& (1)\\ L_i(x)-L_j(x)& =x^T{\Sigma }^{-1}\left({\mu }_i-{\mu }_j\right)-\frac{1}{2}\left({\mu }_i^T{\Sigma }^{-1}{\mu }_i-{\mu }_j^T{\Sigma }^{-1}{\mu }_j\right)+\ln P(i)-\ln P(j)& (2)\\ & =x^{T}w+b& (3)\end{array}$$

（1）、（2）根据上面的决策边界的线性函数形式，可以得到对应的 w 和 b
（3）$P(k)$ 表示 k 分类的先验概率，即 k 分类的样本数除以总样本数
（4）$\mu$ 表示均值向量
（5）$\Sigma$ 为协方差矩阵的估计，因为需要使得每一个分类的协方差矩阵相同，这里将每一个分类的协方差矩阵相加后使用 $N-K$ 来做归一化处理
$$\begin{array}{rlr}w& ={\Sigma }^{-1}\left({\mu }_i-{\mu }_j\right)& (1)\\ b& =-\frac{1}{2}\left({\mu }_i^T{\Sigma }^{-1}{\mu }_i-{\mu }_j^T{\Sigma }^{-1}{\mu }_j\right)+\ln P(i)-\ln P(j)& (2)\\ P(k)& =\frac{N_k}{N}& (3)\\ {\mu }_k& =\frac{1}{N_k}\sum _{i=1}^{N_k}x& (4)\\ \Sigma & =\frac{1}{N-K}\sum _{i=1}^K\left(x-{\mu }_i\right){\left(x-{\mu }_i\right)}^T& (5)\end{array}$$

  当需要判断新样本点的分类时只需将样本点带入该线性函数中，对于二分类来说，最后的结果大于零，则分类为 i，反之则分类为 j。

两种角度下的线性判别分析的关系

  两种角度下的线性判别分析的权重系数的公式如下：
$$\begin{array}{cc}w=S_w^{-1}\left({\mu }_1-{\mu }_2\right)& (1)\\ w={\Sigma }^{-1}\left({\mu }_1-{\mu }_2\right)& (2)\end{array}$$

  当每一个分类的协方差矩阵相同的时候，会有 $S_w$ 等于 $2\Sigma$，说明其权重系数的方向是相同的，当样本服从正态分布时且分类的协方差矩阵相同时，两种角度下线性判别分析的结果是相同的。

三、代码实现

使用 Python 实现线性判别分析（LDA）：

def ldaBayes(X, y):
   """
   线性判别分析（LDA）
   args:
       X - 训练数据集
       y - 目标标签值
   return:
       W - 权重系数
       b - 偏移量
   """
   # 标签值
   y_classes = np.unique(y)
   # 第一类
   c1 = X[y==y_classes[0]][:]
   # 第二类
   c2 = X[y==y_classes[1]][:]
   # 第一类均值向量
   mu1 = np.mean(c1, axis=0)
   # 第二类均值向量
   mu2 = np.mean(c2, axis=0)
   # 第一类协方差矩阵
   sigma1 = c1 - mu1
   sigma1 = sigma1.T.dot(sigma1)
   # 第二类协方差矩阵
   sigma2 = c2 - mu2
   sigma2 = sigma2.T.dot(sigma2)
   # Σ 矩阵
   sigma = (sigma1 + sigma2) / (X.shape[0] - len(y_classes))
   # Σ 逆矩阵
   sigman = np.linalg.pinv(sigma)
   # 求权重系数
   w = sigman.dot(mu2 - mu1)
   b = -0.5 *(mu2.reshape(-1, 1).T.dot(sigman).dot(mu2) - mu1.reshape(-1, 1).T.dot(sigman).dot(mu1)) + np.log(len(c2) / len(c1))
   return w, b

def discriminantBayes(X, w, b):
   """
   判别新样本点
   args:
       X - 数据集
       w - 权重系数
       b - 偏移量
   return:
       分类结果
   """
   r = X.dot(w) + b
   r[r > 0] = 1
   r[r <= 0] = 0
   return r

四、第三方库实现

scikit-learn⁵ 实现线性判别分析：

from sklearn.discriminant_analysis import LinearDiscriminantAnalysis

# 初始化线性判别分析器
lda = LinearDiscriminantAnalysis()
# 拟合线性模型
lda.fit(X, y)
# 权重系数
w = lda.coef_
# 截距
b = lda.intercept_

五、思维导图

六、参考文献

1. https://en.wikipedia.org/wiki/Normal_distribution
2. https://en.wikipedia.org/wiki/Mahalanobis_distance
3. https://en.wikipedia.org/wiki/Bayes%27_theorem
4. https://en.wikipedia.org/wiki/Maximum_a_posteriori_estimation
5. https://scikit-learn.org/stable/modules/generated/sklearn.discriminant_analysis.LinearDiscriminantAnalysis.html

完整演示请点击这里

注：本文力求准确并通俗易懂，但由于笔者也是初学者，水平有限，如文中存在错误或遗漏之处，恳请读者通过留言的方式批评指正

本文首发于——AI导图，欢迎关注