机器学习算法笔记(二十三)：逻辑回归初探

逻辑回归(Logistic Regression)，又被称为“对数几率回归”，是一个在各个行业中使用最多的机器学习算法，应用相当广泛。逻辑回归算法从线性回归算法衍生而来，虽然名叫“回归”，但其实是一种解决分类问题的算法。我们接下来就对逻辑回归的一些基础概念进行简单的探讨。

一、什么是逻辑回归

逻辑回归听名字好像是一个回归算法，但其实它解决的是分类问题。它的基本原理是将样本的特征和样本发生的概率联系起来，并预测样本发生的概率，然后根据概率进行分析。

由于预测得到的概率是一个数，所以是“回归”问题。

机器学习算法的本质是求出一个函数，如果有个样本 x 经过运算后就会得到一个预测值。之前的线性回归的本身就是我们关心的那个预测值（比如房价、成绩），而逻辑回归中，我们得到的值的本质是一个概率值，在这里记作。也就是说，对于逻辑回归来说我们是要得到一个，样本 x 经过运算后会计算出一个概率值，之后我们用这个概率值进行分类：

如果大于等于0.5，我们将设成1，反之设为0（1和0的定义在具体问题中具体看）。由此我们可以看出，逻辑回归既可以看做回归算法，也可以看做是分类算法，通常作为分类算法使用，只可以解决二分类问题。

● 如果模型只预测样本发生的概率，可称该算法为回归算法；如果预测出样本发生的概率后，再根据概率对样本做进一步确认，则可称该算法为分类算法。

● 可以使用一些其他技巧让逻辑回归解决多分类问题。与 kNN算法不同，kNN天生就能解决多分类问题。

下面我们来看一下逻辑回归如何得到一个事件发生概率的。之前的线性回归中，得到的的值域是(-∞, +∞)，可以求得一个任意的值，而对于概率来说，一个事件发生的概率只能在0和1之间，我们要做的就是将一个在(-∞, +∞)上的值“映射”到[0, 1]之间。再拿线性回归为例，我们依然是找一系列的 θ，与 X_b 点乘，我们再将这个结果作为变量，送给一个叫做 σ 的函数，转换成一个[0, 1]之间的值。

那这个 σ 到底是什么呢？我们通常取一个叫 Sigmoid 的函数，它的解析式为，图像大概是这样的：

这个函数能很好的将一个(-∞, +∞)内的值映射到(0, 1)。同时，当t>0时，σ(t)>0.5；当t<0时，σ(t)<0.5，这样的性质也能很好的反映上面分段的要求。将 θ 与 x_b 点乘的结果带入，得：

此时再注意，若如果大于等于0.5，我们将设成1，反之设为0，作为我们的分类结果。这里不妨举个例子：现在对于病人的数据，要预测这个病人患有良性肿瘤还是恶性肿瘤。首先通过训练得到一组 θ，这样每来一组新的病人的数据，我们就用这组数据前面加上一列“1”，与求得的 θ 点乘。点乘的结果得到的数再送给 Sigmoid 函数，得到的结果就是这个病人患有恶性肿瘤的概率。若概率大于等于0.5，我们就说这个病人很可能患有恶性肿瘤；若概率小于等于0.5，我们就说这个病人很可能是一个良性肿瘤的患者。

剩下的问题就变成了，我们给定一组样本数据集 X, y，如何找到参数 θ 使得可以最大程度获得样本数据集 X 对应的分类输出 y（也就是如何建模求出 θ）？

二、逻辑回归的损失函数

在线性回归中，θ 与 x_b 点乘的结果就是估计值，于是我们就可以用估计值减去真值的差来度量估计的好坏，差的平方和再进行平均(MSE)作为损失函数，我们只要找到使这个损失函数最小的 θ 就好了。同理，对于逻辑回归也要这么做，整体的方向是一样的，不过对于逻辑回归来说，找到它的损失函数相对没那么容易，接下来我们就一点一点往这个方向前进。

由于估计的结果有两类：1和0，相应的损失函数也可以分成这两类：

我们进一步用一个更为具体的函数来表示上面的趋势：

绘制这个损失函数的图像如下：

● y = 1（p ≥ 0.5）时，cost = -log(p)，p 越小，样本发生概率越小（最小为 0），则损失函数越大，分类预测值和实际值的偏差越大；相反，p 越大，样本发生概率越大（最大为 0.5），则损失函数越小，则预测值和实际值的偏差越小。

● y = 0（p ≤ 0.5）时，cost = -log(1-p)，p 越小，样本发生概率越小（最小为 0.5），则损失函数越大，分类预测值和实际值的偏差越大；相反，p 越大，样本发生概率越大（最大为 1），则损失函数越小，则预测值和实际值的偏差越小。

进一步将这个分段函数进行等价变换，写成一个函数：

假设一共有 m 个样本，相应的只要将这些损失加在一起即可：

又（其中X_b⁽ⁱ⁾表示全部样本），代入上式得：

这样我们就推导出了逻辑回归算法相应的损失函数。与线性回归不同，这个式子是没有解析解的，但我们仍然可以使用梯度下降法来求解出它的全局最优解。

三、对 J(θ) 的求导过程

由于博客不支持直接打出数学公式，所以下面用一系列截图来推导求导过程。

log(σ(t)) 函数的导数：

变形：

对 J(θ) 前半部分求导：

log(1 - σ(t)) 函数的导数：

变形：

对 J(θ) 后半部分求导：

前后相加：

最终化简得到的导数：

求梯度：

看到这个结果，我们发现，这个结果简直和线性回归太像了：，于是我们就能很方便的进行向量化操作。线性回归梯度的向量化结果为，类似地，就能很方便得出 J(θ) 梯度的向量化结果：

四、编程实现逻辑回归

在MyML包中新增一个 LogisticRegression.py 文件，实现如下代码：

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, 使用梯度下降法训练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"

        #根据推导的公式实现J和dJ
        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.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()"

新建一个main.py文件，实现如下代码：

import numpy as np
from sklearn import datasets

iris = datasets.load_iris()
X = iris.data
y = iris.target

#由于逻辑回归只能处理二分类问题，我们取第0和第1类鸢尾花
X = X[y<2,:2]
y = y[y<2]

from MyML.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(X, y, seed=666)

from MyML.LogisticRegression import LogisticRegression

log_reg = LogisticRegression()
log_reg.fit(X_train, y_train)

print(log_reg.score(X_test, y_test)) #打印预测准确率
print(log_reg.predict_proba(X_test)) #打印预测概率结果向量
print(log_reg.predict(X_test)) #打印分类结果向量
print(y_test) #打印测试数据集的标记

运行结果如下：