机器学习-最全LR逻辑回归原理以及python代码实现

简介

逻辑回归是一种分类模型，多用于二分类，以下讨论二分类情况。
逻辑回归是一种广义线性回归模型，它与传统线性回归的区别是，传统线性回归直接将$wx+b$作为因变量，而逻辑回归将$wx+b$用Sigmoid后激活的值作为因变量，可以将线性函数转化为概率：
$$P(y=1|x)=\frac{1}{1+exp[-(wx+b)]}(1)$$
即$$P(y=1|x)=\frac{exp(wx+b)}{1+exp(wx+b)}(2)$$
若将y视为事件发生的概率，1-y为反例发生的概率，这时，线性函数的值越接近正无穷，概率就越接近1，越接近负无穷，概率值就越接近0，这样的模型就是逻辑回归模型。
注意：为什么使用Sigmoid激活：1. 满足统计的最大熵模型
2. Sigmoid函数是平滑的，而且任意阶可导，一阶二阶导数可以直接由函数值得到不用进行求导，计算速度快，这在实现中很实用。

对于逻辑回归而言，输出Y=1的对数几率是X的线性模型：
$$log\frac{P(Y=1|x)}{1-P(Y=1|x)}=wx+b(3)$$

模型参数估计

由于逻辑回归模型假设样本服从伯努利分布，则可以写出模型的对数似然函数，然后应用极大似然估计求解参数。
假设$P(Y=1|x)=\pi (x)$，则$P(Y=0|x)=1-\pi (x)$
似然函数为：
$$\prod _{i=1}^N[\pi (x_i){]}^{y_i}[1-\pi (x_i){]}^{1-y_i}$$
对数似然函数为：
$$\begin{gathered} L(w,b)=\sum _{i=1}^N[y_{i}log\pi (x_i)+(1-y_i)log(1-\pi (x_i))] \\ =\sum _{i=1}^N[y_{i}log\frac{\pi (x_i)}{1-\pi (x_i)}+log(1-\pi (x_i))] \\ =\sum _{i=1}^N[y_{i}log\frac{\pi (x_i)}{1-\pi (x_i)}+log(1-\pi (x_i))](4) \end{gathered}$$
将（3）带入（4）得：
$$L(w,b)=\sum _{i=1}^N[y_i\ast (w{x}_i+b)-log(1+exp(w{x}_i+b))]$$
求出对数似然函数后，我们对其进行极大似然估计，及对$-L(w,b)$进行梯度下降，即可求出模型参数
有趣的是，我们对$L(w,b)$进行求导：
$$\begin{gathered} \frac{\partial L}{\partial w}=\sum _{i=1}^N[y_i\ast x_i-x_i\frac{exp(w{x}_i+b)}{(1+exp(w{x}_i+b))}] \\ =\sum _{i=1}^N{x}_i\ast [y_i-\pi (x_i)] \\ =-\sum _{i=1}^N{x}_i\ast error \end{gathered}$$
$$\begin{gathered} \frac{\partial L}{\partial b}=\sum _{i=1}^N[y_i-\frac{exp(w{x}_i+b)}{(1+exp(w{x}_i+b))}] \\ =\sum _{i=1}^N[y_i-\pi (x_i)] \\ =-\sum _{i=1}^{N}error \end{gathered}$$
这时会发现其实对数似然对$w$和$b$的导数只与$x$和误差$error=\pi (x)-y$有关，这也是逻辑回归使用sigmoid函数的原因之一，运算速度快，因为不需要进行求导运算。

补充：1. 逻辑回归为什么不用MSE作为损失函数

1. MSE 会有梯度消失现象
2. MSE的导数非凸函数，求解最优解困难

证明：
$$\delta (x)=sigmoid(x)=\frac{1}{1+e^{-wx}}$$
$$\begin{gathered} \frac{\partial \delta (x)}{\partial w}=-\frac{1}{(1+e^{-wx}{)}^2}\ast -x{e}^{-wx} \\ =x\delta (x)(1-\delta (x)) \end{gathered}$$
均方误差的损失函数为$L=(y-\delta (x){)}^2$
$$\begin{gathered} \frac{\partial L}{\partial w}=2(y-\delta (x))\frac{\partial \delta (x)}{\partial w} \\ =2x(\delta (x{)}^3-(1+y)\delta (x{)}^2+y\delta (x)) \end{gathered}$$

