机器学习算法笔记（二）：逻辑回归

在学习机器学习的过程中，结合数学推导和手写实现，可以加深对相关算法的认识。本部分教程将基于python实现机器学习的常用算法，来加强对算法的理解以及coding能力，仅供学习交流使用，请勿随意转载。

本篇继续逻辑回归算法的学习，全文分为三个部分：

• 数学推导
• python实现
• 逻辑回归优缺点分析

一、逻辑回归的数学推导

​ 逻辑回归（LogisticRegression）名为回归，实为分类。逻辑回归可也可称为对数几率回归，是一种广义线性模型。在上篇文章中我们学习到基础的线性回归模型，其针对的是标签为连续值的数据集，而逻辑回归则是通过sigmoid函数，将线性回归模型中连续的标签转化成分类标签的模型。

​ 相较于线性回归的因变量 y 为连续值，逻辑回归的因变量则是一个 0/1 的二分类值，这就需要我们建立一种映射将原先的实值转化为 0/1 值。

公式1-1：Logistic回归模型的估计概率

\[\hat{p}=h_\theta(x)=\sigma(\theta^{\mathrm{T}}x) \]

​ 在此等式中：

• \(\hat{p}\)是回归模型通过sigmoid函数映射得到的估计概率，取值为0-1
• \(h_\theta\)是假设函数，使用模型参数\(\theta\)\(x_i\)是第\({i}\)个特征值
• \(\sigma\)是sigmoid函数
• \(\theta^{\mathrm{T}}\cdot x\)是向量\(\theta\)和\(x\)的点积，其值等于\(\theta_0x_0+\theta_1x_1+\theta_2x_2+...+\theta_nx_n\)

sigmoid函数如公式1-2所示

公式1-2：sigmoid函数

\[\sigma=\frac{1}{1+e^{-x}} \]

图1-1 sigmoid函数

​ 图1-1给出了Sigmoid函数在不同坐标尺度下的两条曲线图。当x为0时，Sigmoid函数值为0.5。随着x的增大，对应的Sigmoid值将逼近于1；而随着x的减少，Sigmoid值将逼近于0.如果横坐标刻度足够大，Sigmoid函数看起来很像一个阶跃函数。

　　解释Sigmoid函数：将任意的输入映射到了 [0, 1]区间，我们在线性回归中可以得到一个预测值，再将该值映射到 Sigmoid函数中这样就完成了由值到概率的转换，也就是分类任务。

公式1-3：逻辑回归模型预测

