《Machine Learning in Action》—— Taoye给你讲讲Logistic回归是咋回事

2020 年 12 月 7 日
笔记
numpy, Python, 数据挖掘, 机器学习, 算法

在手撕机器学习系列文章的上一篇，我们详细讲解了线性回归的问题，并且最后通过梯度下降算法拟合了一条直线，从而使得这条直线尽可能的切合数据样本集，已到达模型损失值最小的目的。

在本篇文章中，我们主要是手撕Logistic回归，这个在李航老师的《统计学习方法》一书中也叫做为逻辑斯谛回归。听到回归一词，有的读者可能会想，上一篇线性回归求解的是拟合问题，这篇文章手撕的是Logistic回归，会不会也是一个拟合问题？只不过使用到的算法原理不同而已，而求解的问题是一致的？？？

其实不然，Logistic回归模型是一种广义的线性回归模型，主要是针对的是分类问题。只不过这个分类模型与上一章中的拟合模型有一些相似之处，或者我们可以这样说，如果你明白了上篇文章中所讲解的线性回归，那么这篇文章中所涉及到的Logistic回归的内容，你将会学的很轻松。这也是为什么Taoye首先肝线性回归的原因。

手撕机器学习系列目前已经更新了八篇，读者可自行根据需要“充电”（持续更新中）：

本文主要包括以下两个部分的内容：

和Logistic回归交个朋友，互相认识一下吧（原理剖析，公式推导）
基于Logistic回归手撕解决二分类问题

一、和Logistic回归交个朋友，互相认识一下吧（原理剖析）

关于机器学习中的分类算法，我们前面已经肝了不少。而Logistic回归是同样也是分类算法中的一员，一般用于二分类问题，比如患者是否患有胃癌、明天是否会下雨等等。当然了，对于多分类问题，Logistic回归也是有应对之法的，毕竟你有张良计，我也有过墙体嘛。本文主要是以二分类问题为例，来剖析Logistic回归中的那些小秘密。

假设现在有一些样本数据点，每个样本包括两个属性，我们可以利用一条直线来对其进行拟合，这个拟合的过程就叫做回归，拟合效果如下所示：

这个也是我们在上一篇文章中所详细讲到的内容，具体可见：《Machine Learning in Action》—— 浅谈线性回归的那些事： //www.zybuluo.com/tianxingjian/note/1761762

而Logistic回归是一种分类算法，它分类的核心思想是建立在线性回归基础之上的，并且对其进行了拓展，主要利用到了Sigmoid函数阈值在[0,1]这么一个特性，而这个特性则正好符合概率所在区间。所以说，Logistic回归的本质就是一个基于条件概率的判别模型（在已知样本属性特征的前提下，判断其属于哪一类别）。

现在就让我们来和Logistic回归互相认识一下吧，看看能不能挖出它的一些小秘密，一睹庐山真面目。

不知道各位看官还记不记得之前在讲解SVM的时候，我们是通过一种间隔最大化的形式来找出最佳的决策面，以此来对数据集进行分类。

而在Logistic回归看来，每一个样本的类别标签都对应着一个概率，哪个类别标签的概率大，那么我就将样本归于哪一类。我们不妨假设单个样本有\(N+1\)个属性特征，我们可以把它看做一个向量形式，且每个样本对应的标签情况只可能有两类，即：

\[\begin{aligned}
& x_i = (x_i^{(0)},x_i^{(1)},x_i^{(2)},…,x_i^{(N)}) \\
&y_i\in (0,1)
\end{aligned}
\]

也就是说，我们现在的目的就是找到一个超平面，在超平面的一侧为一个类别样本，在另一侧为另一类样本。且每个样本只可能存在两种情况，非0即1。对此，我们不妨假设在已知样本属性特征和模型参数\(w=(w_0,w_1,…,w_N)\)的前提下，令该样本标签类别为1的概率为\(h_w(x)\)，因为只可能存在两个类别，所以样本标签为0的概率为\(1-h_w(x)\)，即：

\[\begin{aligned}
& P(y_i=1|x_i;w)=h_w(x_i) \\
& P(y_i=0|x_i;w)=1-h_w(x_i)
\end{aligned}
\]

