机器学习二元分类算法 二元归类

分类任务一直都是机器学习的基础任务，已经被广泛应用在新闻分类、情感分类、主题分类、图片分类、视频分类等领域。

机器学习分类通过训练集进行学习，建立一个从输入空间 X 到输出空间 Y（离散值）的映射。按输出类别（标签）不同，可以分为二元分类（Binary Classification）、多元分类（Multi-Class Classification）。本文以二元分类为例，介绍一下机器学习在分类问题中的应用。

1 分类问题的模型表达

分类问题的模型表达与回归问题基本相同，区别仅在于输出值 Y 为离散值。通常，我们使用$x^{(i)}$表示输入变量；$y^{(i)}$表示输出或目标变量；$(x^{(i)}, y^{(i)})$表示训练集。我们还将使用 X 表示输入值的值域，使用 Y 表示输出值的类别。

监督学习就是给定一个训练集，学习一个函数$h：X→Y$，这样$h(x)$是$y$的相应值的"良好"预测因子。由于历史原因，$h(x)$被称为假设函数。

2 假设函数

在分类问题中，要预测的变量 $y$ 是离散的值。分类问题的例子有：判断一封电子邮件是否是垃圾邮件；判断一次金融交易是否是欺诈；区别一个肿瘤是恶性的还是良性的。这些分类也都被称为二元分类。一般来说，我们将因变量 (dependent variable) 可能属于的两个类分别称为负向类（negative class）和正向类 (positive class) ，则因变量$y\in { 0,1 \\}$

线性回归模型只能预测连续的值，然而对于分类问题，只需要输出0或1。解决这个问题的一种方法是使用 logistic sigmoid 函数将线性函数的输出压缩进区间 (0, 1) 。该值可以解释为概率：

$p(y=1|x; \theta) = \sigma (\theta ^ T x)$

这个方法被称为逻辑回归 (Logistic Regression) 。

逻辑回归模型的假设函数如下，其作用是，对于给定的输入变量，根据选择的参数计算输出变量等于1的可能性：

$h_\theta \left( x \right)=g\left(\theta^{T}X \right) = p(y=1|x; \theta)$

其中，$X$ 代表特征向量，$g$ 代表逻辑函数（logistic function)。logistic function是一个常用的逻辑函数，也被称为Sigmoid function，公式为：

$g\left( z \right)=\frac{1}{1+{{e}^{-z}}}$

例如，如果对于给定的$x$，通过已经确定的参数计算得出$h_\theta \left( x \right)=0.7$，则表示有70%的几率$y$为正向类，相应地$y$为负向类的几率为1-0.7=0.3。

在逻辑回归中：

• 当${h_\theta}\left( x \right)>=0.5$时，预测 $y=1$。
• 当${h_\theta}\left( x \right)<0.5$时，预测 $y=0$

又 $z={\theta^{T}}x$

• ${\theta^{T}}x>=0$ 时，预测 $y=1$
• ${\theta^{T}}x<0$ 时，预测 $y=0$

假设函数和Sigmoid函数的代码表示如下：

import numpy as np

def sigmoid(z):
    return 1 / (1 + np.exp(-z))

def hypothesis_fun(x, theta):
    return sigmoid(theta @ x)

Sigmoid函数的图像如下：

import matplotlib.pyplot as plt

x = np.linspace(-5, 5, 100)
plt.ylim(0, 1)
plt.vlines(x=0, ymin=0, ymax=1, colors='k')
plt.plot(x, sigmoid(x));

3 损失函数

对于线性回归模型，我们定义的损失函数是所有模型误差的平方和。理论上来说，我们也可以对逻辑回归模型沿用这个定义。但是问题在于，当我们将二元分类的假设函数${h_\theta}\left( x \right)=\frac{1}{1+{e^{-\theta^{T}x}}}$带入到均方误差函数时，我们得到的损失函数将是一个非凸函数（non-convexfunction）。

这意味着我们的损失函数有许多局部最小值，这将影响梯度下降算法寻找全局最小值。

因此，我们重新定义逻辑回归的损失函数为：

$J\left( \theta \right)=\frac{1}{m}\sum\limits_{i=1}^{m}{{Cost}\left( {h_\theta}\left( {x}^{\left( i \right)} \right),{y}^{\left( i \right)} \right)}$

其中,

$Cost(h_\theta (x), y) = \begin{cases} -\log{(h_\theta(x))} \quad & y=1 \\ -\log{(1 - h_\theta(x))} \quad & y=0 \end{cases}$

这样构建的$Cost\left( {h_\theta}\left( x \right),y \right)$函数的特点是：当实际的 $y=1$ 且${h_\theta}\left( x \right)$也为 1 时误差为 0，当 $y=1$ 但${h_\theta}\left( x \right)$不为1时误差随着${h_\theta}\left( x \right)$变小而变大；当实际的 $y=0$ 且${h_\theta}\left( x \right)$也为 0 时代价为 0，当$y=0$ 但${h_\theta}\left( x \right)$不为0时误差随着 ${h_\theta}\left( x \right)$的变大而变大。

