机器学习 逻辑回归（Logistic Regression）

在解决分类问题时，如何用解决回归问题的方法来解决一个分类问题呢？那么逻辑回归会是一个很好的方法，逻辑回归将样本的特征和样本发生的概率联系起来，概率表示分类类别。

1 什么是逻辑回归

在线性回归中，我们得到的预测值形如：$$\hat{y}=f(x)$$得到的预测值可以作为回归结果输出，但逻辑回归的目的是为了分类，所以我们就要改造的我们函数：$$\hat{p}=f(x)$$这时我们得到的概率，再将概率映射到分类结果：$$\hat{y}=\{\begin{array}{c}1,\hat{p}\ge 0.5\\ 0,\hat{p}<0.5\end{array}$$这样我们就可以把利用回归方法得到概率转化到分类问题的结果上。所以逻辑回归既可以看做是回归算法，也可以看做分类算法，但从上面的式子可以看出逻辑回归作为分类算法使用时只可以解决二分类问题。

整体思路理解之后，那我们就着手处理第一个问题，我们模型中的预测函数$\hat{y}=f(x)$应该是什么样的，在线性回归中：$$\hat{y}=f(x)\Rightarrow \hat{y}={\theta }^T\cdot X_b$$这个函数的值域在正无穷负无穷之间（-infinity,+infinity），而概率的取值仅在0到1之间，这说明我们也要对上面的函数进行改造：$$\hat{p}=\sigma ({\theta }^T\cdot X_b)\sigma (t)=\frac{1}{1+e^{-t}}$$这里的函数$\sigma (t)$被称作Sigmoid函数，让我们把它画出来看一下这个函数的性质：

import numpy as np
import matplotlib.pyplot as plt

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

x = np.linspace(-10, 10, 500)
y = sigmoid(x)

plt.plot(x, y)
plt.show()

我们可以看到Sigmoid函数的值域为（0，1），当$t>0$时，$p>0.5$, 当$t<0$时，$p<0.5$，可见Sigmoid函数满足我们需要的所有条件（谁想出来的我不知道），于是我们就有了逻辑回归的模型：$$\hat{p}=\sigma ({\theta }^T\cdot X_b)=\frac{1}{1+e^{-{\theta }^T\cdot X_b}}$$$$\hat{y}=\{\begin{array}{c}1,\hat{p}\ge 0.5\\ 0,\hat{p}<0.5\end{array}$$有了模型需要的函数以后，我们接下来的问题就是对于给定的样本数据集$x,y$，我们如何找到参数$theta$，使得用上面这样的方式，可以最大程度获取样本数据集$x$对应的分类输出$y$，也就是像我上一章所讲到的一样，机器学习中参数学习都是要确定代价函数和效用函数，之后在使用不同方法去最优化代价函数与效用函数，那么我们现在的任务就是去确定代价函数。

2 逻辑回归的代价函数（损失函数）

对于分类问题，在一个样本分类错误的情况下，那么损失函数就会去加上一个固定的损失值，逻辑回归中，分类结果是由概率值确定的，也就是说我们有多大的把握确定它是否符合条件，形象的来说，一般简单的分类算法类似于k邻近算法的损失函数就如超市不小心把苹果装到梨的箱子中去卖了，这样放错了几个苹果超市只需要承担苹果和梨差价的几倍就可以了，但在逻辑回归中，如果我们是医院的医生，第一种情况模型告诉我们病人的肿瘤是良性的概率比百分之九十九还要高，我们就会告诉病人你的肿瘤不是恶性的，第二章情况如果模型告诉我们病人的肿瘤是良性的概率在百分之六十左右，我们就会告诉病人你的肿瘤很可能是良性的，不排除恶性的情况，最后两种情况下病人的肿瘤都是恶性的，我们所承受的代价也是不一样的，第二种情况下病人可能经过进一步检查发现了恶性肿瘤及时医治，而第一种情况下病人可能出院导致没有医治而丢掉性命，理解了这些我们就可以先去描述一下逻辑回归的代价函数：$$cost=\{\begin{array}{c}如果y=1,p越小,cost越大\\ 如果y=0,p越大,cost越大\end{array}$$结合上面的小故事理解一下，下面我们就是要找到一个适合去描述这样性质的函数，于是有了:$$cost=\{\begin{array}{cc}-log(\hat{p})，if& y=1\\ -log(1-\hat{p}),if& y=0\end{array}$$让我们画出这两个函数的图像来直观感受一下：

注意代价函数中的$p$也就是我们模型给出的概率估计，那么观察上面的两个函数，在$y=1$时函数$-log(p)$在（0，1）上是递减，也就是我们估计的概率越大，越确定$y=1$时代价就越小，如果当$p$很小，接近于0时，函数就会趋近正无穷，这样说明我们模型预测错误，代价函数会用很大的值去惩罚它，可以观察$y=0$时也是一样的。

了解了代价函数的机制以后，就要想办法去优化代价函数，也就是最小化代价函数，但上面这种形式显然不利于我们优化，于是我们对上面的函数进行改造，得到了下面这个简单的式子（我不知道谁想出来的，就是聪明）：$$cost=-ylog(\hat{p})-(1-y)log(1-\hat{p})$$这个简单的式子与上面的代价函数是等价的，要注意在实际代入每一个样本点的值时，由于y的取值不是0就是1，所以上面式子中是只有一项的，必定有一个系数为0，所以与之前的式子等价，体会一下。

于是我们就得到了我们的损失函数，从而把问题进一步化简：$$J(\theta )=-\frac{1}{m}\sum _{i=1}^m(y^{(i)}log(\sigma (X_b^{(i)}\theta ))+(1-y^i)log(1-\sigma (X_b^{(i)}\theta )))$$遗憾的是，上面这个式子不像线性回归一样，是没有正规方程解的，所有我们只能用梯度下降法来求解。

3 逻辑回归的梯度下降法

有了代价函数使用梯度下降法我们就要对$J(\theta )$求导，由于其中包含$\sigma (t)$，我们一步一步求导：$$\sigma (t)=\frac{1}{1+e^{-t}}{\sigma }^{\prime }(t)=-(1+e^{-t}{)}^{-2}\cdot e^{-t}\cdot (-1)=(1+e^{-t}{)}^{-2}\cdot e^{-t}$$$$\begin{gathered} (log\sigma (t){)}^{\prime }=\frac{1}{\sigma (t)}\cdot {\sigma }^{\prime }(t)=\frac{1}{\sigma (t)}\cdot (1+e^{-t}{)}^{-2}\cdot e^{-t} \\ =(1+e^{-t}{)}^{-1}\cdot e^{-t}=1-\frac{1}{1+e^{-t}}=1-\sigma (t) \end{gathered}$$$$\frac{d(y^{(i)}log\sigma (X_b^{(i)}{\theta }_j))}{d{\theta }_j}=y^{(i)}(1-\sigma (X_b^{(i)}\theta ))\cdot X_j^{(i)}$$$$\begin{gathered} (log(1-\sigma (t)){)}^{\prime }=\frac{1}{1-\sigma (t)}\cdot (-1)\cdot {\sigma }^{-1}(t)=-\frac{1}{1-{\sigma }^{-1}(t)}\cdot (1+e^{-t}{)}^{-2}\cdot e^{-t} \\ =-\frac{1+e^{-t}}{e^{-t}}\cdot (1+e^{-t}{)}^{-2}\cdot e^{-t}=-\sigma (t) \end{gathered}$$$$\frac{d[(1-y^{(i)})log(1-\sigma (X_b^{(i)}\theta ))]}{d{\theta }_j}=(1-y^{(i)})(-\sigma (X_b^{(i)}\theta ))\cdot X_j^{(i)}$$有了上面每一部分的导数，我们就可以$J(\theta )$求导：$$\frac{dJ(\theta )}{d{\theta }_j}=\frac{1}{m}\sum _{i=1}^m\sigma (X_b^{(i)}\theta -y^{(i)})X_j^{(i)}=\frac{1}{m}\sum _{i=1}^m({\hat{y}}^{(i)}-y^{(i)})X_j^{(i)}$$于是我们有了$J(\theta )$的梯度：$$▽J(\theta )=\left(\begin{array}{c}\alpha J/\alpha {\theta }_0\\ \alpha J/\alpha {\theta }_1\\ \alpha J/\alpha {\theta }_2\\ ...\\ \alpha J/\alpha {\theta }_n\end{array}\right)=\frac{1}{m}\left(\begin{array}{c}{\sum }_{i=1}^m({\hat{y}}^{(i)}-y^{(i)})\\ {\sum }_{i=1}^m({\hat{y}}^{(i)}-y^{(i)})X_1^{(i)}\\ {\sum }_{i=1}^m({\hat{y}}^{(i)}-y^{(i)})X_2^{(i)}\\ .....\\ {\sum }_{i=1}^m({\hat{y}}^{(i)}-y^{(i)})X_n^{(i)}\end{array}\right)=\frac{1}{m}\cdot X_b^T\cdot (\sigma (X_b\theta )-y)$$最后给出的结果是向量化后的结果，因为观察$J(\theta )$的梯度最后化简形式与线性回归模型中相似程度十分高，向量化步骤也是一样的，最后直接给出了向量化结果。

4 实现自己的逻辑回归

import numpy as np
from .metrics import accuracy_score


class LogisticRegression:

    def __init__(self):
        """初始化Logistic Regression模型"""
        self.coef_ = None
        self.interception_ = 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,使用批量梯度下降法训练Logistic 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(y)

        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.interception_ = self._theta[0]
        self.coef_ = self._theta[1:]

        return self


    def predict(self, X_predict):
        """给定待预测数据集X_predict, 返回表示X_predict的结果向量"""
        assert self.interception_ 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 predict_proba(self, X_predict):
        """给定待预测数据集X_predict, 返回表示X_predict的结果概率向量"""
        assert self.interception_ 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 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()"

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

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

下面使用一下我们自己的逻辑回归去应用于真实的鸢尾花数据集：

这里为了可视化方便，我只取了两个特征，同时逻辑回归是一个二分类算法，所以我只取了其中的两种鸢尾花：

使用我们自己的逻辑回归：

这里可以以看到我们模型预测的所有样本都正确了，并且可以从结果概率向量看出来逻辑回归分类的原理，同时这里去两个特征为了可视化是方便引出对于分类算法来说，每一个分类算法都有一个决策边界。

5 决策边界

我们回过头来再看逻辑回归：$$\hat{p}=\sigma ({\theta }^T,X_b)=\frac{1}{1+e^{-{\theta }^T{X}_b}}$$对于上式来说，我们每次新来一个样本，样本就会与${\theta }^T$去乘积最后得到概率，简而言之：$$y=\{\begin{array}{c}1,\hat{p}\ge 0.5\Rightarrow {\theta }^T\cdot X_b\ge 0\\ 0,\hat{p}<0.5\Rightarrow {\theta }^T\cdot X_b<0\end{array}$$这样我们就可以看出逻辑回归的决策边界为：$${\theta }^T\cdot X_b=0$$如果X有两个特征：$${\theta }_0+{\theta }_1{X}_1+{\theta }_2{X}_2=0$$代表了一条直线，那么：$$x_2=\frac{-{\theta }_0-{\theta }_1{x}_1}{{\theta }_2}$$下面让我们用代码来看一下决策边界：

我们按照上面的逻辑画出了我们逻辑回归的决策边界，由于我们之前取的两个特征，就是为了这里的可视化，可以看到上面有一个点是分类错误的，可能你还记得上一节最后的分类评价分数是满分，原因是在测试集上的分类是完全正确的：

这里的决策边界都是一条直线，如果决策边界是不规则的情况怎么办呢？我们最笨的思路是将整个平面分成足够多的点，计算每个点是什么类别，然后给它图上颜色，这样只要点足够密集就可以看出决策边界，按照我们的思路：

def plot_decision_boundary(model, axis):
    
    x0, x1 = np.meshgrid(
        np.linspace(axis[0], axis[1], int((axis[1]-axis[0])*100)).reshape(-1, 1),
        np.linspace(axis[2], axis[3], int((axis[3]-axis[2])*100)).reshape(-1, 1),
    )
    X_new = np.c_[x0.ravel(), x1.ravel()]

    y_predict = model.predict(X_new)
    zz = y_predict.reshape(x0.shape)

    from matplotlib.colors import ListedColormap
    custom_cmap = ListedColormap(['#EF9A9A','#FFF59D','#90CAF9'])
    
    plt.contourf(x0, x1, zz, linewidth=5, cmap=custom_cmap)

我们作用于逻辑回归试一下：

为了得到不规则的边界我这里下面选用K邻近算法，knn算法是一个理论度很低的算法，基本没有用到什么数学推导，就是一个根据距离投票的分类算法，让我们来直观的看一下knn算法的决策边界：

这里我们可以看出来在欧式距离下的knn算法二分类的决策边界就明显是一个不规则的决策边界了，让我们试一下多分类：

在多分类情况下我们使用sklearn中的默认参数使得我们的模型决策边界变得特别不规则，实质上这种极度不规则的情况是过拟合的表现，过拟合是模型过度学习训练数据集样本的属性特征，导致在训练数据集上表现极好，但在实际应用表现欠佳，让我们来调整一下超参数：

这样我们把knn的超参数默认值的5变为50，我们就得到了一个比较清晰的决策边界。

6 在逻辑回归中使用多项式回归

上面我们学到的逻辑回归，就是找到一个直线或者超平面去分割样本点，但是缺陷就是只能用于二分类，所以我们将多项式回归用于逻辑回归，尝试得到非直线的决策边界。多项式回归完全应用的是线性回归的思想，但我们给样本添加若干特征，这些特征时样本已有特征的线性组合。

这里我们模拟创建二次曲线围成的决策边界。

我们首先使用之前我们自己实现的逻辑回归去拟合一下新的样本数据。

可见以直线为决策边界的逻辑回归并不能很好完成分类，结果的评分也只有0.6分。

上面我们引入sklearn中的管道去封装一下我们的带有多项式类的逻辑回归，这里要注意管道中的第三项是我们自己实现的逻辑回归而不是sklearn中自带的，没有报错是因为我们完全按照sklearn的规则去实现的逻辑回归，这也是我们爱编写代码时所规定的几个函数来的好处，再我们引入二次项后，我们的模型的评分可以达到0.95，已经是很不错的结果了，让我们来看它的决策边界：

但在实际数据中，我们无法确定我们的数据中到底包含几次项之间的关系，所有我们会加入不同的次方项来拟合数据，下面我们使用一个较高的次项：

可见这时候模型的决策边界在边缘地方已经不是很规则了，这就是过拟合的变现，所以我们要对多项式回归使用正则化。

7 scikt-learn中的逻辑回归

在讲述sklearn中的逻辑回归之前，我们要先引入一个概念，这个概念就是模型正则化（Regularization），模型正则化的作用就是限制参数的大小，在代价函数中我们为了去限制参数的大小会在代价函数之后加上一项正则项，其中会分为L1正则和L2正则：$$Ridge\sum _{i=1}^n{\theta }_i^2{L}_2正则项$$$$LASSO\sum _{i=1}^n|{\theta }_i|L_1正则项$$关于这两个正则项会有不同的限制参数大小的作用，感兴趣的同学可以自己搜索查询一下二者的差别。

sklearn中的使用正则化的方式为：$$C\cdot J(\theta )+L_1或L_2$$下面我们看一下sklearn中如何封装逻辑回归和使用正则化：

这里创建数据去模拟二次曲线，为了模拟正则化作用，我随机选择了20个点将其样本输出置为1，这样相当于手动为我们的数据添加了噪音。

我们使用sklearn中的逻辑回归，可以看到第一个参数C就是上面我们的$J(\theta )$前面的系数C，penalty为正则项，sklearn中默认使用为L2正则。

可见不加上多项式项的逻辑回归预测的效率也不是很高，sklearn中封装的逻辑回归仅有0.8分，同时决策边界为直线。

我们引出多项式项后，明显模型预测评分达到了0.94，让我们来看一下决策边界：

这里画出的决策边界已经是很不错的边界了，让我们继续用高次项去拟合数据，这样再去看正则项的作用，下面是20次方的模型画出的决策边界：

我们去重新封装管道，去调节我们超参数C去看一下效果：

可以看到我们把C变为0.1之后，实质上是正则项所占的比例更大，虽然得到的决策边界并不是特别准确，还是存在过拟合，但是相比上一个决策边界已经有很好的缓和了。

调节完超参数C接着我们去改变正则项尝试一下是否能再次优化：

当我们使用L1正则项时，可以看到模型的决策边界已经非常相似于我们创建样本数据时的二次曲线了。

在实际建立模型时，我们很多时候是不知道我们的超参数应该取什么值的，这时我们就要再引入一个概念了，交叉验证（Cross Validation），下面我要偷一点懒了，放一下笔记把：

交叉验证的步骤实质上就是将数据分为训练集，测试集还有验证数据集。训练集调试参数，验证集调试超参数，测试集衡量我们的模型，同时将其中每个部分交换，相互作为其他的部分而选出最好的参数与超参数。如果还是对交叉验证这部分感触不深可以看周志华的机器学习西瓜书，上面会有更详细的例子。

8 分类算法的评价

在描述分类算法评价之前，我想先介绍一种整体模型评估我们所选择指标，在我学习西瓜书以后大致应该可以给出下面的思路：

上面的方法不一定准确，但是希望在之后的学习中可以帮助到大家。下面我们介绍分类算法的评价，首先我们要从分类准确度问题来引入：这里我们假设一下，我说我可以很轻松的开发一个系统，系统的准确度可高达99.9%，你觉得这种可能性大吗？实际上是很容易的事情，假设我开发一个癌症预测系统，输入体检信息，可以判断是否患有癌症，假设这种癌症的产生率在0.01%，那我只需要直接预测所有人都健康就可以达到99.9%的准确率。怎么样？是不是觉得有道理，但是又觉得这个模型有什么地方不对劲？

原因就出在这个系统的拿到的是极度偏斜（skewed Data）的数据，只使用分类准确度是远远不够的。 这时我们就要请出我们的新工具进行进一步分析：混淆矩阵（Confusion Matrix）

这里的精准率和召回率也被称为查准率和查全率，精准率和召回率是一对矛盾的度量，精准率高时召回率往往偏低，召回率高时精准率往往偏低。被称为P-R反向变动关系：

为了二者都兼顾，我们引入F1 Score：

其实到这里博主也是初学者，仅仅在理解混淆矩阵与掌握性质的情况下并不能很清楚的讲明中间的数学建模，同时ROC曲线和代价敏感错误率曲线都是及其重要的指标，但是如果用通俗的语言来形容就无法描述曲线所代表的数学意义了，建议大家这里一定要去看书或者搜集资料，下面我只放上两种曲线的定义：





书上式子虽然枯燥晦涩，但是这里数学公式是最好学习评价机器学习模型的工具，上面尤其要仔细理解AUC和代价敏感错误率，这里仅凭文字很难理解，要结合资料并且去搜索课程完整的理解上面几个曲线这意义，这并不是仅仅能去评价模型结果，我认为对于整个机器学习的思路和整体把握都由很大的帮助，到这里基本上所有机器学习中的一些概念与评价指标都介绍清楚了，之后就为我们接触更加复杂的算法奠定了基础，大家加油！