一天一个机器学习小知识——逻辑回归

---

前言

前面几小节介绍的线性回归、Lasso、Ridge以及弹性网都是回归模型，但是现实生活中还是会有很多分类问题，因此本文就介绍一个机器学习中最常见的分类模型——逻辑回归。逻辑回归是最经典的分类模型之一，一方面，它保持了线性回归的可解释性，从参数的大小可以知道每个特征对结果的影响程度；另一方面，它的输出具有概率意义，可以为很多决策提供参考。

一、算法推导

李航老师的《统计学习方法》中提到，统计学习方法都是由模型、策略和算法构成的，因此本文在算法推导也主要从这三部分进行展开讨论。

1.模型

虽然逻辑回归带有“回归”二字，但其实是一个分类模型，不过既然它这么起名字，说明肯定跟线性回归有很大的关系。其实逻辑回归只不过是将线性回归的输出结果进行压缩转换，本质上还是一个线性模型。我们知道普通的线性回归如下所示：
$$z=f(x)=w_0^{T}x+b=w^{T}x$$
它的输出范围是$(-\infty ,+\infty )$，不符合分类问题的输出格式，因此我们需要对这个结果进行压缩，也就是做一个映射，而sigmoid函数是最常用的映射函数。
$$y=sigmoid(z)=\frac{1}{1+e^{-z}}=\frac{1}{1+e^{-{\mathbf{w}}^T\mathbf{x}}}$$

通过sigmoid函数，将线性回归的输出从$(-\infty ,+\infty )$映射至$(0,1)$，这个结果可以看成是输出为类别1的概率。如果y>0.5（这里0.5可以改成其他值），说明结果为类别1的概率超过50%，那么将结果判定为类别1，否则判定为类别0。

2.策略

由于经过sigmoid函数的转换，按照线性回归那样直接以均方误差作为损失函数对于逻辑回归不再适用。既然上面提到逻辑回归的输出具有概率意义，所以我们可以尝试从概率的角度推导出逻辑回归的损失函数。
$$y=sigmoid(z)=\frac{1}{1+e^{-z}}=\frac{1}{1+e^{-{\mathbf{w}}^T\mathbf{x}}}$$
输出的是正类的概率，不妨令其为$p$，那么预测为负类的概率为$1-p$,