进而将构建的 $Cost\left( {h_\theta}\left( x \right),y \right)$简化如下：

$Cost\left( {h_\theta}\left( x \right),y \right)=-y\times log\left( {h_\theta}\left( x \right) \right)-(1-y)\times log\left( 1-{h_\theta}\left( x \right) \right)$

带入损失函数得到：

$J\left( \theta \right)=-\frac{1}{m}\sum\limits_{i=1}^{m}{[{{y}^{(i)}}\log \left( {h_\theta}\left( {{x}^{(i)}} \right) \right)+\left( 1-{{y}^{(i)}} \right)\log \left( 1-{h_\theta}\left( {{x}^{(i)}} \right) \right)]}$

这时，损失函数$J(\theta)$会是一个凸函数，并且没有局部最优值。损失函数的代码表示如下：

def loss_fun(x, y, theta):
    y_p = hypothesis_fun(x, theta)
    first = y * np.log(y_p)
    second = (1 - y) * np.log(1 - y_p)
    return -np.sum(first + second)/(y.size)

4 梯度优化

在得到这样一个损失函数以后，我们便可以用梯度下降算法来求得能使损失函数最小的参数了。算法为：

$\theta_j := \theta_j - \alpha \frac{\partial}{\partial\theta_j} J(\theta) \quad \quad j=0,1,2,...n$

求导后得到：

$\theta_j := \theta_j - \alpha \frac{1}{m}\sum\limits_{i=1}^{m}{{\left( {h_\theta}\left( \mathop{x}^{\left( i \right)} \right)-\mathop{y}^{\left( i \right)} \right)}}\mathop{x}_{j}^{(i)} \quad \quad j=0,1,2,...n$

具体的推导过程如下，考虑到：

${h_\theta}\left( {{x}^{(i)}} \right)=\frac{1}{1+{{e}^{-{\theta^T}{{x}^{(i)}}}}}$

则：

$\begin{aligned} J\left( \theta \right)&= -\frac{1}{m}\sum\limits_{i=1}^{m}{[{{y}^{(i)}}\log \left( {h_\theta}\left( {{x}^{(i)}} \right) \right)+\left( 1-{{y}^{(i)}} \right)\log \left( 1-{h_\theta}\left( {{x}^{(i)}} \right) \right)]} \newline &=-\frac{1}{m}\sum\limits_{i=1}^{m}{[{{y}^{(i)}}\log \left( \frac{1}{1+{{e}^{-{\theta^T}{{x}^{(i)}}}}} \right)+\left( 1-{{y}^{(i)}} \right)\log \left( 1-\frac{1}{1+{{e}^{-{\theta^T}{{x}^{(i)}}}}} \right)]} \newline &=-\frac{1}{m}\sum\limits_{i=1}^{m}{[-{{y}^{(i)}}\log \left( 1+{{e}^{-{\theta^T}{{x}^{(i)}}}} \right)-\left( 1-{{y}^{(i)}} \right)\log \left( 1+{{e}^{{\theta^T}{{x}^{(i)}}}} \right)]} \end{aligned}$

所以：

$$
\begin{aligned}
\frac{\partial }{\partial {\theta_{j}}}J\left( \theta \right) &=\frac{\partial }{\partial {\theta_{j}}}[-\frac{1}{m}\sum\limits_{i=1}^{m}{[-{{y}{(i)}}\log \left( 1+{{e}^{-{\theta{T}}{{x}^{(i)}}}} \right)-\left( 1-{{y}^{(i)}} \right)\log \left( 1+{{e}^{{\theta{T}}{{x}^{(i)}}}} \right)]}]

\newline &=-\frac{1}{m}\sum\limits_{i=1}^{m}{[-{{y}{(i)}}\frac{-x_{j}^{(i)}{{e}{-{\theta^{T}}{{x}{(i)}}}}}{1+{{e}^{-{\theta{T}}{{x}^{(i)}}}}}-\left( 1-{{y}^{(i)}} \right)\frac{x_j^{(i)}{{e}{{\theta^T}{{x}{(i)}}}}}{1+{{e}^{{\thetaT}{{x}^{(i)}}}}}}]

\newline &=-\frac{1}{m}\sum\limits_{i=1}^{m}{[{y}{(i)}}\frac{x_j^{{(i)}}{1+{{e}}{{\theta^T}{{x}{(i)}}}}}-\left( 1-{{y}^{(i)}} \right)\frac{x_j^{(i)}{{e}{{\theta^T}{{x}{(i)}}}}}{1+{{e}^{{\thetaT}{{x}^{(i)}}}}}]

\newline &=-\frac{1}{m}\sum\limits_{i=1}^{{m}{\frac{{{y}}{(i)}}x_j^{(i)}-x_j{(i)}{{e}^{{\thetaT}{{x}^{{(i)}}}}+{{y}}{(i)}}x_j^{(i)}{{e}{{\theta^T}{{x}{(i)}}}}}{1+{{e}^{{\thetaT}{{x}^{(i)}}}}}}

