python3 Fisher线性判别分析(LDA)(含详细推导和代码）

最新推荐文章于 2024-04-26 22:00:43 发布

我不爱机器学习

最新推荐文章于 2024-04-26 22:00:43 发布

阅读量1w

点赞数 16

分类专栏：数据分析 python

本文链接：https://blog.csdn.net/mengjizhiyou/article/details/103309372

版权

数据分析同时被 2 个专栏收录

25 篇文章 10 订阅

订阅专栏

python

24 篇文章 5 订阅

订阅专栏

1、线性判别原理

线性判别分析是常用的降维技术,在模式分类和机器学习的预处理步骤中。其目标是将数据集投影到具有良好的类可分性的低维空间中，以避免过度拟合（维数过多）并降低计算成本，如将一个特征空间(一个数据集n维样本)投射到一个更小的子空间k(其中k ≤n-1)上，同时维护类区分信息。

判别式是一个函数，它接受一个输入向量x，并把它赋值给K个类中的一个，记作 C_k 。在这一章中，我们将把注意力限制在线性判别式上，即那些决策面是超平面的判别式。为了简化讨论，我们首先考虑两个类的情况，然后研究K > 2类的扩展。

1.1 二分类

假设判别面的判别式为： y(x)=w^Tx+w_0 （因为是线性所以可以直接这么写）

若 $y(x)\geq 0$ ,则向量x属于类别 C_1 ；若 y(x)< 0 ,则向量x属于类别 C_2 。所以决策边界由 y(x)= 0 决定。

假设有两个点 x_A 和 x_B 都在决策面上，则 $y(x_A)=y(x_B)=0\Rightarrow w^T(x_A-x_B)=0$ ,所以向量 w与决策面上的向量正交且w决定了决策面的方向。

例：假设有三个点 $(x_1,x_2),(x_{11},x_{22}),(x_{111},x_{222})$ 分别在判别面上方，判别面上，判别面下方。则：

$y_1=w_0\times x_0+ w_1\times x1+w_2\times x_2(x_0=1)$

$y_2=w_0\times x_0+ w_1\times x_{11}+w_2\times x_{22}(x_0=1)$

$y_3=w_0\times x_0+ w_1\times x_{111}+w_2\times x_{222}(x_0=1)$

令 y_2=0 , 则 y1>0,y3<0 ,因此决策边界的方程式也可以看作原数据点到决策面的距离。

原点 (x_1,x_2)=(0,0) 到决策面的距离为(用的决策面上的点)：

$||\overrightarrow{(x_{11},x_{22})}||\times \cos \theta$

$=\frac{||\overrightarrow{(x_{11},x_{22})}||\times \cos \theta \times||\overrightarrow{w}||}{||\overrightarrow{w}||}$

$=\frac{\overrightarrow{(x_{11},x_{22})} \times\overrightarrow{w}}{||\overrightarrow{w}||}$

$=\frac{x_{11}\times w_1+x_{22}\times w_2}{||\overrightarrow{w}||} =\frac{w^Tx}{||\overrightarrow{w}||}=\frac{-w_0}{||\overrightarrow{w}||}$

若向量x是决策面外的点，则点到面的距离：

如上图所示，向量x在决策面y=0上的投影点为向量 ${\color{Red} x_\perp}$ ,与决策面正交，所以与权重向量w平行，因此向量x可以表示为：

$x = x_\perp+\gamma \frac{w}{||w||}$

$w^Tx = w^T(x_\perp+\gamma \frac{w}{||w||})$

$w^Tx = w^Tx_\perp+\gamma w^T\frac{w}{||w||}$

因为 $y(x)=w^Tx+w_0,y(x_\perp )=w^Tx_\perp +w_0=0$ ，所以：

$y(x)-w_0=-w_0+\gamma w^T\frac{w}{||w||}$

$y(x)=\gamma \frac{w_0^2+w_1^2+w_2^2}{\sqrt{w_0^2+w_1^2+w_2^2}}$

$y(x)=\gamma ||w||\Rightarrow \gamma=\frac{y(x)}{||w||}$

简化，便于后面计算：

Y(X)=W^TX(W=[w_0, w],X=[x_0,x])

注：这里的向量默认的都是列向量,单位法向量是长度为1的法向量=向量/向量的模

1.2 多分类

现在考虑k>2类的情况：

一对多（one-versus-the-rest），复杂度K，会导致模糊区域，即数据不属于任何类别
一对一（one-versus-one），复杂度K(K - 1)/2，也会导致模糊区域

可以通过考虑一个包含K个线性函数的K类判别器来避免这些问题，则K个类别的判别式形式可以写成：

$y_k(x)=w_k^Tx+w_{k0}$

如果 $y_k(x)>y_j(x)(k\neq j)$ ,则点x属于类别 C_k 。因此类别 C_k 和类别 C_j 的决策面是 y_k(x)=y_j(x) ，所以对于超平面可以定义为：

$(w_k-w_j)^Tx+(w_{k0}-w_{j0})=0$

上图中， x_A,x_B 均属于类别 C_k ，则两点连线上的任意点 $\hat{x}$ 可以表示为：

$\hat{x} = \lambda x_A+(1-\lambda) x_B,0\leq \lambda\leq 1$

对于线性判别式，则有：

$y_k(\hat{x}) = \lambda y_k(x_A)+(1-\lambda) y_k(x_B)$

2 、Fisher线性判别分析方法

两个类别：

通过调整权重向量w，找到类分离的最大投影。假设类别 C_1 有 N_1 个点，类别 C_2 有 N_2 个点，则两个类别的平均向量为：

$\mathbf{m_1} = \frac{1}{N_1}\sum_{n\in C_1}x_n, \mathbf{m_2} = \frac{1}{N_2}\sum_{n\in C_2}x_n$

度量类分离最简单的方法是投影类到w上(看下图理解），说明可以选择w来最大化： $m_2-m_1=w^T(\mathbf{m_2}-\mathbf{m_1})$

左边的图显示了来自两个类的样本(用红色和蓝色表示)，以及在连接类的直线上投影得到的直方图。请注意，在投影空间中有相当多的类重叠。右边的图显示了基于Fisher线性判别式的对应投影，显示了极大地改进了类分离。

类别 C_k 投影数据的均值可以表示为： $m_k=w^T\mathbf{m_k}$ ，这个表达式可以通过简单地增加w的大小而变得任意大，为此可以让向量w的长度为1，则有 $\sum_iw_i^2=1$ 。通过引入拉格朗日乘子(拉格朗日乘子)来控制最大化，发现 $w\propto (\mathbf{m_2}-\mathbf{m_1})$ 。

$\max_{w}w^T(\mathbf{m_2}-\mathbf{m_1})$

$\sum_iw_i^2=1\Rightarrow w^Tw-1=0$

根据上面两个式子以及拉格朗日乘子得到： $(\mathbf{m_2}-\mathbf{m_1})=\lambda w\Rightarrow w\propto (\mathbf{m_2}-\mathbf{m_1})$

然而，这种方法仍然存在一个问题，如上图所示。这显示了两个类在原始的二维空间(x1, x2)中被很好地分离，但是当它们被投影到连接它们的平均值的直线上时，它们有相当大的重叠。

这种困难来自于类分布的强非对角协方差。Fisher提出的思想是最大化一个函数，该函数将在预计的类均值之间提供一个大的分离，同时在每个类中提供一个小的方差，从而最小化类重叠。

类别 C_k 转换后的数据即投影后的数据y的类内方差是： $s_k^2=\sum_{n\in c_k}(y_n-m_k)^2,(y_n=w^T\mathbf{x_n})$

整个数据集的类内协方差可以定义为： s_1^2+s_2^2 , Fisher准则定义类间协方差与类内协方差的比值： $J(w)=\frac{(m_2-m_1)^2}{s_1^2+s_2^2}$

根据与w的关系，可以写为：

$J(w)=\frac{(m_2-m_1)^2}{s_1^2+s_2^2}$

$= \frac{(w^T\mathbf{m_2}-w^T\mathbf{m_1})^2}{\sum_{n\in c_1}(y_n-m_1)^2+\sum_{n\in c_2}(y_n-m_2)^2}$

$= \frac{(w^T\mathbf{m_2}-w^T\mathbf{m_1})^2}{\sum_{n\in c_1}(w^T\mathbf{x_n}-w^T\mathbf{m_1})^2+\sum_{n\in c_2}(w^T\mathbf{x_n}-w^T\mathbf{m_2})^2}$

$= \frac{w^T(\mathbf{m_2}-\mathbf{m_1})(\mathbf{m_2}-\mathbf{m_1})^Tw}{\sum_{n\in c_1}(w^T\mathbf{x_n}-w^T\mathbf{m_1})^2+\sum_{n\in c_2}(w^T\mathbf{x_n}-w^T\mathbf{m_2})^2}$

$= \frac{w^T(\mathbf{m_2}-\mathbf{m_1})(\mathbf{m_2}-\mathbf{m_1})^Tw}{\sum_{n\in c_1}w^T(\mathbf{x_n}-\mathbf{m_1})(\mathbf{x_n}-\mathbf{m_1})^Tw+\sum_{n\in c_2}w^T(\mathbf{x_n}-\mathbf{m_2})(\mathbf{x_n}-\mathbf{m_2})^Tw}$

$= \frac{w^TS_Bw}{w^TS_Ww}$

对 w求导，得：

(w^TS_Bw)S_Ww=(w^TS_Ww)S_Bw

$\lambda S_Ww= S_Bw$

$\lambda S_W^{-1}S_Ww=S_W^{-1} S_Bw$

$\lambda w= S_W^{-1} S_Bw$

lambda 是特征值，w是特征向量。另一种理解方式：

$\min_w -w^TS_Bw$

w^TS_Ww =1

根据拉格朗日乘子： $S_Bw = \lambda S_Ww$ ，其它同上。

S_Bw 与 $(\mathbf{m_2},\mathbf{m_1})$ 同向，此外，因为不考虑w的长度（分子分母均是w的二次项，因此解与w的长度无关），只考虑其方向，所以可以去除伸缩因子 w^TS_Bw,w^TS_Ww ,然后两边同乘 $S_w^{-1}$ ，则有：

$w\propto S_w^{-1}(\mathbf{m_2}-\mathbf{m_1})$

多个类别：

考虑 K>2 的类别，假设数据维度D要远大于K。 S_W 的计算同上。

$S_W=\sum_{k=1}^KS_k$

$S_k=\sum_{n\in C_k}^N(x_n-m_k)(x_n-m_k)^T$

$m_k = \frac{1}{N_k}\sum_{n \in C_k}x_n$

为了找到 S_B 的一般形式，为此引入总协方差矩阵(total covariance matrix):

$S_T=\sum_{n=1}^N(x_n-m)(x_n-m)^T$

$m=\frac{1}{N}\sum_{n=1}^Nx_n = \frac{1}{N}\sum_{k=1}^KN_km_k$

$N=\sum_k^KN_k$

m:所有数据集的平均值， N_k 每个类别的数量，N数据集的总数量,，K是总类别数量，k指某个类别

S_T = S_W+S_B

$\sum_{n=1}^N(x_n-m)(x_n-m)^T=\sum_{k=1}^KS_k+S_W$

$\sum_{n=1}^N(x_n-m)(x_n-m)^T=\sum_{k=1}^K\sum_{n\in C_k}^N(x_n-m_k)(x_n-m_k)^T+S_W$ $(x_1-m)(x_1-m)^T+...+(x_N-m)(x_N-m)^T=\sum_{k=1}^K\sum_{n\in C_k}^N(x_n-m_k)(x_n-m_k)^T+S_W$

$\sum_{k=1}^K\sum_{n\in C_k}^N(x_n-m)(x_n-m)^T=\sum_{k=1}^K\sum_{n\in C_k}^N(x_n-m_k)(x_n-m_k)^T+S_W$

$S_W=\sum_{k=1}^K\sum_{n\in C_k}^N(x_n-m)(x_n-m)^T-\sum_{k=1}^K\sum_{n\in C_k}^N(x_n-m_k)(x_n-m_k)^T$

$=\sum_{k=1}^K\sum_{n\in C_k}^N(x_n-m-x_n+m_k)(x_n-m-x_n+m_k)^T$

$=\sum_{k=1}^KN_k(m_k-m)(m_k-m)^T$

3、Python3实现

先观察每个特征的分布情况，LDA假设数据是正态分布的，特征是统计独立的，并且每个类都有相同的协方差矩阵。但是，这只适用于LDA作为分类器，如果违背了这些假设，那么降维的LDA也可以很好地工作。

from sklearn import datasets
import matplotlib.pyplot as plt
import numpy as np
import math

# prepare the data
iris = datasets.load_iris()
X = iris.data
y = iris.target
names = iris.feature_names
labels = iris.target_names
y_c = np.unique(y)

"""visualize the distributions of the four different features in 1-dimensional histograms"""
fig, axes = plt.subplots(2, 2, figsize=(12, 6))
for ax, column in zip(axes.ravel(), range(X.shape[1])):
    # set bin sizes
    min_b = math.floor(np.min(X[:, column]))
    max_b = math.ceil(np.max(X[:, column]))
    bins = np.linspace(min_b, max_b, 25)

    # plotting the histograms
    for i, color in zip(y_c, ('blue', 'red', 'green')):
        ax.hist(X[y == i, column], color=color, label='class %s' % labels[i],
                bins=bins, alpha=0.5, )
    ylims = ax.get_ylim()

    # plot annotation
    l = ax.legend(loc='upper right', fancybox=True, fontsize=8)
    l.get_frame().set_alpha(0.5)
    ax.set_ylim([0, max(ylims) + 2])
    ax.set_xlabel(names[column])
    ax.set_title('Iris histogram feature %s' % str(column + 1))

    # hide axis ticks
    ax.tick_params(axis='both', which='both', bottom=False, top=False, left=False, right=False,
                   labelbottom=True, labelleft=True)

    # remove axis spines
    ax.spines['top'].set_visible(False)
    ax.spines["right"].set_visible(False)
    ax.spines["bottom"].set_visible(False)
    ax.spines["left"].set_visible(False)

axes[0][0].set_ylabel('count')
axes[1][0].set_ylabel('count')
fig.tight_layout()
plt.show()

这些简单的特征图形表示，花瓣的长度和宽度可能更适合分类。

Step 1: Computing the d-dimensional mean vectors

np.set_printoptions(precision=4)

mean_vector = []  # 类别的平均值
for i in y_c:
    mean_vector.append(np.mean(X[y == i], axis=0))
    print('Mean Vector class %s:%s\n' % (i, mean_vector[i]))

Step 2: Computing the Scatter Matrices

"""within-class scatter matrix"""
S_W = np.zeros((X.shape[1], X.shape[1]))
for i in y_c:
    Xi = X[y == i] - mean_vector[i]
    S_W += np.mat(Xi).T * np.mat(Xi)
print('within-class scatter matrix:\n', S_W)

"""between-class scatter matrix """
S_B = np.zeros((X.shape[1], X.shape[1]))
mu = np.mean(X, axis=0)  # 所有样本平均值
for i in y_c:
    Ni = len(X[y == i])
    S_B += Ni * np.mat(mean_vector[i] - mu).T * np.mat(mean_vector[i] - mu)
print('within-class scatter matrix:\n', S_B)

Step 3: Solving the generalized eigenvalue problem for the matrix $S_W^{-1} S_B$

eigvals, eigvecs = np.linalg.eig(np.linalg.inv(S_W) * S_B)  # 求特征值，特征向量
np.testing.assert_array_almost_equal(np.mat(np.linalg.inv(S_W) * S_B) * np.mat(eigvecs[:, 0].reshape(4, 1)),
                                     eigvals[0] * np.mat(eigvecs[:, 0].reshape(4, 1)), decimal=6, err_msg='',
                                     verbose=True)

Step 4: Selecting linear discriminants for the new feature subspace

# sorting the eigenvectors by decreasing eigenvalues
eig_pairs = [(np.abs(eigvals[i]), eigvecs[:, i]) for i in range(len(eigvals))]
eig_pairs = sorted(eig_pairs, key=lambda k: k[0], reverse=True)
W = np.hstack((eig_pairs[0][1].reshape(4, 1), eig_pairs[1][1].reshape(4, 1)))

Step 5: Transforming the samples onto the new subspace

X_trans = X.dot(W)
assert X_trans.shape == (150, 2)

对比sklearn

plt.figure(figsize=(8, 4))
plt.subplot(121)
plt.scatter(X_trans[y == 0, 0], X_trans[y == 0, 1], c='r')
plt.scatter(X_trans[y == 1, 0], X_trans[y == 1, 1], c='g')
plt.scatter(X_trans[y == 2, 0], X_trans[y == 2, 1], c='b')
plt.title('my LDA')
plt.xlabel('LD1')
plt.ylabel('LD2')
plt.legend(labels, loc='best', fancybox=True)

from sklearn.discriminant_analysis import LinearDiscriminantAnalysis

X_trans2 = LinearDiscriminantAnalysis(n_components=2).fit_transform(X, y)
plt.subplot(122)
plt.scatter(X_trans2[y == 0, 0], X_trans2[y == 0, 1], c='r')
plt.scatter(X_trans2[y == 1, 0], X_trans2[y == 1, 1], c='g')
plt.scatter(X_trans2[y == 2, 0], X_trans2[y == 2, 1], c='b')
plt.title('sklearn LDA')
plt.xlabel('LD1')
plt.ylabel('LD2')
plt.legend(labels, loc='best', fancybox=True)

我不爱机器学习

关注

16
点赞
踩
110

收藏

觉得还不错? 一键收藏
8
评论
python3 Fisher线性判别分析(LDA)(含详细推导和代码）

1、线性判别原理线性判别分析是常用的降维技术,在模式分类和机器学习的预处理步骤中。其目标是将数据集投影到具有良好的类可分性的低维空间中，以避免过度拟合（维数过多）并降低计算成本，如将一个特征空间(一个数据集n维样本)投射到一个更小的子空间k(其中k ≤n-1)上，同时维护类区分信息。判别式是一个函数，它接受一个输入向量x，并把它赋值给K个类中的一个，记作。在这一章中，我们将把注意力限制在线...
复制链接

扫一扫

专栏目录