$$\frac{{\partial }^{2}L}{\partial w^2}=2x^2\delta (x)(1-\delta (x)(3\delta (x{)}^2-2(1+y)\delta (x)+y)$$
因为$\delta (x)\in (0,1)$，所以$2x^2\delta (x)(1-\delta (x)>0$
令$H(x)=3c^2-2(1+y)\delta (x)+y$
因为$y\in [0,1]$，

• y=0时，$H(x)=3\delta (x{)}^2-2\delta (x)=3\delta (x)(\delta (x)-\frac{2}{3})$；
  当$0<\delta (x)<\frac{2}{3}$时，$H(x)<0$
• y=1时，$H(x)=3\delta (x{)}^2-4\delta (x)+1=3(\delta (x)-1)(\delta (x)-\frac{1}{3})$；
  当$\frac{1}{3}<\delta (x)<1$时，$H(x)<0$

根据凸函数的定义：二阶导数非负，则L不为凸函数

2. LR为什么要将特征离散化

1. 离散化后的特征对异常数据有很强的鲁棒性：比如一个特征是年龄>30是1，否则0。如果特征没有离散化，一个异常数据“年龄300岁”会给模型造成很大的干扰。
2. 逻辑回归属于广义线性模型，表达能力受限；单变量离散化为N个后，每个变量有单独的权重，相当于为模型引入了非线性，能够提升模型表达能力，加大拟合。
3. 特征离散化以后，起到了简化了逻辑回归模型的作用，降低了模型过拟合的风险。(当使用连续特征时，一个特征对应于一个权重，那么，如果这个特征权重较大，模型就会很依赖于这个特征，这个特征的一个微小变化可能会导致最终结果产生很大的变化，这样子的模型很危险，当遇到新样本的时候很可能因为对这个特征过分敏感而得到错误的分类结果，也就是泛化能力差，容易过拟合）。而使用离散特征的时候，一个特征变成了多个，权重也变为多个，那么之前连续特征对模型的影响力就被分散弱化了，从而降低了过拟合的风险

---

python代码及详解

定义激活函数Sigmoid

from numpy import *
def sigmoid(x):
    return 1.0/(1+exp(-x))

训练，请注意观察梯度更新的计算，非常简便，只需要求误差，再乘以数据即可：

def train_logistic_regression(data, y):
    data, y = mat(data), mat(y).transpose()
    m, n = shape(data)
    weights = ones((n, 1))  # 初始化权重 n*1
    alpha = 0.001  # 学习率
    epoches = 10000  # 训练轮数
    l2_rate = 0.00001  # L2正则化系数

    for i in range(epoches):
        l2 = linalg.norm(weights, ord=2, keepdims=True)
        pred_y = sigmoid(data*weights)  # m*1
        error = y-pred_y  # m*1
        h = data.transpose()*error  # n*1
        weights -= alpha*h+l2*l2_rate
    return weights

预测

def test_logistic_regression(train_data, y, test_data):
    W = test_logistic_regression(train_data, y)
    pred_y=sigmoid(mat(test_data)*W)
    return pred_y

加载数据：这里使用sklearn自带cancer数据集

from sklearn import datasets
from sklearn.model_selection import train_test_split

def load_data():
    cancer = datasets.load_breast_cancer()
    X = cancer.data
    y = cancer.target
    X_train, X_test, y_train, y_test = train_test_split(X, y)
    return X_train, X_test, y_train, y_test

计算准确率

# 计算准确率
def accuracy(y,pred):
    N = len(y)
    res = sum([1 for i in range(N) if y_test[i] == pred[i]])
    return res/N

函数定义好之后，运行如下代码即可得到结果：

X_train, X_test, y_train, y_test = load_data()
pred = test_logistic_regression(X_train, y_train, X_test)
print('acc',accuracy(y_test,pred))