对于上式，它表示的意思是在已知样本属性特征x和模型参数\(w\)的前提下，该样本标签为1和0的概率分别为\(h_w(x)\)和\(1-h_w(x)\)。根据我们自己意愿，我们当然是希望这两者相差越大越好咯，这样得到的分类结果才更具有说服力。

比如说对于一个样本来讲，该样本标签为0的概率为0.9，为1的概率则等于0.1，我们当然是更加情愿将这个样本归于0那一类别，而且这样的分类结果比较容易服众。那假设该样本标签为0的概率为0.51，为1的概率为0.49，那么这个时候你愿意将这个样本归于哪一类呢？？？是不是很难做出抉择？？？是不是相当的纠结？？？因为这样得出的分类概率结果和我们盲猜的概率也相差不多，完全达不到一个服众的目的。

对此，我们更加希望是这样的：标签的分类概率相差越大越好，这样我们对样本的分类结果更加具有说服力。

上述两个概率其实就是两种情况，且其样本标签为非0即1，根据这个特性，为了方便我们表示分类的概率情况，可以将上述两个概率值合二为一，得到如下：

\[P(y_i|x_i;w)=h_w(x_i)^{y_i}(1-h_w(x_i))^{1-y_i}
\]

上述这种将二者合二为一的处理方式，曾经在SVM那篇文章中也是有谋面过的，如有遗忘，读者可暂且跳转进行复习：《Machine Learning in Action》—— 剖析支持向量机，单手狂撕线性SVM： //www.zybuluo.com/tianxingjian/note/1755051

合并出来的上式，当y等于1时，(1-y)项(第二项)为0；当y等于0时，y项(第一项)为0。也就是说，我们现在的目的就是根据训练样本集，使得上式的值最大化，这样我们对样本集分类的准确性就越高。

此外，上述\(P(y|x_i;w)\)单个样本所表示的概率值，而我们知道，样本集是由大量的单个样本所组成的一个集合，所以此时我们的极大似然法就该出场了，假设样本与样本之间是相互独立的，那么其似然函数为（假定有n个样本）：

\[\prod_{i=1}^nh_w(x_i)^{y_i}(1-h_w(x_i))^{1-y_i}
\]

此时，我们的目的就是要根据数据样本集，来求解上述似然函数的最大值。求最值问题，我们自然会使用到求解导函数，而对于上述的乘积形式，我们直接求导并不简单，而且还会大大提高求解的效率和复杂度。对此我们需要对其进行对数化，从而将乘积的形式转化成求和的形似，这样一来对于求导来讲就比较的友好了。

假设我们对上述式子对数化之后，命名为\(L(w)\)，则其为我们最终的损失函数，或者说是待优化的目标函数，具体形式如下：

\[L(w)=\sum_{i=1}^ny_ilog(h_w(x))+(1-y_i)log(1-h_w(x))
\]

写到这里，我们算是得到了最终所需要求解的目标函数了，现在就是需要将上述的损失函数值最大化，这样才能使得最终分类的结果集的准确性更高。

注意一点：对于上述的损失函数，读者可能在这里有个疑问，按道理将应该是使得损失值最小化才对，为什么会需要使得上式最大化呢？？？其实是这样的，上述的式子说是说损失函数，其实它真正代表的意思是使得整个样本集的分类正确率尽可能的提高。这一点若读者还有疑问，可返回仔细想想该式子的推导过程，及其背后隐藏的真正含义。

这样一来，对于上述的损失函数还有一个\(h_w(x)\)我们是不知道的。通过前面分析，我们也可以知道\(h_w(x)\)表示的是一个概率，其值为0-1。然而，我们通过计算\(w_0x_i^{(0)}+w_1x_i^{(1)}+w_2x_i^{(2)}+…+w_Nx_i^{(N)}\)，可以发现该值的具体范围是不确定的。为此，我们需要对该计算得到的值进行一定的处理，将其转化到0-1范围之间，如此一来，才能符合概率的范围特性。

那么，该如何处理呢？？？

聪明的研究人员就发现，有这么一个函数，无论该值的有多大，或者有多小，都能将其映射到0-1之间。这个函数就是大名鼎鼎的Sigmoid函数，其具体形式和图像如下所示：

\[g(z)=\frac{1}{1+e^{-z}}
\]

