【机器学习】逻辑回归实现与应用

介绍

逻辑回归（Logistic Regression），又叫逻辑斯蒂回归，是机器学习中一种十分基础的分类方法，由于算法简单而高效，在实际场景中得到了广泛的应用。本次实验中，我们将探索逻辑回归的原理及算法实现，并使用 scikit-learn 构建逻辑回归分类预测模型。

Sigmoid分布函数

因为今天要解决的是一个二分类问题，所以先介绍此函数

$$f(z)=\frac{1}{1+e^{-z}}$$
python实现：

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

这个图像呈现出完美的 S 型（Sigmoid 的含义）。它的取值仅介于 $0$ 和 $1$ 之间，且关于 $z=0$ 轴中心对称。同时当 $z$ 越大时，$y$ 越接近于 $1$，而 $z$ 越小时，$y$ 越接近于 $0$。如果我们以 $0.5$ 为分界点，将 $>0.5$ 或 $<0.5$ 的值分为两类，这不就是解决 $0-1$ 二分类问题的完美选择嘛。

逻辑回归模型

如果一组连续随机变量符合 Sigmoid 函数样本分布，就称作为逻辑分布。逻辑分布是概率论中的定理，是一种连续型的概率分布。
在逻辑回归中，定义：
$$z_i=w_0{x}_0+w_1{x}_1+\cdots +w_i{x}_i=w^{T}x$$

$$f(z_i)=\frac{1}{1+e^{-z_i}}$$
即：
$$h_w(x)=f(w^{T}x)=\frac{1}{1+e^{-w^{T}x}}$$

