机器学习入门之逻辑回归——全网最详细（上）

本次博客主要内容是机器学习中非常著名的逻辑回归算法已经相关问题，分为上下两篇博客，并且会附上完整的数学推导和证明，希望大家在能成功地从坑里面爬出来！

一 .背景，起因与展望

（1）引入背景

在前面的机器学习中，我们绝大部分的工作是“预测”，即根据数据拟合模型，对于新进入的变量，根据前面的模型给出相应的预测值，然而，在日常生活中，我们并不需要一个准确的值，往往来说，对一件事的主观感觉才是我们关注点所在，这个时候，一个“可能性”，一个“概率”，也许比那些一大串具体的数字来的更直观有效，那么在机器学习中，解决这种问题的方法就是 “分类” ，所以，本次我们博客学习的主要问题就是分类算法，虽然它叫“逻辑回归”，但是还是逃脱不了它 “分类” 的本质！

（2）应用与适用

应用：

• 用于分类：适合做很多分类算法的基础组件。
• 用于预测：预测事件发生的概率（输出）。
• 用于分析：单一因素对某一个事件发生的影响因素分析（特征参数值）。

适用：

• 基本假设：输出类别服从伯努利二项分布。
• 样本线性可分。
• 特征空间不是很大的情况。
• 不必在意特征间相关性的情景。
• 后续会有大量新数据的情况。

（3）前景

这是来自著名的数据分析网站kaggle网站上 汇总的算法信息表，我们可以看出逻辑回归算法在数据分析领域使用的频率是最高的，虽然同样是分类问题，决策树的使用频率比逻辑回归少了将近20个百分点，在许多前沿科技论文中，使用逻辑回归算法进行模型拟合与数据分析也是非常的常见！

二 .初步理解与算法推导