通过Sigmoid函数的表达式和具体图像，我们也可以发现，其正好满足我们的实际需求。另外，Sigmoid函数在今后的学习过程中还会经常见到的，比如在卷积神经网络中，就经常会使用Sigmoid函数来作为我们的激活函数，以解决非线性问题。

如此一来，我们通过Sigmoid函数来处理样本，得到如下结果：

\[h_w(x)=g(w^Tx)=\frac{1}{1+e^{-w^Tx}}
\]

了解了\(h_w(x)\)具体表达之后，我们就能对损失函数进行进一步的处理变换，处理过程如下：

\[\begin{aligned}
L(w) & =\sum_{i=1}^ny_ilog(h_w(x_i))+(1-y_i)log(1-h_w(x_i)) \\
& = \sum_{i=1}^n[y_i\log\frac{h_w(x_i)}{1-h_w(x_i)}+log(1-h_w(x_i))] \\
& = \sum_{i=1}^n[y_i(w^Tx_i)-w^Tx_i-log(1+e^{-w^Tx_i})]
\end{aligned}
\]

将损失函数处理成如上式子之后，我们可以发现在整个训练数据样本集中，\(x_i和y_i\)都是已知的，唯一不确定的就是\(w\)，其是一个与单个样本属性特征相对应的向量形式。我们对\(w\)进行求导之后可以得到如下结果：

\[\frac{\partial L(w)}{\partial w}=\sum_{i=1}^n(y_i-h_w(x_i))x_i
\]

由此，我们计算得到损失函数关于\(w\)的梯度之后，就可以通过梯度上升算法来不断更新迭代\(w\)参数，从而使得损失函数的值最大化，即使得训练样本集分类的正确性尽可能高。

另外，我们知道\(w\)参数其实是一个向量的形式，其与\(x\)的属性特征是相对应的，对此我们对\(w\)进行更新的时候，是对其内部的每一个元素同时更新。根据上述的求导结果，我们得到每个元素的具体更新如下：

\[\begin{aligned}
& w_0^{new}=w_0^{old}+\alpha\sum_{i=1}^n(y_i-h_w(x_i))x_i^0 \\
& w_1^{new}=w_1^{old}+\alpha\sum_{i=1}^n(y_i-h_w(x_i))x_i^1 \\
& w_2^{new}=w_2^{old}+\alpha\sum_{i=1}^n(y_i-h_w(x_i))x_i^2 \\
& \vdots \\
& w_N^{new}=w_N^{old}+\alpha\sum_{i=1}^n(y_i-h_w(x_i))x_N^2 \\
\end{aligned}
\]

注意，在这里我们是假设每个元素有\(N+1\)个属性特征的，即\(x_i = (x_i^{(0)},x_i^{(1)},x_i^{(2)},…,x_i^{(N)})\)，这一点在前面也是有提到。其中\(old\)表示的是更新迭代之前的\(w\)值，而\(new\)表示的是更新迭代之后的\(w\)值，而\(\alpha\)表示的是一个学习率，也代表着学习的快慢，这个在之前讲解线性回归的时候也是有详细讲到的。

二、基于Logistic回归手撕解决二分类问题

得到Logistic回归模型之后，我们就能根据\(w^{old}\)来不断更新迭代得到\(w^{new}\)，最终使得损失函数值最大化。

接下来，我们尝试着通过Python代码来实现Logistic回归分类，本次主要针对于二分类。数据集依然通过NumPy随机生成，定义一个establish_data方法随机生成数据集的代码如下：

"""
    Author: Taoye
    微信公众号: 玩世不恭的Coder
    Explain: 通过NumPy准备数据集
    Return:
        x_data：样本的属性特征
        y_label：样本属性特征所对应的标签
"""
def establish_data():
    # np.random.seed(1) # 可根据需要是否选择随机种子
    x_data = np.concatenate((np.add(np.random.randn(50, 2), [1.5, 1.5]),       
                                 np.subtract(np.random.randn(50, 2), [1.5, 1.5])),
                                 axis = 0)      # random随机生成数据，+ -1.5达到不同类别数据分隔的目的 
    y_label = np.concatenate((np.zeros([50]), np.ones([50])), axis = 0)    # concatenate合并数据集
    return x_data. y_label

可视化数据之后的分布如下：