\newline &=-\frac{1}{m}\sum\limits_{i=1}^{{m}{\frac{{{y}}{(i)}}\left( 1\text{+}{{e}^{{\thetaT}{{x}^{(i)}}}} \right)-{{e}^{{\thetaT}{{x}^{{(i)}}}}}{1+{{e}}{{\theta^T}{{x}{(i)}}}}}x_j^{(i)}}

\newline &=-\frac{1}{m}\sum\limits_{i=1}^{m}{({{y}{(i)}}-\frac{{{e}^{{\thetaT}{{x}^{{(i)}}}}}{1+{{e}}{{\theta^T}{{x}{(i)}}}}})x_j^{(i)}}

\newline &=-\frac{1}{m}\sum\limits_{i=1}^{m}{({{y}{(i)}}-\frac{1}{1+{{e}^{-{\thetaT}{{x}^{{(i)}}}}})x_j}{(i)}}

\newline &=-\frac{1}{m}\sum\limits_{i=1}^{m}{[{{y}{(i)}}-{h_\theta}\left( {{x}^{(i)}} \right)]x_j^{(i)}}

\newline &=\frac{1}{m}\sum\limits_{i=1}^{m}{[{h_\theta}\left( {{x}^{(i)}} \right)-{{y}^{(i)}}]x_j{(i)}}
\end{aligned}
$$

线性回归中，我们能够通过求解正规方程以找到最佳权重。相比而言，逻辑回归会更困难些。其最佳权重没有闭解，只能通过梯度下降算法最小化损失函数来搜索。具体代码如下：

def gradient_decent(x, y, theta, learning_rate):
    gradient = x @ (hypothesis_fun(x, theta) - y)/ y.size
    return theta - learning_rate * gradient

5 测试数据集

关注公众号，回复20230310下载测试数据。测试数据共三列，前两列为自变量，最后一列为因变量。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

path = 'class.data'
data = pd.read_csv(path)
data.head()

	X	Y	Admitted
0	34.623660	78.024693	0
1	30.286711	43.894998	0
2	35.847409	72.902198	0
3	60.182599	86.308552	1
4	79.032736	75.344376	1

x = data[['X', 'Y']].to_numpy()
x = x.transpose()
y = data['Admitted'].to_numpy()
x.shape, y.shape

((2, 100), (100,))

positive = data[data['Admitted'].isin([1])]
negative = data[data['Admitted'].isin([0])]

fig, ax = plt.subplots(figsize=(12, 8))
ax.scatter(positive['X'], positive['Y'], s=50, 
           c='b', marker='o', label='Admitted')
ax.scatter(negative['X'], negative['Y'], s=50, 
           c='r', marker='x', label='Not Admitted')
ax.legend()
ax.set_xlabel('Exam 1 Score')
ax.set_ylabel('Exam 2 Score')
plt.show()

6 模型训练与预测

二元分类和多元线性回归相同，在正式训练之前需要对输入变量进行了无量纲化。此外，在输入变量x矩阵的第一行前插入全是1的一行，用于表示偏置项。

x2 = feature_scale(x)
x2 = np.vstack((np.ones(x2.shape[1]), x2))

parameters = np.random.rand(x2.shape[0])
learning_rate = 0.1
losses = []

batch_size = 70
epoch_size = 1000

for epoch in range(epoch_size):
    for i in range(y.size//batch_size+1):
        random_samples = np.random.choice(x2.shape[1], batch_size)
        parameters = gradient_decent(x2[:, random_samples], 
                              y[random_samples], 
                              parameters, 
                              learning_rate)
        losses.append(loss_fun(x2, y, parameters))

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 9))
ax1.plot(losses)
ax1.set_xlabel('Iteration number')
ax1.set_ylabel('Loss')

predict_y = hypothesis_fun(x2, parameters)
predict_y[predict_y>=0.5] = 1
predict_y[predict_y<1] = 0
ax2.scatter(x[0][predict_y!=y], x[1][predict_y!=y], s=80, 
            c='r', marker='o', label='Wrong', alpha=1)
ax2.scatter(x[0][predict_y==1], x[1][predict_y==1], s=50, 
            c='b', marker='*', label='Admitted')
ax2.scatter(x[0][predict_y==0], x[1][predict_y==0], s=50, 
            c='b', marker='x', label='Not Admitted')
ax2.legend()
ax2.set_xlabel('X')
ax2.set_ylabel('Y')
plt.show();

此外，和非线性回归的原理类似，在进行二元分类时也可以对输入变量进行变换，以使假设函数更好的逼近预测值。而这种变换同样不需要对模型进行任何改变，只需要改变一下输入即可。如，将输入变为如下形式：

times = 3
x2 = x.copy()
for i in range(2, times+1):
    x2 = np.vstack((x2, x2**i))
x2 = feature_scale(x2)
x2 = np.vstack((np.ones(x2.shape[1]), x2))

此时预测结果如下，可见预测结果更为准确。但需要注意的是，变换之后由于自变量数量增加需要更多的训练次数（epoch）：