无论是在机器学习上是刚入门不久的小可爱，还是资深的大佬，线性回归的问题已经老生常谈了，至少也应当是耳熟能详，那么仅与 “线性” 回归具有两字之差的逻辑回归是啥样的呢？和线性回归之间又有怎样的爱恨情仇呢？
$$\hat{y}=f(x)$$
在传统的线性回归或者是多线性回归中，我们会通过模拟数据集来训练模型从而得到$f(x)$，返回的$\hat{y}$ 是我们需要的直接答案，例如在房价预测模型中，返回的预测值就是未来房价，是个具体的值，而在逻辑回归中，却不是那么的简单了！
$$\hat{p}=f(x)$$
显然，这里的结果不再是 $\hat{y}$ 了，而是 $\hat{p}$ ，它代表的是一个 概率，但是我们最终需要的还是 $\hat{y}$，所以，我们定义了如下规则：
$$\hat{y}=\{\begin{array}{ll}1,& \hat{p}\ge 0.5\\ 0,& \hat{p}\le 0.5\end{array}$$

上面这些简单介绍，间接的告诉了我们逻辑回归可以看做是归类算法，也可以看做是分类算法，并且通常应用在分类算法中的二分类问题上！下面，我们一步一步的来介绍：当输入$x$时，模型是如何计算并得出概率值的！

在线性回归中，通常情况下，$\hat{p}=f(x)={\theta }^T\cdot x_b$，并且这个式子的值域是整个$\mathbb{R}$，算法需要拟合的就是这个${\theta }^T$，但这并不好，因为我们需要将所以的结果都收缩在$\left[0,1\right]$，单单这样一个简单的式子无法满足我们的需求，这个时候，需要的是一个“转换函数”，将结果转换到$\left[0,1\right]$ 之间，其类似于下面这种形式：
$$\hat{p}=\sigma \left({\theta }^T\cdot x_b\right)$$

根据前面所说，这里的$\sigma$是一个能“收缩”特定值 到固定区间内的神奇函数！我们一般把它称为:$Sigmoid$函数，其数学表达式为：
$$\sigma (t)=\frac{1}{1+e^{-t}}$$

也许你会疑惑为啥选择这个函数呢？看了下面的图你就知道了：

首先，这个函数的值域是$\left[0,1\right]$，并且在0附近骤增，当$t$足够大的时候，函数值趋近于1，当$t$ 比较小的时候，函数值趋近于0；问题来了，这有什么用，唉！巧妙之处就在这，如果这里的 $t$ 我们用 ${\theta }^T\cdot x_b$ 来替代的话，就恰好与前面概率$\hat{p}$ 对于$\hat{y}$的要求相匹配！其实原理上类似于这样的函数模型有很多，远远不止这一种，但是这里我们只需要理解这个即可，后面会有越来越多神奇的“激活函数”！

那么我们现在得到的就是：$$\hat{p}=\sigma \left({\theta }^T\cdot x_b\right)=\frac{1}{1+e^{-{\theta }^T\cdot x_b}}$$
再根据上面得到的 $\hat{p}$ 和 $\hat{y}=\{\begin{array}{ll}1,& \hat{p}\ge 0.5\\ 0,& \hat{p}\le 0.5\end{array}$ 我们就可以大致的得到一个二分类问题的结果！

到这里，方向有了，方法也逐渐浮现出来了，下一个急需解决的问题也就跃然纸上了，对于这个${\theta }^T$ 我们该怎么求呢？

三 . 损失函数推导

说一千道一万，评价我们的模型建立及其拟合效果才是王道，再把思绪拉回到我们学习过的线性回归中来，对于拟合后线性模型的准确率问题，我们沿用的一直都是高中就说过的最小二乘法，也就是MSE，虽然这里的研究对象换成了逻辑回归，但是整体的思路和方向没有发生特别大的改变，然而我们又不能像传统的MSE一样，因为逻辑回归等价于二分类问题，无论真实的准确率是多少，最后返回的结果只有 0 和 1 ，所以这里我们有需要引入一个新的方法来解决这个棘手的问题：

雄关漫道真如铁，而今迈步从头越! 解决问题的最好方法就是回到问题出现的原点，之前为了满足需要我们引入了$\hat{y}=\{\begin{array}{ll}1,& \hat{p}\ge 0.5\\ 0,& \hat{p}\le 0.5\end{array}$，也正是因为它，才导致我们的损失函数目前无从下手，所以焦点应当放在这里！

其实我们可以完全模仿这个分段函数来选择我们的 $loss$ $functions$，就像这样：
$$\text{ cost }=\{\begin{array}{l}\text{ 如果y=1, p越小, cost越大 }\\ \text{ 如果y=0, p越大, cost越大 }\end{array}$$

如果这个分类结果的 真实值 为1，那么返回的概率$p$越大，说明该模型鉴别或者说拟合的效果越好，如果得到的$p$越小，甚至趋近0.5，显然这个鉴别的准确度就比较低了，反之其损失函数的值也比较大，同样的道理也自然可以平移到真实值是 0 的时候！了解了这个规律之后，迫在眉睫的就是寻找到符合类似规律的数学函数！
$$\mathrm{cost}=\{\begin{array}{ccc}-\log (\hat{p})& \text{ if }& y=1\\ -\log (1-\hat{p})& \text{ if }& y=0\end{array}$$
我们来看一下这个cost function 对应的函数图像：

结合前面的解释来看这个公式，当 $y$ 的真实值为 1 的时候，随着$x$ 轴趋向于1，$-\log (\hat{p})$ 值也就越来越小，反正，$cost$ 的值就越来越大，道理是一样的，当 $y$ 的真实值为 0 的时，函数就是反过来的！

由于分段函数的不便，我们可以将上面的损失函数化简为：
$$\mathrm{cost}=-y\log (\hat{p})-(1-y)\log (1-\hat{p})$$

当然，这只是一个样本的的损失函数，当输入了m个样本数据之后，得到整体的损失函数为：
$$J(\theta )=-\frac{1}{m}\sum _{i=1}^m{y}^{(j)}\log \left({\hat{p}}^{(i)}\right)+\left(1-y^{(i)}\right)\log \left(1-{\hat{p}}^{(i)}\right)$$

这个时候，我们再把$\hat{p}=\sigma \left({\theta }^T\cdot x_b\right)$ 代入，就得到：

$$J(\theta )=-\frac{1}{m}\sum _{i=1}^m{y}^{(i)}\log \left(\sigma \left(X_b^{(i)}\theta \right)\right)+\left(1-y^{(i)}\right)\log \left(1-\sigma \left(X_b^{(i)}\theta \right)\right)$$

这里需要注意的是，$X_b^{(i)}$ 是数据$x_b$ 组成的矩阵，$m$ 行 $n$列 代表有 $m$ 个样本 ，每个样本有$n$ 个特征,所以$\hat{p}$ 的式子需要简单改一下！那么$\hat{p}$的结果就这个$X_b^{(i)}$ 的第 $i$ 行乘以$\theta$ 得到的结果再传给$\sigma$函数，从而得到概率！

再次确定我们的目的：通过大量的数据计算和模型拟合，最终在$cost$ $function$ 最小的时候求出$\theta$ 向量即可；但是非常遗憾，这个函数我们没有固定的公式套路去求解，需要借助前面学习过的梯度下降法！

我们通过原式对$\theta$ 求偏导，得到的最终结果如下：

$$\nabla J(\theta )=\frac{1}{m}\cdot \left(\begin{array}{c}{\sum }_{i=1}^m\left({\hat{y}}^{(i)}-y^{(i)}\right)\\ {\sum }_{i=1}^m\left({\hat{y}}^{(i)}-y^{(i)}\right)\cdot X_1^{(i)}\\ {\sum }_{i=1}^m\left({\hat{y}}^{(i)}-y^{(i)}\right)\cdot X_2^{(i)}\\ \dots \\ {\sum }_{i=1}^m\left({\hat{y}}^{(i)}-y^{(i)}\right)\cdot X_n^{(i)}\end{array}\right)=\frac{1}{m}\cdot X_b^T\cdot \left(\sigma \left(X_b\theta \right)-y\right)$$

这里就不展开推导了，毕竟我们是在学算法编程，数学推导的话就看小伙伴们高数下册学的怎么样了！

这里我们直接给出算法的源码：

import numpy as np
from metrics import accuracy_score

class LogisticRegression:

    def __init__(self):
        """初始化Logistic Regression模型"""
        self.coef_ = None
        self.intercept_ = None
        self._theta = None

    def _sigmoid(self,t):
        return 1./(1.+ np.exp(-t))
        

    def fit(self, X_train, y_train, eta=0.01, n_iters=1e4):
        """根据训练数据集X_train, y_train, 使用梯度下降法训练Linear Regression模型"""
        assert X_train.shape[0] == y_train.shape[0], \
            "the size of X_train must be equal to the size of y_train"

        def J(theta, X_b, y):
            y_hat = self._sigmoid(X_b.dot(theta))
            
            try:
                return -np.sum(y*np.log(y_hat)+(1-y)*np.log(1 - y_hat)) / len(y)
            except:
                return float('inf')

        def dJ(theta, X_b, y):

            return X_b.T.dot(self._sigmoid(X_b.dot(theta)) - y) / len(X_b)

        def gradient_descent(X_b, y, initial_theta, eta, n_iters=1e4, epsilon=1e-8):

            theta = initial_theta
            cur_iter = 0

            while cur_iter < n_iters:
                gradient = dJ(theta, X_b, y)
                last_theta = theta
                theta = theta - eta * gradient
                if (abs(J(theta, X_b, y) - J(last_theta, X_b, y)) < epsilon):
                    break

                cur_iter += 1

            return theta

        X_b = np.hstack([np.ones((len(X_train), 1)), X_train])
        initial_theta = np.zeros(X_b.shape[1])
        self._theta = gradient_descent(X_b, y_train, initial_theta, eta, n_iters)

        self.intercept_ = self._theta[0]
        self.coef_ = self._theta[1:]

        return self

    def predict_proba(self, X_predict):
        """给定待预测数据集X_predict，返回表示X_predict的结果概率向量"""
        assert self.intercept_ is not None and self.coef_ is not None, \
            "must fit before predict!"
        assert X_predict.shape[1] == len(self.coef_), \
            "the feature number of X_predict must be equal to X_train"

        X_b = np.hstack([np.ones((len(X_predict), 1)), X_predict])
        return self.sigmoid(X_b.dot(self._theta))

    def predict(self, X_predict):
        """给定待预测数据集X_predict，返回表示X_predict的结果向量"""
        assert self.intercept_ is not None and self.coef_ is not None, \
            "must fit before predict!"
        assert X_predict.shape[1] == len(self.coef_), \
            "the feature number of X_predict must be equal to X_train"

        proba = self.predict_proba(X_predict)
        return np.array(proba>=0.5,dtype = 'int')

    def score(self, X_test, y_test):
        """根据测试数据集 X_test 和 y_test 确定当前模型的准确度"""

        y_predict = self.predict(X_test)
        return accuracy_score(y_test, y_predict)

    def __repr__(self):
        return "LogisticRegression()"

这里的metrics文件也一并给大家：

import numpy as np
from math import sqrt


def accuracy_score(y_true, y_predict):
    """计算y_true和y_predict之间的准确率"""
    assert len(y_true) == len(y_predict), \
        "the size of y_true must be equal to the size of y_predict"

    return np.sum(y_true == y_predict) / len(y_true)


def mean_squared_error(y_true, y_predict):
    """计算y_true和y_predict之间的MSE"""
    assert len(y_true) == len(y_predict), \
        "the size of y_true must be equal to the size of y_predict"

    return np.sum((y_true - y_predict)**2) / len(y_true)


def root_mean_squared_error(y_true, y_predict):
    """计算y_true和y_predict之间的RMSE"""

    return sqrt(mean_squared_error(y_true, y_predict))


def mean_absolute_error(y_true, y_predict):
    """计算y_true和y_predict之间的RMSE"""
    assert len(y_true) == len(y_predict), \
        "the size of y_true must be equal to the size of y_predict"

    return np.sum(np.absolute(y_true - y_predict)) / len(y_true)


def r2_score(y_true, y_predict):
    """计算y_true和y_predict之间的R Square"""

    return 1 - mean_squared_error(y_true, y_predict)/np.var(y_true)

小伙伴们可以根据前面的数学推导解释，对照着代码一步一步的理解，其实并不是很难！

四 .逻辑回归 $VS$ 线性回归

最大区别在于它们的因变量不同，Logistic回归与现行回归属于同一家族：广义线性模型！

• 线性回归处理的是回归问题，逻辑回归处理的是分类问题，这是本质区别；
• 线性回归要求因变量是连续型的，Logistic回归要求因变量是分类变量（离散变量）；
• 线性回归要求变量服从正态分布，Logistic回归对变量分布无要求；
• 线性回归要求自变量和因变量呈线性关系，Logistic并不要求；
• Logistic回归是分析因变量取某个值的概率与自变量的关系，而线性回归是直接分析因变量与自变量的关系