\[\hat{y}=\left\{ \begin{aligned} 0, if\quad\hat{p}<0.5 \\ 1, if\quad\hat{p}\geq0.5 \end{aligned} \right. \]

​ 在此等式中：

• \(\hat{y}\)是通过估计概率得到的预测值，为0或者1

• \(\hat{p}\)是回归模型通过sigmoid函数映射得到的估计概率，取值为0-1

（一）损失函数——交叉熵损失

​ 关于交叉熵的推导，可以根据极大似然法/后验概率获得

• 极大似然法

​ 假设有N个样本，样本的标签只有0和1两类，可以用极大似然估计法估计模型参数，从而得到逻辑回归模型。

​ 由于sigmoid函数的特性，我们可以做出如下假设：设\(y_i=1\)的概率为\(p_i\)，\(y_i=0\)的概率则为\(1-p_i\)，那么观测的概率为：

公式1-4：样本i观测概率

\[p(y_i)=p_i^{y_i}\times(1-p_i)^{(1-y_i)} =(h_\theta(x_i)^{y_i}(1-h_\theta(x_i))^{(1-y_i)}) \]

​ 假定样本间相互独立，则似然函数为式1-5：

公式1-5：预测值似然函数

\[\begin{align} &L(\theta)=\prod_{i=1}^{n}p(y_i|x_i;\theta)\\ &=\prod_{i=1}^{n}(h_\theta(x_i)^{y_i}(1-h_\theta(x_i))^{(1-y_i)}) \end{align} \]

​ 为了求解最大似然函数，我们对式1-5取对数，得到交叉熵损失：

公式1-6：交叉熵损失

\[\begin{align} l(\theta)&=\log L(\theta)\\ &=\sum_{i=1}^{n}y_i\log h(x_i)+(1-y_i)\log (1-h(x_i))\\ J(\theta)&=-\frac{1}{m}l(\theta) \end{align}\\ \]

• 后验概率

  ​ 由公式1-1以及1-2可知逻辑回归模型的基本形式为：

公式1-7：逻辑回归基本形式

\[y=\frac{1}{1+e^{(-\theta^{\mathrm{T}}x)}} \]

​ 对上式做一下变换：

\[\ln \frac{y}{1-y}=\theta^{\mathrm{T}}x \]

​ 将y视为类后验概率\(p(y=1|x)\)，则上式可以写为：

\[\ln \frac{p(y=1|x)}{p(y=0|x)}=\theta^{\mathrm{T}}x \]

​ 则有：

\[p(y=1|x)=\frac{e^{\theta^{\mathrm{T}}x}}{1+e^{\theta^{\mathrm{T}}x}}=\hat{y}\\ p(y=0|x)=\frac{1}{1+e^{\theta^{\mathrm{T}}x}}=1-\hat{y} \]

​ 将上式进行简单综合，可写成如下形式：

\[p(y|x)=\hat{y}^y(1-\hat{y})^{1-y} \]

​ 这与极大似然法中预测概率一致，将其写成对数形式，就可以得到交叉熵损失。

（二）求解损失函数

​ 为了得到使损失函数最小的\(\theta\)值，可以使用梯度下降法进行迭代调整参数从而使成本函数最小化。

​ 首先对损失函数进行求导

\[\begin{align} l(\theta)&=\sum_{i=1}^{n}y_i\log h(x_i)+(1-y_i)\log (1-h(x_i))\\ \frac{\partial l(\theta)}{\partial\theta_j}&=\frac{\partial l(\theta)}{\partial\sigma(\theta^{\mathrm{T}}x)}\cdot\frac{\partial\sigma(\theta^{\mathrm{T}}x)}{\partial\theta^{\mathrm{T}}x}\cdot\frac{\partial\theta^{\mathrm{T}}x}{\partial\theta_j}\\ 其中：\\ \frac{\partial l(\theta)}{\partial\sigma(\theta^{\mathrm{T}}x)}&=y\cdot\frac{1}{\sigma(\theta^{\mathrm{T}}x)}+(1-y)\cdot\frac{1}{1-\sigma(\theta^{\mathrm{T}}x)}\cdot(-1)\\ 而由sigmoid函数：\\ \sigma^{'}(x)&=\frac{\mathrm{d}}{\mathrm{d}x}\frac{1}{1+e^{-x}}\\ &=\frac{1}{(1+e^{-x})^2}\cdot(e^{-x})\\ &=\frac{1}{1+e^{-x}}\cdot(1-\frac{1}{1+e^{-x}})\\ &=\sigma(x)\cdot(1-\sigma(x))\\ 则有：\\ \frac{\partial\sigma(\theta^{\mathrm{T}}x)}{\partial\theta^{\mathrm{T}}x}&=\sigma(\theta^{\mathrm{T}}x)\cdot(1-\sigma(\theta^{\mathrm{T}}x))\\ 最后：\\ \frac{\partial\theta^{\mathrm{T}}x}{\partial\theta_j}&=\frac{\partial(\theta_1x_1+\theta_2x_2+\cdots+\theta_nx_n)}{\partial\theta_j}=x_j \end{align} \]

公式1-8：逻辑回归损失函数求导

\[\frac{\partial l(\theta)}{\partial \theta_j}=(y-h_\theta(x))x_j\\ \frac{\partial J(\theta)}{\partial \theta_j}=\frac{1}{n}\sum_{i=1}^{n}(\sigma(\theta^{\mathrm{T}}x^{(i)})-y^{(i)})x_j^{(i)} \]

注：当\(j=0\)时，\(x_0=1\)

梯度下降法的伪代码：

for 每一步 in 所有的训练批次:
使用整批训练数据获得梯度

​ 进行参数更新

二、Python实现

​ 在对线性回归的数学原理进行大致了解之后，我们开始进行算法的编写。我们首先构造数据集，然后根据梯度下降法的思路，我们需要在每个批次中对参数进行更新，这就需要我们在训练之初对参数进行初始化，然后根据损失函数获得损失函数的参数求导结果——即梯度，基于梯度对参数进行更新。最后我们使用交叉验证获得更加稳健的参数估计值。

​ 导入需要的相关包

import numpy as np
from sklearn.utils import shuffle
from sklearn.datasets import load_diabetes

（一）数据集构建

​ 本文基于sklearn构建分类数据集，进行简单说明。

    def createData():
        """
        :return:
        X_train-(90, 2), y_train-(90, 1), X_test-(10, 2) y_test-(10, 1)
        """
        X, labels = make_classification(n_samples=100, n_features=2, n_redundant=0, n_informative=2, random_state=42, n_clusters_per_class=2)
        labels = labels.reshape((-1, 1))
        offset = int(X.shape[0]*0.9)
        X_train, y_train = X[:offset], labels[:offset]
        X_test, y_test = X[offset:], labels[offset:]
        return X_train, y_train, X_test, y_test

​ 将分布结果可视化：

（二）初始化参数

​ 根据输入数据的维度，返回维度分别为（dims，1）的参数w和（1，）的截距项。

def initializeParams(dims):
    """
    :param dims: X的维度2
    :return:
     w-(2, 1)
     b-(1,)
     """
    w = np.zeros((dims, 1))
    b = 0
    return w, b

（三）sigmoid函数

​ 定义sigmoid函数，实现由回归值到分类概率的映射。

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

（四）损失函数

​ 通过sigmoid函数得到\(\hat{y}\)，并计算损失函数，求出对应参数偏导数。

def logistic(X, y, w, b):
    num_train = X.shape[0]
    num_feature = X.shape[1]

    a = self.sigmoid(np.dot(X, w) + b)
    cost = - 1 / num_train * np.sum(y * np.log(a) + (1 - y) * np.log(1 - a))
    dw = np.dot(X.T, (a - y)) / num_train
    db = np.sum(a - y) / num_train
    cost = np.squeeze(cost)
    return a, cost, dw, db

（五）模型训练

​ 基于梯度下降对模型进行训练，不断更新参数。

def logisticTrain(X, y, learning_rate, epochs):
    w, b = self.initializeParams(X.shape[1])
    costs = []
    for epoch in range(1, epochs+1):
        a, cost, dw, db = self.logistic(X, y, w, b)
        w += -learning_rate * dw
        b += -learning_rate * db
        if epoch % 100 ==0:
            costs.append(cost)
            print(f'epoch: {epoch} cost: {cost}')

            params = {
                'w': w,
                'b': b
            }

            grads = {
                'dw': dw,
                'db': db
            }
	return costs, params, grads

（六）完整代码

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# author: GeekThomas
# time: 2021/3/22
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets.samples_generator import make_classification

class LogisticRegression:
    def __init__(self):
        pass

    def createData(self):
        """
        :return:
        X_train-(90, 2), y_train-(90, 1), X_test-(10, 2) y_test-(10, 1)
        """
        X, labels = make_classification(n_samples=100, n_features=2, n_redundant=0, n_informative=2, random_state=42, n_clusters_per_class=2)
        labels = labels.reshape((-1, 1))
        offset = int(X.shape[0]*0.9)
        X_train, y_train = X[:offset], labels[:offset]
        X_test, y_test = X[offset:], labels[offset:]
        return X_train, y_train, X_test, y_test

    def initializeParams(self, dims):
        """
        :param dims: X的维度-2
        :return:
        w-(2, 1), b-(1,)
        """
        w = np.zeros((dims, 1))
        b = 0
        return w, b

    def sigmoid(self, x):
        return 1 / (1 + np.exp(-x))

    def logistic(self, X, y, w, b):
        num_train = X.shape[0]
        num_feature = X.shape[1]

        a = self.sigmoid(np.dot(X, w) + b)
        cost = - 1 / num_train * np.sum(y * np.log(a) + (1 - y) * np.log(1 - a))
        dw = np.dot(X.T, (a - y)) / num_train
        db = np.sum(a - y) / num_train
        cost = np.squeeze(cost)
        return a, cost, dw, db

    def logisticTrain(self, X, y, learning_rate, epochs):
        w, b = self.initializeParams(X.shape[1])
        costs = []
        for epoch in range(1, epochs+1):
            a, cost, dw, db = self.logistic(X, y, w, b)
            w += -learning_rate * dw
            b += -learning_rate * db
            if epoch % 100 ==0:
                costs.append(cost)
                print(f'epoch: {epoch} cost: {cost}')

            params = {
                'w': w,
                'b': b
            }

            grads = {
                'dw': dw,
                'db': db
            }

        return costs, params, grads

    def predict(self, X, params):
        w, b = params['w'], params['b']
        y_pred = self.sigmoid(np.dot(X, w) + b)
        for i in range(len(y_pred)):
            if y_pred[i] < 0.5:
                y_pred[i] = 0
            else:
                y_pred[i] = 1
        return y_pred

    def accuracy(self, y_test, y_pred):
        return np.sum(y_pred == y_test) / len(y_test)

    def plotLogistic(self, X_train, y_train, params):
        fig, axes = plt.subplots(figsize=(5, 5))
        axes.scatter(X_train[y_train.squeeze()==0, 0], X_train[y_train.squeeze()==0, 1], s=100, c='red')
        axes.scatter(X_train[y_train.squeeze()==1, 0], X_train[y_train.squeeze()==1, 1], s=100, c='green')
        x = np.arange(-1.5, 3, 0.1)
        y = (-params['b'] - params['w'][0] * x) / params['w'][1]
        axes.plot(x, y)
        plt.xlabel('x1')
        plt.xlabel('x2')
        plt.show()


if __name__ == '__main__':
    model = LogisticRegression()
    X_train, y_train, X_test, y_test = model.createData()
    print(X_train.shape, y_train.shape, X_test.shape, y_test.shape)

    costs, params, grads = model.logisticTrain(X_train, y_train, 0.01, 1000)
    print(params)

    y_train_pred = model.predict(X_train, params)
    accuracy_score_train = model.accuracy(y_train, y_train_pred)
    print(f'train accuracy is {round(accuracy_score_train, 4)}')

    y_test_pred = model.predict(X_test, params)
    accuracy_score_test = model.accuracy(y_test, y_test_pred)
    print(f'test accuracy is {round(accuracy_score_test, 4)}')

    model.plotLogistic(X_train, y_train, params)

三、逻辑回归模型优缺点分析

（一）优点

1. 实现简单，广泛的应用于工业问题上；
2. 分类时计算量非常小，速度很快，存储资源低；
3. 便利的观测样本概率分数；
4. 对逻辑回归而言，多重共线性并不是问题，它可以结合L2正则化来解决该问题；
5. 计算代价不高，易于理解和实现；

（二）缺点

1. 当特征空间很大时，逻辑回归的性能不是很好；
2. 容易欠拟合，一般准确度不太高
3. 不能很好地处理大量多类特征或变量；
4. 只能处理二分类问题（在此基础上衍生出来的softmax可以用于多分类），且必须线性可分；
5. 对于非线性特征，需要进行转换；

---

​ 以上就是本文的全部内容——基于Python实现逻辑回归。

参考资料：

1.机器学习实验室：机器学习公式推导与代码实现