机器学习（九）——线性判别分析(LDA)

1 线性判别原理

1.1 基于贝叶斯公式的理解

基于贝叶斯即基于分类思想，一个常见的LDA分类基本思想是：假设各个类别的样本数据符合高斯分布，这样利用LDA进行投影后，可以利用极大似然估计计算各个类别投影数据的均值和方差，进而得到该类别高斯分布的概率密度函数。当一个新的样本到来后，我们可以将它投影，然后将投影后的样本特征分别带入各个类别的高斯分布概率密度函数，计算它属于这个类别的概率，最大的概率对应的类别即为预测类别。

对于多分类问题数据集为 X = { x 1 , x 2 , … , x n } X=\{x_1,x_2,\dots,x_n\} X={x1​,x2​,…,xn​}，类别标签为 Y = { y 1 , y 2 , … , y n } Y=\{y_1,y_2,\dots,y_n\} Y={y1​,y2​,…,yn​}，共有$K$类，参考贝叶斯公式我们可以将贝叶斯公式进一步写成：
$$P(Y=k|X=x)=\frac{P(Y=k)P(X=x|Y=k)}{{\sum }_{y\in Y}P(y)P(X=x|y)}$$
$P(Y=k)$为观测样本来自第$k$类的先验概率，即$P(Y=k)=\frac{n(Y=k)}{n}$，$n(Y=k)$表示所有样本中属于第$k$类的样本的数量，$n$为所有样本的数量，也可以理解为第$k$类样本占所有样本的比例。$P(X=x|Y=k)=\frac{n(X=x,Y=k)}{n(Y=k)}$表示当类别为$k$时，观测样本为$x$的概率密度，即在类别为$k$的所有样本中样本为$x$所占的比例。因此分子可以化简为$\frac{n(X=x,Y=k)}{n}$。对于分母，可化简为$P(X=x)=\frac{n(X=x)}{n}$，也就是所有观测样本中样本$x$所占的比例。贝叶斯公式整体可以化简为$\frac{n(X=x,Y=k)}{n(X=x)}$，即观测样例$x$中类别属于$k$的样例所占的比例。

该公式给出了在给定样本条件下，$Y=k$这个类别的概率。但是当样本给定时，分母$P(X=x)$就是一个与类别$k$无关的常数，因此问题可以简化为计算分子，进而比较哪个类别最大就可以判断出样本属于哪个类别。因此对于分类问题，计算贝叶斯定理的分子，分子最大的一类为样本所属类别。