$$\begin{array}{rl}& P(y|\mathbf{x})=\{\begin{array}{rl}p,y& =1\\ 1-p,y& =0\end{array}\\ & \\ & 可以合并成：P\left(y_i|{\mathbf{x}}_i\right)=p^{y_i}(1-p{)}^{1-y_i}\end{array}$$
其中$P(y_i|x_i)$是事件（样本）$i$ 发生的概率，当$y_i=1$时，$P(y_i|x_i)=p_i$刚好等于正类发生的概率，$y_i=0$时，$P(y_i|x_i)=1-p_i$刚好等于负类发生的概率。整个事件一共有N的样本，那么整个事件的发生概率为（似然函数）：
$$\begin{array}{rl}L(w)=P_{总}& =P\left(y_1|{\mathbf{x}}_1\right)P\left(y_2|{\mathbf{x}}_2\right)P\left(y_3|{\mathbf{x}}_3\right)\dots P\left(y_N|{\mathbf{x}}_N\right)\\ & =\prod _{n=1}^N{p}^{y_n}(1-p{)}^{1-y_n}\end{array}$$我们应该选择一个参数$w$，使得事件出现的概率最大，即可以通过最大化上面的似然函数来求得$w$的最优解。但是似然函数涉及连乘，直接计算会出现“下溢”，并且不好求最小值。因此可以对似然函数做一个线性转换，转换后并不会改变其最优解，比如说取对数，得到对数似然：
$$\begin{array}{rl}ln(L(w))=\ln \left(P_{{\xi }_1}\right)& =\ln \left(\prod _{n=1}^N{p}_n^{y_n}(1-p_n{)}^{1-y_n}\right)\\ & =\sum _{n=1}^N\ln \left(p_n^{y_n}(1-p_n{)}^{1-y_n}\right)\\ & =\sum _{n=1}^N\left(y_n\ln (p_n)+\left(1-y_n\right)\ln (1-p_n)\right)\end{array}$$其中$y_n$是已知的，而$p_n$则与参数$w$有关，只需要使用梯度下降法求$ln(L(w))$最大的时候对应的$w$即为最优$w^{\ast }$。
但是到这里还没有结束，因为我们其实想求的是损失函数，然后最小化损失函数得到最优参数解，而目前推导出来的是似然函数。事实上，最大似然和最小损失是对应的，我们只需要在对数似然前面加上一个负号，就可以转化为求最小值问题。即损失函数可以写成
$$J(w)=-\sum _{n=1}^N\left(y_n\ln (p_n)+\left(1-y_n\right)\ln (1-p_n)\right)$$你可以把它看作是交叉熵（cross-entropy）(其实不完全等价于交叉熵，除非$y_i$对应的是真实分布)，但是为了理解直观理解，你可以先把他看成是交叉熵。交叉熵是分类中最常用的损失函数，没有之一。
因此最大化似然函数的问题就可以转变为我们所熟悉的最小化损失函数的问题
$$w^{\ast }=argmin(J(w))=argmin(-\sum _{n=1}^N\left(y_n\ln (p_n)+\left(1-y_n\right)\ln (1-p_n)\right))$$

3.算法

算法部分就是通过最小化$J(w)$来求出参数$w$。这里采用的梯度下降法。其实前面几个小节也提到过梯度下降法，但是一直没有对它的原理进行一个详细的说明，所以本小节将对这个优化算法进行深入的讨论。
梯度下降法是一个迭代算法， 通过不断的朝着负梯度更新参数，最终使得$J(w)$达到最小值。这里有个关键点：为什么一定要朝着负梯度更新参数？其实是因为负梯度方向是下降最快的方法。下面通过图形和简单公式推导来证明这个假设。
（1）图形理解

引用网友的一段话“首先来看看梯度下降的一个直观的解释。比如我们在一座大山上的某处位置，由于我们不知道怎么下山，于是决定走一步算一步，也就是在每走到一个位置的时候，求解当前位置的梯度，沿着梯度的负方向，也就是当前最陡峭的位置向下走一步，然后继续求解当前位置梯度，向这一步所在位置沿着最陡峭最易下山的位置走一步。这样一步步的走下去，一直走到觉得我们已经到了山脚。当然这样走下去，有可能我们不能走到山脚，而是到了某一个局部的山峰低处。从上面的解释可以看出，梯度下降不一定能够找到全局的最优解，有可能是一个局部最优解。当然，如果损失函数是凸函数，梯度下降法得到的解就一定是全局最优解。”
（2）公式简单推导
我这里使用公式来对负梯度下降最快这一说法进行简单推导，不一定很严谨。我们对$J(w)$进行泰勒展开$$J(w^{i+1})=J(w^i)+J^{{}^{\prime }}(w^i)(w^{i+1}-w^i)+o(w^{i+1}-w^i)$$ 在时刻 $t=i+1$,$w^i$ 是已知的，我们需要做的就是通过 $w^i$ 来迭代求出 $w^{i+1}$ 。对上式进行简单的移项，并省略掉余项 （$o(w^{i+1}-w^i)$） 。
$$J(w^i)-J(w^{i+1})=-J^{{}^{\prime }}(w^i)(w^{i+1}-w^i)$$ 要使损失函数 $J(w)$ 下降的最快，那么 $J(w^i)-J(w^{i+1})$ 达到最大，等价于 $-J^{{}^{\prime }}(w^i)(w^{i+1}-w^i)$ 最大，即 $J^{{}^{\prime }}(w^i)(w^{i+1}-w^i)$ 最小。要知道 $J^{{}^{\prime }}(w^i)$和$(w^{i+1}-w^i)$ 都是向量，都是有方向的，我们通过向量的乘法可以知道 $$J^{{}^{\prime }}(w^i)(w^{i+1}-w^i)=|J^{{}^{\prime }}(w^i)|\ast |(w^{i+1}-w^i)|\ast cos\theta$$所以在 $|J^{{}^{\prime }}(w^i)|和|(w^{i+1}-w^i)|$ 固定的情况下，要使 $J^{{}^{\prime }}(w^i)(w^{i+1}-w^i)$ 最小， $cos\theta$ 必须等于-1，即 $J^{{}^{\prime }}(w^i)$和$(w^{i+1}-w^i)$ 方向相反： $(w^{i+1}-w^i)=-J^{{}^{\prime }}(w^i)$ ，这就证明了沿着负梯度方向更新参数会使得 $J(w)$ 下降最快。
所以 $w^{i+1}=w^i-\lambda J^{{}^{\prime }}(w^i)$ 这个就是梯度下降的最终公式。
下面对梯度下降法进行一个简单的模拟：

import numpy as np
import matplotlib.pyplot as plt
def f(x):
    return x**2 + 1

x = np.random.randint(low=-100,high=100,size=1000)
x.sort()
y = f(x)

def Gradient_descent(r):
    w_i = 10  # 初始值
    while True:
        w_i_1 = w_i - r * (2* w_i)  # 更新公式：w[i+1] = w[i] - r * f'(w[i])
        if abs(w_i_1-w_i)<= 0.001 or abs(f(w_i_1) - f(w_i)) <= 0.001:  # 停止条件
            break 
        w_i = w_i_1
    return w_i

w_ = Gradient_descent(0.1)
plt.plot(x,y)
plt.scatter(w_, f(w_),color='r',label='argmin(loss)')
plt.legend(loc=9)

对于这种简单的凸优化问题，梯度下降很快就能找到最优解。

二、应用场景

逻辑回归的应用大致分为三类
（1）用于分类：适合做很多分类算法的基础组件。比如说和LGB一起组成Stacking模型。
（2）用于预测：预测事件发生的概率。比如说预测一个客户的信用卡违约的概率。
（3）用于分析：某一因素对事件发生的影响分析。比如说在（2）的场景中，不仅分析了客户违约的概率，还可以通过模型的系数来确定哪个因素（特征）对违约的可能性影响最大。

三、代码实现

这里主要使用sklearn这个库来实现机器学习算法。

1.导入相关库

from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.model_selection import cross_val_score
from sklearn.datasets import load_breast_cancer
from sklearn.metrics import accuracy_score
import pandas as pd
import numpy as np

2.读取样例数据

data = load_breast_cancer() # 获取样例数据，这里的数据乳腺癌检测,y=0代表没有乳腺癌，y=1代表有
X = pd.DataFrame(data.data,columns=data.feature_names)
y = data.target

3.划分训练集和测试集

X_train, X_test, y_train, y_test =train_test_split(X,y,test_size=0.3,random_state=0)

4.建立模型

model = LogisticRegression().fit(X_train,y_train)
y_pred = model.predict(X_test)
最终模型参数

5.评估模型

由于是分类问题，所以我们可以直接使用准确率来评估模型

print('ACC:%.3f'%(accuracy_score(y_test,y_pred)))
ACC:0.953

其实上面用的方法是留出法，我们也可以使用交叉验证法来计算模型误差。这样就把划分训练集和测试集、建立模型以及评估模型这几步合并在一起。

acc = np.mean(cross_val_score(LogisticRegression(),X,y,cv=10,scoring='accuracy'))
print('ACC:%.3f'%(acc))
ACC:0.946

可以看到两者比较接近。

四、优缺点

1.优点

（1）输出具有概率意义
（2）参数代表每个特征对输出的影响，可解释性强
（3）可以增加惩罚项来解决过拟合和多重共线性问题

2.缺点

（1）容易欠拟合（精度不高），所以最好结合其他模型一起使用
（2）本质上还是一个线性分类器，不好处理非线性问题