线性分类器及Python实现

线性分类实际问题

样本描述

二分类问题诸如下图：

有两类坐标点（图中用红点和蓝点来表示），我们需要寻找一条最优直线来将它们分离开。这些样本点也可能是高维空间中的点，因此需要一个超平面来将两类样本点分离。
这里我们以二维平面上的点集$(x_i^1,x_i^2)$来表示，每个点拥有一个特征，由此将它们分成两类。这个分类函数我们可以如下表示：
$$f(x_i^1,x_i^2)=\{\begin{array}{cc}红色& (x_i^1,x_i^2)满足某条件1\\ 蓝色& (x_i^1,x_i^2)满足某条件2\end{array}$$
在这幅图中可以很清楚地认定红色的点位置靠上，蓝色的点位置靠下。因此我们可以用$"1"$和$"0"$来区分这两类样本点，即：
$$f(x_i^1,x_i^2)=\{\begin{array}{cc}1& i\in D\\ 0& i\notin D\end{array}$$
其中$D$为红色样本点序号所构成的集合，由此完成了用$"1"$和$"0"$来区分这两类样本点，其中$f(x_i^1,x_i^2)=1$表示靠近上方位置的点，$f(x_i^1,x_i^2)=0$表示靠近下方位置的点。
设我们要寻找的直线为$w_0+w_1{x}^1+w_2{x}^2=0$，我们需要保证这两类样本点都尽可能地远离这条直线。由空间位置我们可以得到
$$\{\begin{array}{cc}{w}_0+w_1{x}_i^1+w_2{x}_i^2>0& when:f(x_i^1,x_i^2)=1\\ w_0+w_1{x}_i^1+w_2{x}_i^2<0& when:f(x_i^1,x_i^2)=0\end{array}$$
这里只需要保证$w_2>0$就能保证$w_0+w_1{x}_i^1+w_2{x}_i^2>0$表示的是上半平面部分。

爱因斯坦求和约定

我们引入爱因斯坦求和约定：当一个求和式子角标同时包含上下两个指标时，求和符号可以省略不写，例如记：
$$\sum _{\alpha }{u}_{\alpha }{v}^{\alpha }=u_{\alpha }{v}^{\alpha }$$
当求和式子有多个指标时也同样适用，例如：
$$\sum _{\alpha ,\beta }{u}_{\alpha \beta }{v}_{\gamma }^{\alpha \beta }=u_{\alpha \beta }{v}_{\gamma }^{\alpha \beta }$$

需要注意的是，求和符号$\sum$内部需要至少两个乘积项，而且求和指标必须分别在两个乘积项的上标和下标。记$x_i^0=1$，于是上式我们可以简写为：
$$\{\begin{array}{cc}{w}_j{x}_i^j>0& when:f(x_i^1,x_i^2)=1\\ w_j{x}_i^j<0& when:f(x_i^1,x_i^2)=0\end{array}(w_2>0,j=0,1,2)$$

寻找分类直线（超平面）的方法

根据之前的分析，我们需要保证这两类样本点都尽可能地远离目标直线$w_0+w_1{x}^1+w_2{x}^2=0$，对于满足$f(x_i^1,x_i^2)=1$的点有$w_0+w_1{x}_i^1+w_2{x}_i^2>0$；满足$f(x_i^1,x_i^2)=0$的点有$w_0+w_1{x}_i^1+w_2{x}_i^2<0$
因此$w_0+w_1{x}_i^1+w_2{x}_i^2=w_j{x}_i^j$可以作为一个指标，对于满足$f(x_i^1,x_i^2)=1$的点我们希望$w_j{x}_i^j$越大越好；满足$f(x_i^1,x_i^2)=0$的点我们希望$w_j{x}_i^j$越小越好。
写到这里，读者应该能想到用阶跃函数来描述:

但是如果我们需要用该直线去预测新的样本点应该属于哪部分的概率，单位阶跃函数就失效了，因为单位阶跃函数$u(x)$值域只有$0$和$1$，没有办法表征样本点的置信程度。例如将第一类样本点映射到右半部分，无论在哪个区域所对应的值均为$1$，输出量没有表征分类直线将样本点区分的决策优化程度。

LOGISTIC分析

下面我们提供一个性质良好的函数：逻辑斯蒂（LOGISTIC）函数，记为$g(x)=\frac{1}{1+e^{-x}}$

LOGISTIC函数有非常良好的性质：

• 性质一：$0<g(x)<1$，其值域在$0$到$1$之间，可以表示概率。
• 性质二：$g(0)=0.5$，可以作为上述两类样本点的分界，且该函数以此点为对称中心。在横轴左半轴区域的事件不太可能发生，在横轴右半轴区域的事件很有可能发生。
• 性质三：$g(x)$单调递增，可以作为两类样本点的概率划分指标。
• 性质四：微分性质，$g^{\prime }(x)=g(x)(1-g(x))$，事实上，LOGISTIC函数就是通过该微分方程解得的，它为我们提供了求$g(x)$微分的初等方法，可以大大简化计算机的运算。

通过LOGISTIC函数，我们可以将满足$f(x_i^1,x_i^2)=1$的点映射到输入的正半轴，将满足$f(x_i^1,x_i^2)=0$的点映射到输入的负半轴，中央$x=0$的点就是两类样本点的分界，两类样本点都需尽量远离中心。根据之前的讨论，这个距离指标恰恰可以用$w_j{x}_i^j$来描述。

1. 满足$f(x_i^1,x_i^2)=1$点的映射：
   构造概率函数$p_i=g(w_j{x}_i^j)$，我们需要让样本点$(x_i^1,x_i^2)$尽量远离直线，即$w_j{x}_i^j$尽可能地大，由性质三LOGISTIC函数单增性质，即需要$p_i$尽可能地大。
2. 满足$f(x_i^1,x_i^2)=0$点的映射：
   构造概率函数$p_i=1-g(w_j{x}_i^j)$，我们需要让样本点$(x_i^1,x_i^2)$尽量远离直线，即$w_j{x}_i^j$尽可能地小，于是需要$p_i$尽可能地大。

对于两类样本点，可以将$p_i$合写为
$$p_i=g(w_j{x}_i^j{)}^{f(x_i^1,x_i^2)}\cdot (1-g(w_j{x}_i^j){)}^{1-f(x_i^1,x_i^2)}$$
这个函数即表征分类之依据，我们需要让样本点$(x_i^1,x_i^2)$尽量远离直线，即需要$p_i$尽可能地大。

极大似然估计的梯度上升法

由之前的讨论，我们构造分类函数
$$L(w_0,w_1,w_2)=\prod _i{p}_i=\prod _{i}g(w_j{x}_i^j{)}^{f(x_i^1,x_i^2)}\cdot (1-g(w_j{x}_i^j){)}^{1-f(x_i^1,x_i^2)}$$
取对数得$$l(w_0,w_1,w_2)=lnL(w_0,w_1,w_2)=\sum _i[f(x_i^1,x_i^2)ln(g(w_j{x}_i^j))+(1-f(x_i^1,x_i^2))ln(1-g(w_j{x}_i^j))]$$
于是目标转化为求$l(w_0,w_1,w_2)$的最大值，这里我们采用梯度上升法，
（梯度上升（下降）法可以参考我之前的文章：梯度下降法线性回归模拟）

对$w_0,w_1,w_2$逐个求偏导得：
$$\frac{\partial l}{\partial w_k}=\sum _i(f(x_i^1,x_i^2)-g(w_j{x}_i^j))x_i^k$$
若设$n$为迭代次数，求得迭代公式如下：
$$\{\begin{array}{c}{w}_k^{(n+1)}=w_k^{(n)}+\alpha {\sum }_i(f(x_i^1,x_i^2)-g(w_j^{(n)}{x}_i^j))x_i^k\\ l^{(n)}={\sum }_i[f(x_i^1,x_i^2)ln(g(w_j^{(n)}{x}_i^j))+(1-f(x_i^1,x_i^2))ln(1-g(w_j^{(n)}{x}_i^j))]\end{array}$$
其中$k=0,1,2;x_i^0=1$
写成展开形式为
$$\{\begin{array}{c}{w}_0^{(n+1)}=w_0^{(n)}+\alpha {\sum }_i(f(x_i^1,x_i^2)-\frac{1}{1+e^{-(w_0^{(n)}+w_1^{(n)}{x}_i^1+w_2^{(n)}{x}_i^2)}})\\ w_1^{(n+1)}=w_1^{(n)}+\alpha {\sum }_i{x}_i^1(f(x_i^1,x_i^2)-\frac{1}{1+e^{-(w_0^{(n)}+w_1^{(n)}{x}_i^1+w_2^{(n)}{x}_i^2)}})\\ w_2^{(n+1)}=w_2^{(n)}+\alpha {\sum }_i{x}_i^2(f(x_i^1,x_i^2)-\frac{1}{1+e^{-(w_0^{(n)}+w_1^{(n)}{x}_i^1+w_2^{(n)}{x}_i^2)}})\\ l^{(n)}={\sum }_i[f(x_i^1,x_i^2)ln(\frac{1}{1+e^{-(w_0^{(n)}+w_1^{(n)}{x}_i^1+w_2^{(n)}{x}_i^2)}})+(1-f(x_i^1,x_i^2))ln(1-\frac{1}{1+e^{-(w_0^{(n)}+w_1^{(n)}{x}_i^1+w_2^{(n)}{x}_i^2)}})]\end{array}$$
其中$\alpha$为学习参数，选择合适的$\alpha$，当经过足够多的迭代次数之后，$w_i$收敛到$\widehat{w_i}$，$l(w_0,w_1,w_2)$收敛到$l_M$，$l_M$即为分类函数$l(w_0,w_1,w_2)$的最大值。在这种情况下收敛的$\widehat{w_i}$所满足的直线$\widehat{w_j}{x}_i^j=\widehat{w_0}+\widehat{w_1}{x}^1+\widehat{w_2}{x}^2=0$可以使得两类样本点分别尽可能远离这条直线，因此该直线即为最后分类两类样本点的目标直线。
在初值的选择上$w_i^{(0)}=0$，这样在梯度上升的作用下能保证$w_2^{(k)}>0$，由此满足要求。

样本点正确率的计算和预测

对于已有样本点$(x_i^1,x_i^2)$，我们可以通过比较$\widehat{w_0}+\widehat{w_1}{x}^1+\widehat{w_2}{x}^2$与$0$的关系来实现，也就是计算$p_i=g(\widehat{w_j}{x}_i^j{)}^{f(x_i^1,x_i^2)}\cdot (1-g(\widehat{w_j}{x}_i^j){)}^{1-f(x_i^1,x_i^2)}$是否大于$0.5$来判断是否分类成功。
若给定新的样本点$(x_{ia}^1,x_{ia}^2)$， $p_i$即为表征置信概率的指标。
若把新的样本点$(x_{ia}^1,x_{ia}^2)$划分到红色点集，即令$f(x_{ia}^1,x_{ia}^2)=1$，由此算出的$p_i$即为这种划分模式的可信度。
其中置信函数的LOGISTIC形式为：
$$g(x^1,x^2)=\frac{1}{1+e^{-(\widehat{w_0}+\widehat{w_1}{x}^1+\widehat{w_2}{x}^2)}}$$

代码实现

数据导入

本例中训练集和测试集的数据如下：
百度网盘分享：Test_Train_Set.xls，提取码：mr6d

数据的训练集和测试集

训练集样本$(x_i,y_i)$，用该样本来计算分类直线：

测试集样本$(x_i^{\prime },y_i^{\prime })$，用之前训练出的分类直线来验证用训练集算出的分类直线的正确率：

Python源代码

代码如下：

import xlrd as xd
import numpy as np
from matplotlib import cm
import matplotlib.pyplot as plt
import math
from mpl_toolkits.mplot3d import Axes3D

# matplotlib画图中中文显示会有问题，需要这两行设置默认字体可以显示中文
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# -----------------------------------------------------------------------------------
# # 读取数据
# 打开excel表所在路径

data = xd.open_workbook('Test_Train_Set.xls')
sheet = data.sheet_by_name('Sheet1')
Data_All = []
test_Desired = []
test_Input1 = []
test_Input2 = []
train_Desired = []
train_Input1 = []
train_Input2 = []

# 将表中数据按行逐步添加到列表中，最后转换为list结构
for r in range(sheet.nrows):
    data1 = []
    for c in range(sheet.ncols):
        data1.append(sheet.cell_value(r, c))
    Data_All.append(list(data1))

    test_Desired.append(data1[0])
    test_Input1.append(data1[1])
    test_Input2.append(data1[2])
    train_Desired.append(data1[3])
    train_Input1.append(data1[4])
    train_Input2.append(data1[5])

# 数据筛选，把表头去掉
test_Desired = test_Desired[1:]
test_Input1 = test_Input1[1:]
test_Input2 = test_Input2[1:]
train_Desired = train_Desired[1:]
train_Input1 = train_Input1[1:]
train_Input2 = train_Input2[1:]

# 数据转化为矩阵
test_DesiredM = np.c_[test_Desired[1:]]
test_Input1M = np.c_[test_Input1[1:]]
test_Input2M = np.c_[test_Input2[1:]]
train_DesiredM = np.c_[train_Desired[1:]]
train_Input1M = np.c_[train_Input1[1:]]
train_Input2M = np.c_[train_Input2[1:]]

# -----------------------------------------------------------------------------------
# 定义函数

# Logistic函数: lgstc'(x) = lgstc(x)(1 - lgstc(x))
def lgstc (x):
    return 1 / (1 + math.exp(-x))

# 目标方程：1 + w1x1 + w2x2 = 0

# w0增量
def Del_w0(Input1, Input2, Desired, w0, w1, w2):
    sumDel = 0
    for i in range(len(Input1)):
        lc_p = lgstc(w0 + w1 * Input1[i] + w2 * Input2[i])
        sumDel += Desired[i] - lc_p
    return sumDel

# w1增量
def Del_w1(Input1, Input2, Desired, w0, w1, w2):
    sumDel = 0
    for i in range(len(Input1)):
        lc_p = lgstc(w0 + w1 * Input1[i] + w2 * Input2[i])
        sumDel += Input1[i] * (Desired[i] - lc_p)
    return sumDel

# w2增量
def Del_w2(Input1, Input2, Desired, w0, w1, w2):
    sumDel = 0
    for i in range(len(Input1)):
        lc_p = lgstc(w0 + w1 * Input1[i] + w2 * Input2[i])
        sumDel += Input2[i] * (Desired[i] - lc_p)
    return sumDel

# Loss函数每次迭代后的损失函数
def Loss(Input1, Input2, Desired, w0, w1, w2):
    loss = 0
    for i in range(len(Input1)):
        lc_p = lgstc(w0 + w1 * Input1[i] + w2 * Input2[i])
        loss += Desired[i] * math.log(lc_p) + (1 - Desired[i]) * math.log(1 - lc_p)
    return loss

# 梯度上升法求极值
def GDM_LC(Input1, Input2, Desired, w0 = 0, w1 = 0, w2 = 0, Alpha = 0.1, err = 1e-9):
    w0_list = []
    w1_list = []
    w2_list = []
    loss_list = []
    IterTime = 0

    while 1:
        w0, w1, w2 =\
            w0 + Alpha * Del_w0(Input1, Input2, Desired, w0, w1, w2),  \
            w1 + Alpha * Del_w1(Input1, Input2, Desired, w0, w1, w2),  \
            w2 + Alpha * Del_w2(Input1, Input2, Desired, w0, w1, w2)
        w0_list.append(w0)
        w1_list.append(w1)
        w2_list.append(w2)
        loss_list.append(Loss(Input1, Input2, Desired, w0, w1, w2))
        IterTime = IterTime + 1

        if IterTime > 3:
            if abs(w1_list[IterTime-1] - w1_list[IterTime - 2]) < err:
                if abs(w2_list[IterTime - 1] - w2_list[IterTime - 2]) < err:
                    break

    return w0, w1, w2, w0_list, w1_list, w2_list, loss_list, IterTime

# 正确率验证函数
def bool_rate(Input1, Input2, Desired, w0, w1, w2):

    boolnum = len(Input1)
    boolT = 0

    for i in range(boolnum):
        if Desired[i] == 0:
            if w0 + w1 * Input1[i] + w2 * Input2[i] < 0:
                boolT += 1
        elif Desired[i] == 1:
            if w0 + w1 * Input1[i] + w2 * Input2[i] >= 0:
                boolT += 1

    return boolT / boolnum

# -----------------------------------------------------------------------------------
#
Alpha = 0.1         # 学习参数
err = 1e-9          # 误差

w0, w1, w2, w0_list, w1_list, w2_list, loss_list, IterTime =\
    GDM_LC(test_Input1, test_Input2, test_Desired, 0, 0, 0, Alpha, err)

print('测试集样本数量：'+str(len(test_Desired)))
print('训练集样本数量：'+str(len(train_Desired)))
print('迭代次数为：'+str(IterTime))
print('学习参数为：'+str(Alpha))
print('误差不超过：'+str(err))
str_line = '线性分类的直线为：'+str(round(w0, 3))+'+'+str(round(w1, 3))+'x+'+str(round(w2, 3))+'y=0'
str_test = '训练集的正确率：'+str(round(100 * bool_rate(test_Input1, test_Input2, test_Desired, w0, w1, w2), 3))+'%'
str_train = '测试集的正确率：'+str(round(100 * bool_rate(train_Input1, train_Input2, train_Desired, w0, w1, w2), 3))+'%'
print(str_line)
print(str_test)
print(str_train)

X = test_Input1M
Y = - w0 / w2 - w1 * X / w2

# -----------------------------------------------------------------------------------
# 绘图
plt.figure('线性分类求解（训练集）')
plt.plot(X, Y, label=str_line)
for i in range(len(test_Desired)):
    if test_Desired[i] == 1:
        plt.scatter(test_Input1[i], test_Input2[i], c='r', marker='+')
    else:
        plt.scatter(test_Input1[i], test_Input2[i], c='b', marker='.')
plt.title('训练集样本')
plt.xlabel('训练集坐标x')
plt.ylabel('训练集坐标y')
plt.legend()

plt.figure('线性分类验证（测试集）')
plt.plot(X, Y, label=str_line)
for i in range(len(train_Desired)):
    if train_Desired[i] == 1:
        plt.scatter(train_Input1[i], train_Input2[i], c='r', marker='+')
    else:
        plt.scatter(train_Input1[i], train_Input2[i], c='b', marker='.')
plt.title('测试集样本')
plt.xlabel('测试集坐标x')
plt.ylabel('测试集坐标y')
plt.legend()

plt.figure('指标收敛图像')

plt.subplot(1, 2, 1)
plt.title('概率函数收敛曲线')
plt.semilogx(loss_list, label='指标函数收敛到'+str(round(loss_list[len(loss_list) - 1], 3)))
plt.legend()

plt.subplot(1, 2, 2)
plt.title('直线的三个指标收敛曲线')
plt.semilogx(w0_list, label='w0收敛到'+str(round(w0_list[len(w0_list) - 1], 3)))
plt.semilogx(w1_list, label='w1收敛到'+str(round(w1_list[len(w1_list) - 1], 3)))
plt.semilogx(w2_list, label='w2收敛到'+str(round(w2_list[len(w2_list) - 1], 3)))
plt.legend()

plt.show()

# 概率预测曲面图(三维图)
fig = plt.figure('概率预测曲面图')
ax = fig.add_subplot(111, projection='3d')

x_ax = np.arange(-1, 1, 0.005)
y_ax = np.arange(0, 1, 0.005)
x_ax, y_ax = np.meshgrid(x_ax, y_ax)
z_ax = 1 / (1 + np.exp(- (w0 + w1 * x_ax + w2 * y_ax)))

ax.plot_surface(x_ax, y_ax, z_ax, rstride=1, cstride=1, cmap=cm.coolwarm)

ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.set_title('概率预测曲面$z=(1+e^{-8.937+2.590x+20.537y})^{-1}$')

x_line = np.linspace(-1, 1, 10)
y_line = - w0 / w2 - w1 * x_line / w2
z_line = 1 / (1 + np.exp(w0 + w1 * x_line + w2 * y_line))
ax.plot(x_line, y_line, z_line, c='k', label='在曲面上的分类直线：'+str(round(w0, 3))+'+'+str(round(w1, 3))+'x+'+str(round(w2, 3))+'y=0')
ax.legend()
plt.show()

输出结果为：

• 测试集样本数量：$30$
• 训练集样本数量：$30$
• 迭代次数为：$13549$
• 学习参数为：$0.1$
• 误差不超过：$1\times 10^{-9}$
• 线性分类的直线为：$-8.937+2.590x+20.537y=0$
• 训练集的正确率：$86.667\%$
• 测试集的正确率：$73.333\%$

$w_0,w_1,w_2,l$的收敛曲线如下：

于是预测样本点的概率函数图像如下：
$$z=\frac{1}{1+e^{8.937-2.590x-20.537y}}$$

例如，随意增添样本点$(0.6,0.6)$，我们希望将该点认定为从点集$D$中取出，即把该点分到$f(0.6,0.6)=1$的情况，可信度为
$$p=\frac{1}{1+e^{8.937-2.590\times 0.6-20.537\times 0.6}}=99.29\%$$

注意问题

直线的一般方程为$w_0+w_{1}x+w_{2}y=0$，有三个未知量，理论上有无数组$(w_0,w_1,w_2)$可以表示同一条直线。但是如果在本例中我们在求解分类直线时只设置两个未知量$1+w_{1}x+w_{2}y=0$，再利用上述方法求解得不到正确的收敛值。因为在LOGISTIC函数中还应在设置一个参数$a$，$$g(a;x)=\frac{1}{1+e^{-ax}}$$
用一般方程$w_0+w_{1}x+w_{2}y=0$去拟合直线只不过把$a$含盖在了里面。
$a$参数也是表征样本性质的参数，比较以下两个样本

很显然样本1的两类点集的分离情况要好于样本2，因此在寻求最优分类直线的时候，样本1两类点集中间的空间有更多的置信空间，再求得分类的最优直线后再给定样本点，其所在区域分类到某一侧的确定度非常高。相比之下样本2的两类点集有交叉，因此样本2的最优分类直线两侧的区域置信概率有一定的模糊性，在选择分类标准时不容易区分。
因此我们可以得到LOGISTIC函数在样本1和样本2分别求出的$a$参数$a_1$和$a_2$有如下关系：
$$a_1>a_2$$