下面我们基于此来推导LDA算法。此时我们令${\pi }_k=P(Y=k)$，$f_k(x)=P(X=x|Y=k)$，令分子为$g_k=\pi f_k(x)$。
（1）我们先来推导一下自变量只有一个的情况，即$x$是一维的，只有一个特征。我们假设各个类别的样本数据服从正态分布，即$f_k(x)\sim N(\mu ,{\sigma }_k^2)$，并且每个${\sigma }_k^2={\sigma }^2$，同方差假设。那么
$$f_k(x)=\frac{1}{\sqrt{2\pi {\sigma }^2}}{e}^{-\frac{1}{2{\sigma }^2}|x-{\mu }_k{|}^2}$$
那么分子最终为：
$$g_k={\pi }_k\frac{1}{\sqrt{2\pi {\sigma }^2}}{e}^{-\frac{1}{2{\sigma }^2}|x-{\mu }_k{|}^2}$$
由于此式较为复杂，难以求解，因此我们对其取对数：
$$\begin{gathered} {\delta }_k=ln(g_k)=ln({\pi }_k)-\frac{1}{2}ln(2\pi {\sigma }^2)-\frac{1}{2{\sigma }^2}|x-{\mu }_k{|}^2 \\ =ln({\pi }_k)-ln(\sigma )-\frac{x^2}{2{\sigma }^2}+\frac{{\mu }_k}{{\sigma }^2}x-\frac{{\mu }_k^2}{2{\sigma }^2} \\ =ln({\pi }_k)+\frac{{\mu }_k}{{\sigma }^2}x-\frac{{\mu }_k^2}{2{\sigma }^2} \end{gathered}$$
至此，我们需要做的就是将${\mu }_k$和$\sigma$估计出来：
$$\begin{gathered} \widehat{{\mu }_k}=\frac{1}{n_k}\sum _{i,y_i=k}{x}_i \\ \hat{\sigma }=\frac{1}{n-K}\sum _{k=1}^K\sum _{i,y_i=k}(x_i-{\mu }_k{)}^2 \end{gathered}$$
即求$y=k$类中$x$的均值，和每一类的方差再求均值。
因此LDA最终可表示为：
$$\{\begin{array}{l}{\delta }_k(x)=ln(g_k(x))=ln{\pi }_k+\frac{{\mu }_k}{{\sigma }^2}x-\frac{{\mu }_k^2}{2{\sigma }^2}\\ {\hat{\mu }}_k=\frac{1}{n_k}\sum _{i:y_i=k}{x}_i\\ {\hat{\sigma }}^2=\frac{1}{n-K}\sum _{k=1}^K\sum _{i:y_i=k}(x_i-{\hat{\mu }}_k{)}^2\end{array}$$
（2）对于多个自变量的情况，即$x$是d维的。只需要将将一元正态分布扩展为多元正态分布$f_k(x)=\frac{1}{(2\pi {)}^{\frac{d}{2}}|\Sigma {|}^{\frac{1}{2}}}{e}^{[-\frac{1}{2}(x-{\mu }_k{)}^T{\Sigma }^{-1}(x-{\mu }_k)]}$
$\widehat{{\mu }_k}=({\mu }_{k1},{\mu }_{k2},......,{\mu }_{kd}),\hat{\Sigma }=\frac{1}{d-1}\sum _{j=1}^d(x_j-\overline{x})(x_j-\overline{x}{)}^T$
${\delta }_k(x)=ln({\pi }_k{f}_k(x))=ln({\pi }_k)-(\frac{d}{2}ln(2\pi )+\frac{1}{2}ln(|\Sigma |))-\frac{1}{2}(x-{\mu }_k{)}^T{\Sigma }^{-}1(x-{\mu }_k)=x^T\hat{\Sigma }{\hat{\mu }}_k-\frac{1}{2}{\hat{\mu }}_k^T{\hat{\Sigma }}^{-1}{\hat{\mu }}_k+ln{\hat{\pi }}_k$

1.2 基于降维分类思想的理解

LDA是一种监督学习的降维技术，也就是说它的数据集的每个样本是有类别输出的。其思想非常简单：对于给定数据集，设法将其投影到一条直线上，使得同类样本的投影点尽可能接近，异类样本的投影点尽可能远离；对新的测试样本，将其投影到同样的直线上，在根据投影点的位置来判断测试样本的类别。“投影后类内方差最小，类间方差最大”

由于协方差不仅可以反映变量之间的相关性，同样可以反映多维样本分布的离散程度（一维样本使用方差），协方差越大（对于负相关来说是绝对值越大），表示数据的分布越分散。
$$J(w)=\frac{w^T|{\mu }_1-{\mu }_2^{}{|}^2}{s_1^2+s_2^2}$$
分子为投影数据后的均值之差，分母为方差之和，LDA的目的就是使得$J$值最大化。可以理解为最大化分子，即使得类别之间的距离越远，同时最小化分母，使得每个类别内部的方差越小。这样就能使得每个类类别的数据可以在投影矩阵$w$的映射下分得越开。

数据集特征为$X=(x_1,x_2,\dots ,x_n{)}^T$，类别标签为$Y=(y_1,y_2,\dots ,y_n{)}^T$，其中 y i ∈ { + 1 , − 1 } y_i \in \{+1,-1\} yi​∈{+1,−1}，类别$c_1$的特征为 X c 1 = { x i ∣ y i = + 1 } X_{c_1} = \{x_i|y_i=+1 \} Xc1​​={xi​∣yi​=+1}，同理$c_2$的特征为 X c 2 = { x i ∣ y i = − 1 } X_{c_2} = \{x_i|y_i=-1 \} Xc2​​={xi​∣yi​=−1}，属于$c_1$的样本个数为$N_1$，属于$c_2$的样本个数为$N_2$，且$N=N_1+N_2$。
投影前样本的均值为
$\bar{x}=\frac{1}{N}{\sum }_{i=1}^N{x}_i$；

${\bar{x}}_1=\frac{1}{N_1}{\sum }_{i=1}^{N_1}{x}_i$；

${\bar{x}}_2=\frac{1}{N_2}{\sum }_{i=1}^{N_2}{x}_i$；

投影前样本方差为：
$\sum ={\bar{s}}^2=\frac{1}{N}{\sum }_{i=1}^N(x_i-\bar{x})(x_i-\bar{x}{)}^T$

${\sum }_1={\bar{s}}_1^2=\frac{1}{N_1}{\sum }_{i=1}^{N_1}(x_i-{\bar{x}}_1)(x_i-{\bar{x}}_1{)}^T$

${\sum }_2={\bar{s}}_2^2=\frac{1}{N_2}{\sum }_{i=1}^{N_2}(x_i-{\bar{x}}_2)(x_i-{\bar{x}}_2{)}^T$

设投影矩阵为$w$，则特征经过投影后$z=w^{T}x$。那么:

全体样本投影后的均值为$\mu =\frac{1}{N}{\sum }_{i=1}^N{w}^T{x}_i=w^T\left(\frac{1}{N}{\sum }_{i=1}^N{x}_i\right)=w^T\bar{x}$

全体样本投影后的协方差为$s^2=\frac{1}{N}{\sum }_{i=1}^N(w^T{x}_i-\mu )(w^T{x}_i-\mu {)}^T=\frac{1}{N}{\sum }_{i=1}^N(w^T{x}_i-w^T\bar{x})(w^T{x}_i-w^T\bar{x}{)}^T=\frac{1}{N}{\sum }_{i=1}^N{w}^T(x_i-\bar{x})(x_i-\bar{x}{)}^{T}w=w^T\left(\frac{1}{N}{\sum }_{i=1}^N(x_i-\bar{x})(x_i-\bar{x}{)}^T\right)w=w^T\sum w$

类别为$c_1$的样本投影后的均值为${\mu }_1=\frac{1}{N_1}{\sum }_{i=1}^{N_1}{w}^T{x}_i=w^T{\bar{x}}_1$

类别为$c_1$的样本投影后的协方差为$s_1^2=\frac{1}{N_1}{\sum }_{i=1}^{N_1}(w^T{x}_i-{\mu }_1)(w^T{x}_i-{\mu }_1{)}^T=w^T{\sum }_{1}w$

类别为$c_2$的样本投影后的均值为${\mu }_2=\frac{1}{N_2}{\sum }_{i=1}^{N_2}{w}^T{x}_i=w^T{\bar{x}}_2$

类别为$c_2$的样本投影后的协方差为$s_2^2=\frac{1}{N_2}{\sum }_{i=1}^{N_2}(w^T{x}_i-{\mu }_2)(w^T{x}_i-{\mu }_2{)}^T=w^T{\sum }_{2}w$

投影后总体类间方差为$({\mu }_1-{\mu }_2{)}^2=(w^T{\bar{x}}_1-w^T{\bar{x}}_2{)}^2=w^T({\bar{x}}_1-{\bar{x}}_2)({\bar{x}}_1-{\bar{x}}_2{)}^{T}w$

总体类内离散度为$s_1^2+s_2^2=w^T{\sum }_{1}w+w^T{\sum }_{2}w=w^T({\sum }_1+{\sum }_2)w$

由于我们要最大类间方差，最小化类内离散度，因此可以定义优化函数如下：

$J(w)=\frac{({\mu }_1-{\mu }_2{)}^2}{s_1^2+s_2^2}=\frac{w^T({\bar{x}}_1-{\bar{x}}_2)({\bar{x}}_1-{\bar{x}}_2{)}^{T}w}{w^T({\sum }_1+{\sum }_2)w}$

其中${\bar{x}}_1,{\bar{x}}_2,{\bar{s}}_1^2,{\bar{s}}_2^2$分别为未投影样本的均值和协方差。

当类间方差变大时分子增大，类内离散度减小时，分母减小。因此我们需要最大化$J(w)$。

$\hat{w}=\underset{w}{\mathrm{argmax}}J(w)$

让$J(w)$对$w$求导，然后等于0，求出$\hat{w}$：

$\hat{w}=({\sum }_1+{\sum }_2{)}^{-1}({\bar{x}}_1-{\bar{x}}_2)$

由于b不在目标函数中，所以手动去找，一般经验是：
$b=-\frac{1}{2}(w^T{\bar{x}}_1+w^T{\bar{x}}_2)$

2 代码实践

现在我们手动计算一个LDA模型：
生成两类样本数据

import numpy as np
# 2dim vectors
x_0 = np.array([[2.95,6.63],[2.45,6.13],[3.57,7.79],[3.16,5.57]])
x_1 = np.array([[2.58,4.46],[2.16,6.22],[3.27,3.52]])
print('x_0:\n',x_0)
print('x_1:\n',x_1)

x_0:
[[2.95 6.63]
[2.45 6.13]
[3.57 7.79]
[3.16 5.57]]
x_1:
[[2.58 4.46]
[2.16 6.22]
[3.27 3.52]]

然后分别计算它们的均值和方差

# mean vectors
u_0 = np.mean(x_0,axis=0)
u_1 = np.mean(x_1,axis=0)
print('u_0:{}'.format(u_0))
print('u_1:{}'.format(u_1))

u_0:[3.0325 6.53 ]
u_1:[2.67 4.73333333]

s_0 = np.cov(x_0,rowvar=False)
s_1 = np.cov(x_1,rowvar=False)
print('s_0:\n',s_0)
print('s_1:\n',s_1)

s_0:
[[0.21709167 0.25986667]
[0.25986667 0.89306667]]
s_1:
[[ 0.3141 -0.7308 ]
[-0.7308 1.87853333]]

根据推导的计算w的公式计算w，并设定b：

w = np.matmul(np.linalg.inv(s_0+s_1),(u_0-u_1))
print('w:',w)
b = -1*np.dot((u_0+u_1)/2,w)
print('b:',b[0])

print('投影直线为：y={0:.2f}x+{1:.2f}'.format(w,b))

w: [1.48009256 0.89972997]
b: -4.220113917448374
投影直线为：y=y=[1.48009256 0.89972997]x-9.287093210053445

现在我们就可以利用计算好的模型做一个分类问题了

test = np.array([[2.95,6.63],[2.53,7.79],[3.57,5.65],[3.16,5.47],[2.58,4.46],[2.16,3.11],[3.27,3.33]])
label =  [1,   1,   1,   1,   0,   0,   0]

print('决策值\tlabel')
for x in test:
    print('{:.2f}'.format(np.dot(w,x)+b),end='\t')
    if np.dot(w,x)+b>0:
        print(1)
    else:
        print(0)

决策值 label
1.04    \quad\; 1
1.47    \quad\; 1
1.08    \quad\; 1
0.31    \quad\; 1
-1.46 \quad 0
-3.29 \quad 0
-1.45 \quad 0

参考
Kaggle开源内容
https://mp.weixin.qq.com/s/pOOL-AOrD94p89CfFTzlxA