从上图我们可以看出数据的大致分布情况，并且能够通过一条直线将两类数据分割开来。在这里，我们假设Sigmoid函数的输入记为\(z\)，那么\(z=w_0x_0+w_1x_1+w_2x_2\)，即可将数据分割开。其中，为了体现出直线的截距，我们将\(x_0\)当作是数值为1的固定值，\(x_1\)为数据集的第一个属性特征，\(x_2\)为数据集的第二个属性特征。另z=0，则得到该直线的一般表达式\(w_0+w_1x_1+w_2x_2=0\)。

对于这个方程，我们已知的是样本数据，也就是横坐标为\(x_1\)，纵坐标为\(x_2\)，表示样本的两个属性。而未知的参数为\(w_0，w_1，w_2\)，也就是我们需要求的回归系数(最优参数)，也正是需要通过梯度上升算法训练的模型参数。

在开始训练模型参数之前，我们把迭代更新的\(w\)式子再次搬出来看看：

\[\frac{\partial L(w)}{\partial w}=\sum_{i=1}^n(y_i-h_w(x_i))x_i
\]

即：

我们知道，对于求和的式子我们可以将其转换成矩阵或向量的形式来表示，这也是矢量化的一种操作，比如\(\sum_{i=1}^nw_ix_1\)可以转化成\(w^Tx\)。因此，对于上式，我们同样可以对其进行矢量化，得到如下结果：

\[w^{new}=w^{old}+\alpha x^T(y-h_w(x))
\]

这里对上式再解释一下，比如说我们有100个数据样本，每个样本含有2个属性，那么此时的x代表整个样本集，其\(x.shape=(100,2)\)，所以\(x^T.shape=(2,100)\)，而\((y-h_w(x)).shape=(100, 1)\)，所以两者相乘之后得到的shape是\((2,1)\)，正好与所需要的\(w\)向量的维数一致。

根据上述矢量化之后的结果，我们定义一个gradient_ascent方法来通过代码实现这个功能。

"""
    Author: Taoye
    微信公众号: 玩世不恭的Coder
    Explain: Sigmoid函数
    Parameters:
        in_data: sigmoid处理的输入数据
    Return:
        sigmoid_result：sigmoid函数处理之后的结果
"""
def sigmoid(in_data):
    return 1 / (1 + np.exp(-in_data))

"""
    Author: Taoye
    微信公众号: 玩世不恭的Coder
    Explain: Logistic回归的核心方法，主要使用到了梯度上升算法
    Parameters:
        x_data：样本集的属性特征
        y_label：样本集的标签
    Return:
        weights: 模型训练完成之后的w参数
"""
def gradient_ascent(x_data, y_label):
    x_data, y_label = np.mat(x_data), np.mat(y_label).T    # 转换成矩阵形式，方便操作
    data_number, attr_number = x_data.shape                # 获取样本数量以及属性特征的数量
    learning_rate, max_iters, weights = 0.001, 500, np.ones([attr_number, 1])    # 一些超参数和参数的初始化
    loss_list = list()
    for each_iter in range(max_iters):         # 更新迭代max_iters次
        sigmoid_result = sigmoid(np.matmul(x_data, weights))        # sigmoid处理 x*w
        difference = y_label - sigmoid_result                             # 计算损失值
        weights = weights + learning_rate * np.matmul(x_data.T, difference)    # 更新权重w向量
        loss = np.matmul(y_label.T, np.log(sigmoid_result)) + np.matmul((1 - y_label).T, np.log(1 - sigmoid_result))
        loss_list.append(loss.tolist()[0][0])
    return weights.getA(), loss_list

获取到了模型的最终训练的参数之后，就可以对结果进行可视化了，以便直观感受下Logistic回归的分类结果。为此，定义一个show_result方法来实现结果的可视化：