由于目标值 $y$ 只有 0 和 1 两个值，那么如果记 $y=1$ 的概率为 $h_w(x)$，则此时 $y=0$ 的概率为 $1-h_w(x)$。那么，我们可以记作逻辑回归模型条件概率分布：
$$P(Y=y|x)=\{\begin{array}{rl}{h}_w(x)& ,y=1\\ 1-h_w(x)& ,y=0\end{array}$$

其可等价写为似然函数：
$$P(y|x;w)={\left(h_w(x)\right)}^y{\left(1-h_w(x)\right)}^{1-y}$$

对于 $i$ 个样本的总概率而言实际上可以看作单样本概率的乘积，记为 $L(w)$：
$$L(w)=\prod _{i=1}^m{\left(h_w\left(x^{(i)}\right)\right)}^{y^{(i)}}{\left(1-h_w\left(x^{(i)}\right)\right)}^{1-y^{(i)}}$$

由于连乘表示起来非常复杂，我们应用数学技巧，即两边取对数将连乘转换为连加的形式，即：
$$\log L(w)=\sum _{i=1}^m\left[y^{(i)}\log h_w\left(x^{(i)}\right)+\left(1-y^{(i)}\right)\log \left(1-h_w\left(x^{(i)}\right)\right)\right]\text{\hspace{1em}(1)}$$

对数损失函数

函数(1)衡量了事件发生的总概率。根据最大似然估计原理，只需要通过对 𝐿(𝑤) 求最大值，即得到 𝑤 的估计值。而在机器学习问题中，我们需要一个损失函数，并通过求其最小值来进行参数优化。所以，对数似然函数取负数就可以被作为逻辑回归的对数损失函数：
$$J(w)=-\frac{1}{m}\sum _{i=1}^m\left[y^{(i)}\log h_w\left(x^{(i)}\right)+\left(1-y^{(i)}\right)\log \left(1-h_w\left(x^{(i)}\right)\right)\right]$$

python实现:

def loss(h, y):
    loss = (-y * np.log(h) - (1 - y) * np.log(1 - h)).mean()
    return loss

梯度下降法

之前我已经写过一期梯度下降法相关的blog，名为【机器学习】梯度下降法求解线性回归参数，但是其中对于数学原理的介绍并不是很详尽，这次我将结合上述本次公式对梯度下降法再次讲解，其实梯度的概念是高等数学中的内容，很明显重点就是求导，因为我们想知道的是如何下降最快，这样只需要沿着上升最快的反方向走就可以了，也就是沿着梯度的反向走，按照步长一点一点走，不然就有可能跳过最小值。

梯度是一个向量，它表示某一函数在该点处的方向导数沿着该方向取得最大值，即函数在该点处沿着该方向（此梯度的方向）变化最快，变化率最大（为该梯度的模）。简而言之，对于一元函数而言，梯度就是指在某一点的导数。而对于多元函数而言，梯度就是指在某一点的偏导数组成的向量。

我们先针对公式（1）化简：
$$\begin{array}{rl}J(w)& =-\frac{1}{m}\sum _{i=1}^m\left[y^{(i)}(-\log (1+e^{-w^T{x}^{(i)}}))+(1-y^{(i)})(-w^T{x}^{(i)}-\log (1+e^{-w^T{x}^{(i)}}))\right]\\ & =-\frac{1}{m}\sum _{i=1}^m\left[-y^{(i)}\log (1+e^{-w^T{x}^{(i)}})-w^T{x}^{(i)}-\log (1+e^{-w^T{x}^{(i)}})+w^T{x}^{(i)}{y}^{(i)}+y^{(i)}\log (1+e^{-w^T{x}^{(i)}})\right]\\ & =-\frac{1}{m}\sum _{i=1}^m\left[-w^T{x}^{(i)}-\log (1+e^{-w^T{x}^{(i)}})+w^T{x}^{(i)}{y}^{(i)}\right]\\ & =-\frac{1}{m}\sum _{i=1}^m\left[-\log (e^{w^T{x}^{(i)}})-\log (1+e^{-w^T{x}^{(i)}})+w^T{x}^{(i)}{y}^{(i)}\right]\\ & =-\frac{1}{m}\sum _{i=1}^m\left[-\log (e^{w^T{x}^{(i)}}(1+e^{-w^T{x}^{(i)}}))+w^T{x}^{(i)}{y}^{(i)}\right]\\ & =-\frac{1}{m}\sum _{i=1}^m\left[-\log (e^{w^T{x}^{(i)}}+1)+w^T{x}^{(i)}{y}^{(i)}\right]\end{array}\text{\hspace{1em}(2)}$$

接下来，我们针对公式（2）求导：
$$\begin{array}{rl}\frac{\partial J}{\partial w}& =-\frac{1}{m}\sum _{i=1}^m\left[-\frac{x^{(i)}{e}^{w^T{x}^{(i)}}}{e^{w^T{x}^{(i)}}+1}+x^{(i)}{y}^{(i)}\right]\\ & =-\frac{1}{m}\sum _{i=1}^m\left[-\frac{e^{w^T{x}^{(i)}}}{e^{w^T{x}^{(i)}}+1}+y^{(i)}\right]x^{(i)}\\ & =-\frac{1}{m}\sum _{i=1}^m\left[-\frac{e^{-w^T{x}^{(i)}}{e}^{w^T{x}^{(i)}}}{e^{-w^T{x}^{(i)}}(e^{w^T{x}^{(i)}}+1)}+y^{(i)}\right]x^{(i)}\\ & =-\frac{1}{m}\sum _{i=1}^m\left[-\frac{1}{1+e^{-w^T{x}^{(i)}}}+y^{(i)}\right]x^{(i)}\end{array}$$

即：
$$\frac{\partial J}{\partial w}=-\frac{1}{m}\sum _{i=1}^m\left[-h_w(x^{(i)})+y^{(i)}\right]x^{(i)}$$

向量的形式表达为：
$$\frac{\partial J}{\partial w}=\frac{1}{m}{x}^T(h_w(x)-y)$$

当我们得到梯度的方向，然后乘以一个常数 $\alpha$，就可以得到每次梯度下降的步长(上图箭头的长度)。最后，通过多次迭代，找到梯度变化很小的点，也就对应着损失函数的极小值了。其中，常数 $\alpha$ 往往也被称之为学习率 Learning Rate。执行权重更新的过程为：
$$w\leftarrow w-\alpha \frac{\partial J}{\partial w}$$
python实现：

def gradient(X, h, y):
    gradient = np.dot(X.T, (h - y)) / y.shape[0]
    return gradient

逻辑回归实现

加载数据

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
df = pd.read_csv(
    "https://labfile.oss.aliyuncs.com/courses/1081/course-8-data.csv", header=0)  # 加载数据集
df.head()  # 预览前 5 行数据

函数代码汇总

def sigmoid(z):
    # Sigmoid 分布函数
    sigmoid = 1 / (1 + np.exp(-z))
    return sigmoid


def loss(h, y):
    # 损失函数
    loss = (-y * np.log(h) - (1 - y) * np.log(1 - h)).mean()
    return loss


def gradient(X, h, y):
    # 梯度计算
    gradient = np.dot(X.T, (h - y)) / y.shape[0]
    return gradient


def Logistic_Regression(x, y, lr, num_iter):
    # 逻辑回归过程
    intercept = np.ones((x.shape[0], 1))  # 初始化截距为 1
    x = np.concatenate((intercept, x), axis=1)
    w = np.zeros(x.shape[1])  # 初始化参数为 0

    for i in range(num_iter):  # 梯度下降迭代
        z = np.dot(x, w)  # 线性函数
        h = sigmoid(z)  # sigmoid 函数

        g = gradient(x, h, y)  # 计算梯度
        w -= lr * g  # 通过学习率 lr 计算步长并执行梯度下降

        l = loss(h, y)  # 计算损失函数值

    return l, w  # 返回迭代后的梯度和参数

逻辑回归

x = df[['X0', 'X1']].values
y = df['Y'].values
lr = 0.01  # 学习率
num_iter = 30000  # 迭代次数

# 训练
L = Logistic_Regression(x, y, lr, num_iter)
L

逻辑回归 scikit-learn 实现

LogisticRegression(penalty='l2', dual=False, tol=0.0001, C=1.0,
fit_intercept=True, intercept_scaling=1, class_weight=None,
random_state=None, solver='liblinear', max_iter=100,
multi_class='ovr', verbose=0, warm_start=False, n_jobs=1) 

介绍其中几个常用的参数，其余使用默认即可：

• penalty: 惩罚项，默认为 $L_2$ 范数。
• dual: 对偶化，默认为 False。
• tol: 数据解算精度。
• fit_intercept: 默认为 True，计算截距项。
• random_state: 随机数发生器。
• max_iter: 最大迭代次数，默认为 100。

另外，solver 参数用于指定求解损失函数的方法。默认为 liblinear（0.22 开始默认为
lbfgs），适合于小数据集。除此之外，还有适合多分类问题的 newton-cg, sag, saga 和 lbfgs
求解器。这些方法都来自于一些学术论文，有兴趣可以自行搜索了解。

from sklearn.linear_model import LogisticRegression

model = LogisticRegression(tol=0.001, max_iter=10000, solver='liblinear')  # 设置数据解算精度和迭代次数
model.fit(x, y)
model.coef_, model.intercept_