"""
    Author: Taoye
    微信公众号: 玩世不恭的Coder
    Explain: 可视化分类结果，也就是Logistic回归的可视化
    Parameters:
        x_data：样本集的属性特征
        y_label：样本集的标签
        weights：模型所需要的参数，也就是权重
"""
def show_result(x_data, y_label, weights):
    from matplotlib import pyplot as plt
    w_1, w_2, w_3 = weights[0][0], weights[1][0], weights[2][0]    # 获取权重参数
    min_x_1, min_x_2 = np.min(x_data, axis = 0)[:-1]               # 获取属性特征的最小值
    max_x_1, max_x_2 = np.max(x_data, axis = 0)[:-1]               # 获取属性特征的最大值
    line_x_1 = np.linspace(min_x_1 - 0.2, max_x_1 + 0.2, 1000)     # 决策直线的横坐标
    line_x_2 = (-w_3 - w_1 * line_x_1) / w_2                       # 决策直线的纵坐标
    plt.scatter(x_data[:, 0], x_data[:, 1], c = y_label)           # 绘制数据的散点图
    plt.plot(line_x_1, line_x_2)                                   # 绘制分类的决策直线

可视化分类结果如下所示：

上图主要包括俩个部分，一个是通过Logistic回归分类的结果，另一个是每次迭代之后损失值的变换情况。

从可视化的结果来看，这个分类效果相当不错，基本上所有的数据点都能够被正确分类，读者可根据需要来决定是否使用随机种子，从而观察不同数据集的分类效果。

从第二张图中，我们可以看出通过Logistic回归进行训练的过程中，损失函数的值是不断增大。而且我们还可以发现，其增大的斜率逐渐减小，尤其是在前面几次迭代过程中，这种现象尤为明显。当损失函数值提高到一定程度之后，这个时候基本属于饱和了，也就是说分类的过程基本结束了。

上图是损失函数的值的可视化结果，我们也可以从其每次迭代之后的具体值来观察这一变化趋势：

总共500次迭代，共输出500次损失值，起初的损失值为300多，每次迭代之后的损失值逐渐增大，且增大的速度在不断减小，最终的损失值停留在3左右达到饱和。这就是梯度上升算法所体现出来的效果，也就是说我们的损失函数值越小，我们梯度上升法优化的效果也就越明显。

完整代码：

import numpy as np

"""
    Author: Taoye
    微信公众号: 玩世不恭的Coder
    Explain: 通过NumPy准备数据集
    Return:
        x_data：样本的属性特征
        y_label：样本属性特征所对应的标签
"""
def establish_data():
#     np.random.seed(1)
    x_data = np.concatenate((np.add(np.random.randn(50, 2), [1.5, 1.5]),       
                                 np.subtract(np.random.randn(50, 2), [1.5, 1.5])),
                                 axis = 0)      # random随机生成数据，+ -1.5达到不同类别数据分隔的目的
    x_data = np.concatenate((x_data, np.ones([100, 1])), axis = 1)
    y_label = np.concatenate((np.zeros([50]), np.ones([50])), axis = 0)    # concatenate合并数据集
    return x_data, y_label

"""
    Author: Taoye
    微信公众号: 玩世不恭的Coder
    Explain: Sigmoid函数
    Parameters:
        in_data: sigmoid处理的输入数据
    Return:
        sigmoid_result：sigmoid函数处理之后的结果
"""
def sigmoid(in_data):
    return 1 / (1 + np.exp(-in_data))

"""
    Author: Taoye
    微信公众号: 玩世不恭的Coder
    Explain: Logistic回归的核心方法，主要使用到了梯度上升算法
    Parameters:
        x_data：样本集的属性特征
        y_label：样本集的标签
    Return:
        weights: 模型训练完成之后的w参数
"""
def gradient_ascent(x_data, y_label):
    x_data, y_label = np.mat(x_data), np.mat(y_label).T    # 转换成矩阵形式，方便操作
    data_number, attr_number = x_data.shape                # 获取样本数量以及属性特征的数量
    learning_rate, max_iters, weights = 0.001, 500, np.ones([attr_number, 1])    # 一些超参数和参数的初始化
    loss_list = list()
    for each_iter in range(max_iters):         # 更新迭代max_iters次
        sigmoid_result = sigmoid(np.matmul(x_data, weights))        # sigmoid处理 x*w
        difference = y_label - sigmoid_result                             # 计算损失值
        weights = weights + learning_rate * np.matmul(x_data.T, difference)    # 更新权重w向量
        loss = np.matmul(y_label.T, np.log(sigmoid_result)) + np.matmul((1 - y_label).T, np.log(1 - sigmoid_result))
        loss_list.append(loss.tolist()[0][0])
    return weights.getA(), loss_list

"""
    Author: Taoye
    微信公众号: 玩世不恭的Coder
    Explain: 可视化分类结果，也就是Logistic回归的可视化
    Parameters:
        x_data：样本集的属性特征
        y_label：样本集的标签
        weights：模型所需要的参数，也就是权重
"""
def show_result(x_data, y_label, weights):
    from matplotlib import pyplot as plt
    w_1, w_2, w_3 = weights[0][0], weights[1][0], weights[2][0]    # 获取权重参数
    min_x_1, min_x_2 = np.min(x_data, axis = 0)[:-1]               # 获取属性特征的最小值
    max_x_1, max_x_2 = np.max(x_data, axis = 0)[:-1]               # 获取属性特征的最大值
    line_x_1 = np.linspace(min_x_1 - 0.2, max_x_1 + 0.2, 1000)     # 决策直线的横坐标
    line_x_2 = (-w_3 - w_1 * line_x_1) / w_2                       # 决策直线的纵坐标
    plt.scatter(x_data[:, 0], x_data[:, 1], c = y_label)           # 绘制数据的散点图
    plt.plot(line_x_1, line_x_2)                                   # 绘制分类的决策直线

if __name__ == "__main__":
    x_data, y_label = establish_data()
    weights, loss_list = gradient_ascent(x_data, y_label)
    show_result(x_data, y_label, weights)
#     from matplotlib import pyplot as plt
#     plt.plot(np.arange(len(loss_result)), loss_list)

以上就是本篇文章中Logistic回归的全部内容了，我们在来总结一下Logistic回归实现的过程：

首先是通过分析推导，得到Logistic回归的损失函数\(L(w)\)（使用到了极大似然法）:

\[L(w)=\sum_{i=1}^ny_ilog(h_w(x))+(1-y_i)log(1-h_w(x))
\]

其次，为了将\(w\)和\(x\)内积之后的结果映射到0-1范围之内，以体现出概率的特性，我们引入了Sigmoid函数对内积结果进行处理：

\[h_w(x)=g(w^Tx)=\frac{1}{1+e^{-w^Tx}}
\]

引入了Sigmoid函数之后，将损失函数化简得到：

最后，因为我们需要不断地对\(w\)参数通过梯度上升算法进行更新迭代，所以我们需要对\(w\)进行求导，求导结果如下：

\[\frac{\partial L(w)}{\partial w}=\sum_{i=1}^n(y_i-h_w(x_i))x_i
\]

如此一来，我们就能通过该梯度值对\(w\)参数一步步的优化，优化方法主要是通过梯度上升算法：

\[w^{new}=w^{old}+\alpha x^T(y-h_w(x))
\]

即对\(w\)向量内部的每个元素进行更新：

当迭代次数达到一定程度时，最终更新迭代得到的\(w\)向量就是我们得到的参数结果，根据该参数就能构造一个模型将数据集分割开来，从而实现了数据的分类。

Logistic回归主要运用的是梯度上升算法，在上面案例实战的过程中，我们的随机生成的样本数据集不是特别的多，所以训练的速度还挺快。但是假定我们的训练样本数量比较多，这个时候的训练效率就比较的低下了。这个时候可能就需要对梯度上升算法进行一定的优化了，而该部分优化常用的方式是使用随机梯度上升算法，限于篇幅和时间原因，我们后面有机会再来肝。

这是手撕机器学习系列文章的第九篇了，也差不多接近尾声了，初步计划这周完成吧，因为后面还有很多任务到现在还没有开始，想到这，Taoye的眼角有不知觉地。。。

当然了，吾生也有涯，而知也无涯，以有涯随无涯，殆已。学习本身就是一个无止境的过程，目前所讲到的机器学习算法以及所涉及到的知识也只是该领域的冰山一角，我们最重要的是要保持一颗不断积极进取的心。

我是Taoye，爱专研，爱分享，热衷于各种技术，学习之余喜欢下象棋、听音乐、聊动漫，希望借此一亩三分地记录自己的成长过程以及生活点滴，也希望能结实更多志同道合的圈内朋友，更多内容欢迎来访微信公主号：玩世不恭的Coder。

我们下期再见，拜拜~~~

参考资料：

[1] 《机器学习实战》：Peter Harrington 人民邮电出版社
[2] 《统计学习方法》：李航第二版清华大学出版社

Tags: numpy Python 数据挖掘机器学习算法