机器学习算法整理(三)

原创
2021/09/10 16:37
阅读数 3.7K

机器学习算法整理(二)

逻辑回归

什么是逻辑回归(Logistic Regression)

逻辑回归是解决分类问题的,那回归问题怎么解决分类问题呢?将样本的特征和样本发生的概率联系起来,概率是一个数。

对于机器学习的本质就是,进来一个x,经过f(x)的运算,就得到一个预测值,对于之前无论是线性回归也好,多项式回归也好,看我们要预测的是什么,如果我们要预测的是房价,那么的值就是房价。如果我们要预测成绩,那么的值就是成绩。但是在逻辑回归中,这个的值是一个概率值。

,进来一个x,经过f(x)的运算,会得到一个概率值。之后我们根据这个概率值来进行分类

如果有50%以上的概率,我们就让的值为1;如果在50%以下概率的话,我们就让的值为0。这里的1和0在实际的场景中可能代表不同的意思,比如说1代表恶性的肿瘤患者,0代表良性的肿瘤患者;或者说1代表银行发给某人信用卡有一定的风险,0代表没有风险等等。

逻辑回归既可以看作是回归算法,也可以看作是分类算法。如果我们不进行最后的一步根据的值进行分类的操作,那么它就是一个回归算法。我们计算的是根据样本的特征来拟合计算出一个事件发生的概率。比如给我一个病人的信息,我计算出他患有恶性肿瘤的概率。给我一个客户的信息,我计算出发给他信用卡产生风险的概率。我们根据这个概率进一步就可以进行分类。不过通常我们使用逻辑回归还是当作分类算法用,只可以解决二分类问题。如果对于多分类问题,逻辑回归本身是不支持的。当然我们可以使用一些其他的技巧进行改进,使得我们用逻辑回归的方法,也可以解决多分类的问题。但是对于KNN算法来说,它天生就可以支持多分类的问题。

逻辑回归使用一种什么方式可以得到一个事件概率的值?对于线性回归来说,它的值域是(-∞,+∞)的。对于线性回归来说它可以求得一个任意的值。但是对于概率来说,它的值域只能是[0,1],所以我们直接使用线性回归的方式,没办法在这个值域内。为此,我们的解决方案就是将线性回归的预测值,放入一个函数内,使得其值域在[0,1]之间,这样就获得了我们需要的概率

那么我们的函数是什么呢?

现在我们用代码来看一下这个函数的图像

import numpy as np
import matplotlib.pyplot as plt

if __name__ == "__main__":

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

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

运行结果

根据这根曲线,我们可以看到它的最左端趋近于0,但是达不到0,最右端趋近于1,但是达不到1。说明这根曲线的值域是(0,1),这是因为=0,=1。当我们画上纵轴

我们可以看到,当t>0时,p>0.5;t<0时,p<0.5。

那么将我们线性回归的预测值代入该函数后就变成了

现在我们已经知道了逻辑回归的基本原理,现在的问题就是:对于给定的样本数据集X,y,我们如何找到参数θ,使得用这样的方式,可以最大程度获得样本数据集X,对应分类出y?

逻辑回归的损失函数

我们先来看一下线性回归的损失函数,而是真值。我们只需要找到让这个损失函数最小的θ值就好了。对于逻辑回归来说,它的预测值要么是1,要么是0,而我们是根据估计出来的,它代表一个概率,来决定到底是1还是0,它分成了两类,相应的逻辑回归的损失函数也相应的分成了两类。如果真值y=1,那么p越小,我们估计出来的就越有可能是0,那么我们的损失(cost)就越大;如果真值y=0,那么p越大,我们估计出来的就越有可能是1,那么我们的损失(cost)就越大。那么损失的趋势为

对于这种趋势,我们可以使用这样的函数来表示

我们把损失函数定义成这两种情况,我们先来看一下y=log(x)的函数图像

而我们的第一个式子中是一个-log(x),那么-log(x)的图像就为

由于估计值的值域为[0,1],所以上面这个图像超过1的部分都是没有意义的。

通过这个图,我们可以发现当取0的时候,-log()趋近于+∞。按照我们之前分类的方式,当=0的时候,预测值应该为0,但样本真实值y实际是1,显然我们分错了,此时我们对其进行惩罚,这个惩罚是+∞的。随着逐渐的增大,我们的损失越来越小,当到达1的时候,根据分类标准,预测值为1,此时它跟样本真值y是一致的。此时-log()为0,也就是没有任何损失。这样一根曲线就描述了如果我们给了一个样本,这个样本真实的标记输出是1的时候,相应用我们的方法估计出了以后,代入-log()得到的损失函数

由于y=-log(x)是上图的样子,那么y=-log(-x)就是根据y轴对称的样子

那么对于-log(1-x)就是将-log(-x)的曲线向右平移一个单位。

由于的值域为[0,1],所以y=-log(1-x)小于0的部分是没有意义的。

在这根曲线上,如果给定的=1的话,那么此时是趋近于+∞的。当=1的时候,预测值=1,但是y的真值为0,我们完全分错了,所以我们给它一个+∞的惩罚,随着的逐渐减小,这个惩罚值会越来越低,直到当=0的时候,=0,而y的真值为0,所以此时分类正确,在这种情况下一点惩罚都没有。所以我们用这两根曲线来作为我们的损失函数。

由于这样写非常不方便,所以我们把它合并成一个函数

,这个函数和上面的分类函数是等价的,原因也很简单,当y=1的时候,该式就等于,当y=0的时候,该式就等于

当我们面对一个样本X的时候,相应的我们也知道这个样本对应的真实的分类y。现在我们使用逻辑回归的方式就可以估计出来对于这个样本X,它相应的概率是多少,于是我们就用这个式子得到了它相应的损失是多少。我们会来m个样本,相应的我们只需要将这些损失加在一起就可以了。

这就是我们逻辑回归相应的损失函数。

由于

所以最终我们的损失函数就为

这里我们要给定样本相应的特征,用表示,以及对应的输出标记,用表示,这些都是已知的。我们真正要求的是里面的θ。之后我们要做的事情就是找到一组θ,使得我们的J(θ)达到一个最小值。对于这个式子,它不能像线性回归使用最小二乘法一样推导出一个正规方程解,它是没有数学的解析解的。但是它可以使用梯度下降法来求解。这个损失函数是一个凸函数,它是没有局部最优解的,只存在唯一的一个全局最优解。

逻辑回归损失函数的梯度

对于这个损失函数

要求它的梯度,其实就是对每一个θ求偏导数

我们先来对函数求导

这是一个复合函数,根据复合函数求导法则,令a=-t,b=e^-t,c=,则a'=-1,b'=e^-t,c'=

所以

我们再来看一下的导数。这也是一个复合函数。log'x=1/x,则

那么对于损失函数前半部分的导数就为(这里同样也是复合函数求导,跟随的特征就是)

对于后半部分,我们先看一下的导数,它同样是一个复合函数,a=logx,b=1-x,c=,则a'=1/x,b'=-1,则

由于,代入上式,可得

的导数就为(这里同样也是复合函数求导,跟随的特征就是)

由于前半部分的导数=

后半部分的导数=

合并后就为

最终J(θ)对某一个θ求偏导的话就为(前面的负号放到了括号中)

就是我们逻辑回归的预测值,所以上式可以写为

则J(θ)的梯度就为

我们再来对比一下线性回归的梯度

我们会发现,逻辑回归和线性回归的梯度其实是一致的,只不过线性回归的预测值=;而逻辑回归的预测值=。线性回归是均方误差,所以前面多了一个2,而逻辑回归没有这个平方,所以前面没有这个2。对于线性回归进行向量化处理,它的梯度可以写成

那么对逻辑回归的梯度进行向量化处理,就有

实现逻辑回归算法

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之间的MAE"""

    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)
import numpy as np
from .metrics import accuracy_score


class LogisticRegression:

    def __init__(self):
        # 初始化LogisticRegression模型
        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], \
            "X_train的列数必须等于y_train的长度"

        def J(theta, X_b, y):
            # 构建损失函数
            p_hat = self._sigmoid(X_b.dot(theta))
            try:
                return - np.sum(y * np.log(p_hat) + (1 - y) * np.log(1 - p_hat)) / len(y)
            except:
                return float('inf')

        def dJ(theta, X_b, y):
            # 对theta求偏导数,获取梯度向量
            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):
            """
            梯度下降算法
            :param X_b: 带虚拟特征1的自变量特征矩阵
            :param y: 因变量向量
            :param initial_theta: 初始的常数向量,这里需要注意的是真正待求的是常数向量,求偏导的也是常数向量
            :param eta: 迭代步长、学习率
            :param n_iters: 最大迭代次数
            :param epsilon: 误差值
            :return:
            """
            theta = initial_theta
            # 真实迭代次数
            i_iter = 0
            while i_iter < n_iters:
                # 获取梯度
                gradient = dJ(theta, X_b, y)
                last_theta = theta
                # 迭代更新theta,不断顺着梯度方向寻找新的theta
                theta = theta - eta * gradient
                # 计算前后两次迭代后的损失函数差值的绝对值
                if abs(J(theta, X_b, y) - J(last_theta, X_b, y)) < epsilon:
                    break
                # 更新迭代次数
                i_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_proba(self, X_predict):
        # 给定待预测数据集X_predict,返回表示X_predict的结果概率向量
        assert self.interception is not None and self.coef is not None, \
            "开始预测前必须fit"
        assert X_predict.shape[1] == len(self.coef), \
            "预测的特征数必须与训练的特征数相等"
        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.interception is not None and self.coef is not None, \
            "开始预测前必须fit"
        assert X_predict.shape[1] == len(self.coef), \
            "预测的特征数必须与训练的特征数相等"
        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()"

现在我们来使用逻辑回归判别鸢尾花数据集,先获取鸢尾花数据集

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets

if __name__ == "__main__":

    iris = datasets.load_iris()
    X = iris.data
    y = iris.target
    # 鸢尾花数据集有3个分类,而逻辑回归只能进行二分类,所以我们去掉2这个特征
    X = X[y < 2, :2]
    y = y[y < 2]
    print(X.shape)
    print(y.shape)
    plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
    plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    plt.show()

运行结果

(100, 2)
(100,)

此时我们看到鸢尾花的数据集有100个样本数,2个特征。现在我们就用我们自己写的逻辑回归进行分类

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from playLA.model_selection import train_test_split
from playLA.LogisticRegression import LogisticRegression

if __name__ == "__main__":

    iris = datasets.load_iris()
    X = iris.data
    y = iris.target
    # 鸢尾花数据集有3个分类,而逻辑回归只能进行二分类,所以我们去掉2这个特征
    X = X[y < 2, :2]
    y = y[y < 2]
    print(X.shape)
    print(y.shape)
    # plt.scatter(X[y == 0, 0], X[y == 1, 1], color='red')
    # plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    # plt.show()

    X_train, X_test, y_train, y_test = train_test_split(X, y, seed=666)
    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(y_test)
    # 查看预测结果
    print(log_reg.predict(X_test))

运行结果

(100, 2)
(100,)
1.0
[0.92972035 0.98664939 0.14852024 0.01685947 0.0369836  0.0186637
 0.04936918 0.99669244 0.97993941 0.74524655 0.04473194 0.00339285
 0.26131273 0.0369836  0.84192923 0.79892262 0.82890209 0.32358166
 0.06535323 0.20735334]
[1 1 0 0 0 0 0 1 1 1 0 0 0 0 1 1 1 0 0 0]
[1 1 0 0 0 0 0 1 1 1 0 0 0 0 1 1 1 0 0 0]

通过训练和查看分类准确度,我们可以看到,我们对鸢尾花数据集的测试数据全部都正确的进行了分类。并且可以查看测试数据集元素概率值。并通过概率值获取最终的分类结果和预测结果。

决策边界

对于逻辑回归,它同样也有系数和截距,我们来看一下鸢尾花数据集的系数和截距

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from playLA.model_selection import train_test_split
from playLA.LogisticRegression import LogisticRegression

if __name__ == "__main__":

    iris = datasets.load_iris()
    X = iris.data
    y = iris.target
    # 鸢尾花数据集有3个分类,而逻辑回归只能进行二分类,所以我们去掉2这个特征
    X = X[y < 2, :2]
    y = y[y < 2]
    print(X.shape)
    print(y.shape)
    # plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
    # plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    # plt.show()

    X_train, X_test, y_train, y_test = train_test_split(X, y, seed=666)
    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(y_test)
    # 查看预测结果
    print(log_reg.predict(X_test))
    # 查看逻辑回归的系数
    print(log_reg.coef)
    # 查看逻辑回归的截距
    print(log_reg.interception)

运行结果

(100, 2)
(100,)
1.0
[0.92972035 0.98664939 0.14852024 0.01685947 0.0369836  0.0186637
 0.04936918 0.99669244 0.97993941 0.74524655 0.04473194 0.00339285
 0.26131273 0.0369836  0.84192923 0.79892262 0.82890209 0.32358166
 0.06535323 0.20735334]
[1 1 0 0 0 0 0 1 1 1 0 0 0 0 1 1 1 0 0 0]
[1 1 0 0 0 0 0 1 1 1 0 0 0 0 1 1 1 0 0 0]
[ 3.01796521 -5.04447145]
-0.6937719272911228

这些系数和截距就组成了我们之前一直说的θ值。对于θ值它有没有几何意义呢,我们应该如何看待这个θ值呢?

之前我们在说逻辑回归原理的时候,有这样的式子

我们通过训练,得出了这个θ向量,如果我们的样本有n个特征(维度)的话,那么θ就有n+1个元素,我们要多添加一个。每当新来一个样本的时候,样本与θ进行点乘,点乘之后的结果送给函数,得到的这个值,我们就称之为是这个样本发生我们定义的事件的概率值。如果概率值≥0.5的话,我们就分类这样样本属于1这一类;如果概率值<0.5的话,我们就分类这个样本属于0这一类。我们再来看一下函数

当t>0时,p>0.5;t<0时,p<0.5,分界点在t=0的时候。将代入t,则有

≥0的时候,我们的概率估计值≥0.5,此时我们将新来的样本x分类为y=1;当<0的时候,我们的概率估计值<0.5,此时我们将新来的样本x分类为y=0。换句话说,我们为新来的样本分类为1还是0,这个边界点,这个位置就被称为决策边界

两个向量进行点乘,它同时代表的是一条直线。如果X有两个特征,那么就可以写成

这样的一个式子是一个直线的表达式。在二维坐标系中,横轴是x1这个特征,纵轴是x2这个特征。经过转换,就有

推导出x2,我们是为了把这根直线给画出来,真实感性的看一下它的样子。

def x2(x1):
    return (- log_reg.coef[0] * x1 - log_reg.interception) / log_reg.coef[1]

x1_plot = np.linspace(4, 8, 1000)
x2_plot = x2(x1_plot)
plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
plt.plot(x1_plot, x2_plot)
plt.show()

运行结果

图中的这根直线就是我们的决策边界。这根决策边界大体上将红色点和蓝色的点分成了两部分。如果我们新来一个样本的话,把样本的每一个特征和θ相乘,如果大于等于0的话,我们就给它分类为1,就在这根直线的下面。这就是这根直线的几何意义。如果得到的结果是小于0的话,我们就给它分类为0,就在这根直线的上方。如果等于0的话,它正好落在这一根决策边界上,此时无论我们把这个样本分类成0还是分类成1都是正确的,只不过我们这里分类成1。这里我们会发现在上图中有一个红色的点在决策边界的下方,显然这是一个分类错误的情况。而我们之前在测试用例中的分类准确率是100%,说明这个红色的点是在训练数据集中,我们可以来验证一下

def x2(x1):
    return (- log_reg.coef[0] * x1 - log_reg.interception) / log_reg.coef[1]

x1_plot = np.linspace(4, 8, 1000)
x2_plot = x2(x1_plot)
# plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
# plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
# plt.plot(x1_plot, x2_plot)
# plt.show()

plt.scatter(X_test[y_test == 0, 0], X_test[y_test == 0, 1], color='red')
plt.scatter(X_test[y_test == 1, 0], X_test[y_test == 1, 1], color='blue')
plt.plot(x1_plot, x2_plot)
plt.show()

运行结果

由上图可见,我们的测试数据集只有这么几个点,它完全是分类正确的。

不规则的决策边界的绘制方法

如果我们的分类样本本身没有一条这样的泾渭分明的可以用数学表达式来表达的决策边界

那我们如何来画出这种不规则的决策边界呢?我们依然还是用之前的逻辑回归的样本数据来进行决策边界的绘制

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from playLA.model_selection import train_test_split
from playLA.LogisticRegression import LogisticRegression
from matplotlib.colors import ListedColormap

if __name__ == "__main__":

    iris = datasets.load_iris()
    X = iris.data
    y = iris.target
    # 鸢尾花数据集有3个分类,而逻辑回归只能进行二分类,所以我们去掉2这个特征
    X = X[y < 2, :2]
    y = y[y < 2]
    # print(X.shape)
    # print(y.shape)
    # plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
    # plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    # plt.show()

    X_train, X_test, y_train, y_test = train_test_split(X, y, seed=666)
    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(y_test)
    # # 查看预测结果
    # print(log_reg.predict(X_test))
    # # 查看逻辑回归的系数
    # print(log_reg.coef)
    # # 查看逻辑回归的截距
    # print(log_reg.interception)

    # def x2(x1):
    #     return (- log_reg.coef[0] * x1 - log_reg.interception) / log_reg.coef[1]
    #
    # x1_plot = np.linspace(4, 8, 1000)
    # x2_plot = x2(x1_plot)
    # plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
    # plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    # plt.plot(x1_plot, x2_plot)
    # plt.show()

    # plt.scatter(X_test[y_test == 0, 0], X_test[y_test == 0, 1], color='red')
    # plt.scatter(X_test[y_test == 1, 0], X_test[y_test == 1, 1], color='blue')
    # plt.plot(x1_plot, x2_plot)
    # plt.show()

    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)
        custom_cmap = ListedColormap(['#EF9A9A', '#FFF59D', '#90CAF9'])
        plt.contourf(x0, x1, zz, linewidth=5, cmap=custom_cmap)

    plot_decision_boundary(log_reg, axis=[4, 7.5, 1.5, 4.5])
    plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
    plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    plt.show()

运行结果

现在我们来看一下KNN算法的决策边界,KNN算法的决策边界也是没有数学表达式的

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from playLA.model_selection import train_test_split
from playLA.LogisticRegression import LogisticRegression
from matplotlib.colors import ListedColormap
from sklearn.neighbors import KNeighborsClassifier

if __name__ == "__main__":

    iris = datasets.load_iris()
    X = iris.data
    y = iris.target
    # 鸢尾花数据集有3个分类,而逻辑回归只能进行二分类,所以我们去掉2这个特征
    X = X[y < 2, :2]
    y = y[y < 2]
    # print(X.shape)
    # print(y.shape)
    # plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
    # plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    # plt.show()

    X_train, X_test, y_train, y_test = train_test_split(X, y, seed=666)
    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(y_test)
    # # 查看预测结果
    # print(log_reg.predict(X_test))
    # # 查看逻辑回归的系数
    # print(log_reg.coef)
    # # 查看逻辑回归的截距
    # print(log_reg.interception)

    # def x2(x1):
    #     return (- log_reg.coef[0] * x1 - log_reg.interception) / log_reg.coef[1]
    #
    # x1_plot = np.linspace(4, 8, 1000)
    # x2_plot = x2(x1_plot)
    # plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
    # plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    # plt.plot(x1_plot, x2_plot)
    # plt.show()

    # plt.scatter(X_test[y_test == 0, 0], X_test[y_test == 0, 1], color='red')
    # plt.scatter(X_test[y_test == 1, 0], X_test[y_test == 1, 1], color='blue')
    # plt.plot(x1_plot, x2_plot)
    # plt.show()

    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)
        custom_cmap = ListedColormap(['#EF9A9A', '#FFF59D', '#90CAF9'])
        plt.contourf(x0, x1, zz, linewidth=5, cmap=custom_cmap)

    # plot_decision_boundary(log_reg, axis=[4, 7.5, 1.5, 4.5])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
    # plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    # plt.show()

    # KNN的决策边界
    knn_clf = KNeighborsClassifier()
    knn_clf.fit(X_train, y_train)
    print(knn_clf.score(X_test, y_test))
    plot_decision_boundary(knn_clf, axis=[4, 7.5, 1.5, 4.5])
    plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
    plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    plt.show()

运行结果

1.0

通过结果,我们可以看到KNN算法的预测数据集的分类准确率也是100%,但是KNN算法的决策边界就不再是一条直线,而是一条弯弯曲曲的曲线。在这条曲线的上方,用K-近邻的思路得到的点就是一个红色的点,而在曲线的下方得到的点用K-近邻的思路就是一个蓝色的点。我们知道KNN算法支持多分类,我们来看一下KNN多分类的决策边界

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from playLA.model_selection import train_test_split
from playLA.LogisticRegression import LogisticRegression
from matplotlib.colors import ListedColormap
from sklearn.neighbors import KNeighborsClassifier

if __name__ == "__main__":

    iris = datasets.load_iris()
    X = iris.data
    y = iris.target
    # 鸢尾花数据集有3个分类,而逻辑回归只能进行二分类,所以我们去掉2这个特征
    X = X[y < 2, :2]
    y = y[y < 2]
    # print(X.shape)
    # print(y.shape)
    # plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
    # plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    # plt.show()

    X_train, X_test, y_train, y_test = train_test_split(X, y, seed=666)
    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(y_test)
    # # 查看预测结果
    # print(log_reg.predict(X_test))
    # # 查看逻辑回归的系数
    # print(log_reg.coef)
    # # 查看逻辑回归的截距
    # print(log_reg.interception)

    # def x2(x1):
    #     return (- log_reg.coef[0] * x1 - log_reg.interception) / log_reg.coef[1]
    #
    # x1_plot = np.linspace(4, 8, 1000)
    # x2_plot = x2(x1_plot)
    # plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
    # plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    # plt.plot(x1_plot, x2_plot)
    # plt.show()

    # plt.scatter(X_test[y_test == 0, 0], X_test[y_test == 0, 1], color='red')
    # plt.scatter(X_test[y_test == 1, 0], X_test[y_test == 1, 1], color='blue')
    # plt.plot(x1_plot, x2_plot)
    # plt.show()

    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)
        custom_cmap = ListedColormap(['#EF9A9A', '#FFF59D', '#90CAF9'])
        plt.contourf(x0, x1, zz, linewidth=5, cmap=custom_cmap)

    # plot_decision_boundary(log_reg, axis=[4, 7.5, 1.5, 4.5])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
    # plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    # plt.show()

    # KNN的决策边界
    knn_clf = KNeighborsClassifier()
    knn_clf.fit(X_train, y_train)
    # print(knn_clf.score(X_test, y_test))
    # plot_decision_boundary(knn_clf, axis=[4, 7.5, 1.5, 4.5])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
    # plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    # plt.show()

    # KNN多分类的决策边界
    knn_clf_all = KNeighborsClassifier()
    knn_clf_all.fit(iris.data[:, :2], iris.target)
    plot_decision_boundary(knn_clf_all, axis=[4, 8, 1.5, 4.5])
    plt.scatter(iris.data[iris.target == 0, 0], iris.data[iris.target == 0, 1], color='red')
    plt.scatter(iris.data[iris.target == 1, 0], iris.data[iris.target == 1, 1], color='blue')
    plt.scatter(iris.data[iris.target == 2, 0], iris.data[iris.target == 2, 1], color='green')
    plt.show()

运行结果

通过图中KNN的决策边界,可以看出是非常不规则的,黄蓝相间的地方图形非常的奇怪,甚至在黄色中还存在一些蓝色的部分,这样一种弯弯曲曲的形状显然就是过拟合的表现。这是因为KNN默认的超参数k太小导致的,我们需要调大这个参数。

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from playLA.model_selection import train_test_split
from playLA.LogisticRegression import LogisticRegression
from matplotlib.colors import ListedColormap
from sklearn.neighbors import KNeighborsClassifier

if __name__ == "__main__":

    iris = datasets.load_iris()
    X = iris.data
    y = iris.target
    # 鸢尾花数据集有3个分类,而逻辑回归只能进行二分类,所以我们去掉2这个特征
    X = X[y < 2, :2]
    y = y[y < 2]
    # print(X.shape)
    # print(y.shape)
    # plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
    # plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    # plt.show()

    X_train, X_test, y_train, y_test = train_test_split(X, y, seed=666)
    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(y_test)
    # # 查看预测结果
    # print(log_reg.predict(X_test))
    # # 查看逻辑回归的系数
    # print(log_reg.coef)
    # # 查看逻辑回归的截距
    # print(log_reg.interception)

    # def x2(x1):
    #     return (- log_reg.coef[0] * x1 - log_reg.interception) / log_reg.coef[1]
    #
    # x1_plot = np.linspace(4, 8, 1000)
    # x2_plot = x2(x1_plot)
    # plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
    # plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    # plt.plot(x1_plot, x2_plot)
    # plt.show()

    # plt.scatter(X_test[y_test == 0, 0], X_test[y_test == 0, 1], color='red')
    # plt.scatter(X_test[y_test == 1, 0], X_test[y_test == 1, 1], color='blue')
    # plt.plot(x1_plot, x2_plot)
    # plt.show()

    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)
        custom_cmap = ListedColormap(['#EF9A9A', '#FFF59D', '#90CAF9'])
        plt.contourf(x0, x1, zz, linewidth=5, cmap=custom_cmap)

    # plot_decision_boundary(log_reg, axis=[4, 7.5, 1.5, 4.5])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
    # plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    # plt.show()

    # KNN的决策边界
    knn_clf = KNeighborsClassifier()
    knn_clf.fit(X_train, y_train)
    # print(knn_clf.score(X_test, y_test))
    # plot_decision_boundary(knn_clf, axis=[4, 7.5, 1.5, 4.5])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
    # plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    # plt.show()

    # KNN多分类的决策边界
    knn_clf_all = KNeighborsClassifier(n_neighbors=50)
    knn_clf_all.fit(iris.data[:, :2], iris.target)
    plot_decision_boundary(knn_clf_all, axis=[4, 8, 1.5, 4.5])
    plt.scatter(iris.data[iris.target == 0, 0], iris.data[iris.target == 0, 1], color='red')
    plt.scatter(iris.data[iris.target == 1, 0], iris.data[iris.target == 1, 1], color='blue')
    plt.scatter(iris.data[iris.target == 2, 0], iris.data[iris.target == 2, 1], color='green')
    plt.show()

运行结果

这个样子的决策边界比之前的样子要规整了很多。整体分成了三大块。

在逻辑回归中使用多项式特征

虽然逻辑回归只能进行二分类,但是直线这种分类方式还是太简单了。有很多情况并不是一根直线就能划分的。

对于上图中的样本点,我们是不可能用一根直线将它们分成两部分的。这个分布是一个非线性的分布。事实上,对于这张图,我们可以用一个圆形边界来将它分割成两部分。

对于一个圆形来说,它的方程是(r为半径),对于这个样本来说,它的决策边界就应该是这个样子。虽然我们之前讲的逻辑回归是不可能产生这样的一个决策边界的,但是类比于多项式回归,我们只需要引入多项式项就可以解决这个问题。我们在边界中,将整体看作一个特征,整体看作一个特征。前面的系数为θ1=1,前面的系数为θ2=1,而=。此时我们得到的决策边界针对来说还是一个线性的决策边界,但是针对x1和x2来说,它就是一个非线性的圆形了。我们可以从线性回归转换到多项式回归这个思路为我们的逻辑回归这个算法添加多项式项。基于这样的方式,我们可以对非线性的数据进行比较好的分类。我们得到的决策边界也可以是一个曲线的形状。我们可以添加x1^3、x2^3或者x1、x2等等得到一个非常复杂的边界。对于我们添加的多项式的值,那个degree可以任意的大,我们几乎可以使用这种方式构造出任意形状的决策边界。现在我们就用代码来模拟这样的多项式分类,先构建这样的样本数据集。

import numpy as np
import matplotlib.pyplot as plt

if __name__ == "__main__":

    np.random.seed(666)
    X = np.random.normal(0, 1, size=(200, 2))
    y = np.array(X[:, 0]**2 + X[:, 1]**2 < 1.5, dtype='int')
    plt.scatter(X[y == 0, 0], X[y == 0, 1])
    plt.scatter(X[y == 1, 0], X[y == 1, 1])
    plt.show()

运行结果

现在我们来使用逻辑回归,不添加多项式项看看是什么样的

import numpy as np
import matplotlib.pyplot as plt
from playLA.LogisticRegression import LogisticRegression

if __name__ == "__main__":

    np.random.seed(666)
    X = np.random.normal(0, 1, size=(200, 2))
    y = np.array(X[:, 0]**2 + X[:, 1]**2 < 1.5, dtype='int')
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()
    # 使用逻辑回归,不添加多项式项
    log_reg = LogisticRegression()
    log_reg.fit(X, y)
    print(log_reg.score(X, y))

运行结果

0.605

通过结果,我们可以看出,对所有数据进行分类(这里没有将样本数据集分成训练数据集和测试数据集),得到的分类准确度只有60.5%,显然这个准确度是比较低的。现在我们使用之前说的划分不规则边界的形式来看看它的决策边界。

import numpy as np
import matplotlib.pyplot as plt
from playLA.LogisticRegression import LogisticRegression
from matplotlib.colors import ListedColormap

if __name__ == "__main__":

    np.random.seed(666)
    X = np.random.normal(0, 1, size=(200, 2))
    y = np.array(X[:, 0]**2 + X[:, 1]**2 < 1.5, dtype='int')
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()
    # 使用逻辑回归,不添加多项式项
    log_reg = LogisticRegression()
    log_reg.fit(X, y)
    print(log_reg.score(X, y))

    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)
        custom_cmap = ListedColormap(['#EF9A9A', '#FFF59D', '#90CAF9'])
        plt.contourf(x0, x1, zz, linewidth=5, cmap=custom_cmap)

    plot_decision_boundary(log_reg, axis=[-4, 4, -4, 4])
    plt.scatter(X[y == 0, 0], X[y == 0, 1])
    plt.scatter(X[y == 1, 0], X[y == 1, 1])
    plt.show()

运行结果

由于我们的逻辑回归是由线性方式来处理的,所以导致这个决策边界非常的差,分类准确度只有60.5%。现在我们就为我们的逻辑回归算法添加多项式项。

import numpy as np
import matplotlib.pyplot as plt
from playLA.LogisticRegression import LogisticRegression
from matplotlib.colors import ListedColormap
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures
from sklearn.preprocessing import StandardScaler

if __name__ == "__main__":

    np.random.seed(666)
    X = np.random.normal(0, 1, size=(200, 2))
    y = np.array(X[:, 0]**2 + X[:, 1]**2 < 1.5, dtype='int')
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()
    # 使用逻辑回归,不添加多项式项
    log_reg = LogisticRegression()
    log_reg.fit(X, y)
    print(log_reg.score(X, y))

    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)
        custom_cmap = ListedColormap(['#EF9A9A', '#FFF59D', '#90CAF9'])
        plt.contourf(x0, x1, zz, linewidth=5, cmap=custom_cmap)

    # plot_decision_boundary(log_reg, axis=[-4, 4, -4, 4])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()

    # 使用逻辑回归,添加多项式项
    def PolynomialLogisticRegression(degree):
        return Pipeline([
            ('poly', PolynomialFeatures(degree)),
            ('std_scaler', StandardScaler()),
            ('log_reg', LogisticRegression())
        ])

    poly_log_reg = PolynomialLogisticRegression(degree=2)
    poly_log_reg.fit(X, y)
    print(poly_log_reg.score(X, y))

    plot_decision_boundary(poly_log_reg, axis=[-4, 4, -4, 4])
    plt.scatter(X[y == 0, 0], X[y == 0, 1])
    plt.scatter(X[y == 1, 0], X[y == 1, 1])
    plt.show()

运行结果

0.605
0.95

这样我们就绘制出了一个圆形的决策边界,而这个添加了多项式项的逻辑回归的识别准确率达到了95%。当然我们在实际中的数据集肯定是千奇百怪的,不可能是这样的一个圆形形状,我们可以通过调整超参数degree来看一下

import numpy as np
import matplotlib.pyplot as plt
from playLA.LogisticRegression import LogisticRegression
from matplotlib.colors import ListedColormap
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures
from sklearn.preprocessing import StandardScaler

if __name__ == "__main__":

    np.random.seed(666)
    X = np.random.normal(0, 1, size=(200, 2))
    y = np.array(X[:, 0]**2 + X[:, 1]**2 < 1.5, dtype='int')
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()
    # 使用逻辑回归,不添加多项式项
    log_reg = LogisticRegression()
    log_reg.fit(X, y)
    print(log_reg.score(X, y))

    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)
        custom_cmap = ListedColormap(['#EF9A9A', '#FFF59D', '#90CAF9'])
        plt.contourf(x0, x1, zz, linewidth=5, cmap=custom_cmap)

    # plot_decision_boundary(log_reg, axis=[-4, 4, -4, 4])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()

    # 使用逻辑回归,添加多项式项
    def PolynomialLogisticRegression(degree):
        return Pipeline([
            ('poly', PolynomialFeatures(degree)),
            ('std_scaler', StandardScaler()),
            ('log_reg', LogisticRegression())
        ])

    poly_log_reg = PolynomialLogisticRegression(degree=2)
    poly_log_reg.fit(X, y)
    print(poly_log_reg.score(X, y))

    # plot_decision_boundary(poly_log_reg, axis=[-4, 4, -4, 4])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()

    poly_log_reg2 = PolynomialLogisticRegression(degree=20)
    poly_log_reg2.fit(X, y)
    print(poly_log_reg2.score(X, y))

    plot_decision_boundary(poly_log_reg2, axis=[-4, 4, -4, 4])
    plt.scatter(X[y == 0, 0], X[y == 0, 1])
    plt.scatter(X[y == 1, 0], X[y == 1, 1])
    plt.show()

运行结果

0.605
0.95
0.955

我们可以看到在这个边界外面,决策边界变的非常的奇怪,显然这个外围的边界并不是我们想要的,出现这样的边界是因为degree取的20太大了,导致这个边界的形状非常的不规则。此时显然发生了过拟合的现象。我们要通过模型的正则化来限制过拟合。

scikit-learn中的逻辑回归

逻辑回归中使用正则化

之前我们在讲多项式回归中,正则化有L2岭回归以及L1 LASSO回归    ,总结下来就为

现在我们来看另外一种正则化的方式

这里的C就是一个新的超参数,如果C越大,我们在优化损失函数的时候,越会将J(θ)相应的减小到最小。但是如果C非常小的话,此时L1或L2的正则项就更加重要,我们在优化损失函数的时候,就会相应的让L1或者L2正则项更加的小。这就是为了限制θ尽可能的小。我们在J(θ)前加的这个C从某种程度上讲可以理解成在正则项前面加的α的倒数。我们依然有这两项,只不过为了平衡这两项,我们在J(θ)前加了这个系数。在逻辑回归中,在对模型进行正则化的时候,更偏好使用这种方式。在scikit-learn中的正则化的实现,就是使用这种方式。这是因为在使用这种方式的时候,你就必须使用正则化,因为L1、L2前面的系数不能为0,至少为1,如果要想让L1、L2不重要的话,就得让C非常的大。逻辑回归本身比较复杂,所以相应的加上正则化是比较重要的。现在我们就来用代码实现这个过程,先获取数据。

import numpy as np
import matplotlib.pyplot as plt

if __name__ == "__main__":

    np.random.seed(666)
    X = np.random.normal(0, 1, size=(200, 2))
    # 生成抛物线的数据
    y = np.array(X[:, 0]**2 + X[:, 1] < 1.5, dtype='int')
    # 添加噪音
    for _ in range(20):
        y[np.random.randint(200)] = 1
    plt.scatter(X[y == 0, 0], X[y == 0, 1])
    plt.scatter(X[y == 1, 0], X[y == 1, 1])
    plt.show()

运行结果

现在我们使用scikit-learn中的逻辑回归来进行分类

import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from matplotlib.colors import ListedColormap

if __name__ == "__main__":

    np.random.seed(666)
    X = np.random.normal(0, 1, size=(200, 2))
    # 生成抛物线的数据
    y = np.array(X[:, 0]**2 + X[:, 1] < 1.5, dtype='int')
    # 添加噪音
    for _ in range(20):
        y[np.random.randint(200)] = 1
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()

    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=666)
    log_reg = LogisticRegression()
    log_reg.fit(X_train, y_train)
    print(log_reg.score(X_train, y_train))
    print(log_reg.score(X_test, y_test))

    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)
        custom_cmap = ListedColormap(['#EF9A9A', '#FFF59D', '#90CAF9'])
        plt.contourf(x0, x1, zz, linewidth=5, cmap=custom_cmap)

    plot_decision_boundary(log_reg, axis=[-4, 4, -4, 4])
    plt.scatter(X[y == 0, 0], X[y == 0, 1])
    plt.scatter(X[y == 1, 0], X[y == 1, 1])
    plt.show()

运行结果

0.7933333333333333
0.86

这个跟之前是一样的,由于逻辑回归本身是线性性质的,而我们的数据集是一个抛物线形式的,所以它的决策边界不准确。现在我们要加入多项式项

import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from matplotlib.colors import ListedColormap
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures
from sklearn.preprocessing import StandardScaler

if __name__ == "__main__":

    np.random.seed(666)
    X = np.random.normal(0, 1, size=(200, 2))
    # 生成抛物线的数据
    y = np.array(X[:, 0]**2 + X[:, 1] < 1.5, dtype='int')
    # 添加噪音
    for _ in range(20):
        y[np.random.randint(200)] = 1
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()

    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=666)
    log_reg = LogisticRegression()
    log_reg.fit(X_train, y_train)
    print(log_reg.score(X_train, y_train))
    print(log_reg.score(X_test, y_test))

    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)
        custom_cmap = ListedColormap(['#EF9A9A', '#FFF59D', '#90CAF9'])
        plt.contourf(x0, x1, zz, linewidth=5, cmap=custom_cmap)

    # plot_decision_boundary(log_reg, axis=[-4, 4, -4, 4])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()

    # 使用逻辑回归,添加多项式项
    def PolynomialLogisticRegression(degree):
        return Pipeline([
            ('poly', PolynomialFeatures(degree)),
            ('std_scaler', StandardScaler()),
            ('log_reg', LogisticRegression())
        ])

    poly_log_reg = PolynomialLogisticRegression(degree=2)
    poly_log_reg.fit(X_train, y_train)
    print(poly_log_reg.score(X_train, y_train))
    print(poly_log_reg.score(X_test, y_test))
    plot_decision_boundary(poly_log_reg, axis=[-4, 4, -4, 4])
    plt.scatter(X[y == 0, 0], X[y == 0, 1])
    plt.scatter(X[y == 1, 0], X[y == 1, 1])
    plt.show()

运行结果

0.7933333333333333
0.86
0.9066666666666666
0.94

现在我们提高多项式项的阶数degree

import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from matplotlib.colors import ListedColormap
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures
from sklearn.preprocessing import StandardScaler

if __name__ == "__main__":

    np.random.seed(666)
    X = np.random.normal(0, 1, size=(200, 2))
    # 生成抛物线的数据
    y = np.array(X[:, 0]**2 + X[:, 1] < 1.5, dtype='int')
    # 添加噪音
    for _ in range(20):
        y[np.random.randint(200)] = 1
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()

    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=666)
    log_reg = LogisticRegression()
    log_reg.fit(X_train, y_train)
    print(log_reg.score(X_train, y_train))
    print(log_reg.score(X_test, y_test))

    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)
        custom_cmap = ListedColormap(['#EF9A9A', '#FFF59D', '#90CAF9'])
        plt.contourf(x0, x1, zz, linewidth=5, cmap=custom_cmap)

    # plot_decision_boundary(log_reg, axis=[-4, 4, -4, 4])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()

    # 使用逻辑回归,添加多项式项
    def PolynomialLogisticRegression(degree):
        return Pipeline([
            ('poly', PolynomialFeatures(degree)),
            ('std_scaler', StandardScaler()),
            ('log_reg', LogisticRegression())
        ])

    poly_log_reg = PolynomialLogisticRegression(degree=2)
    poly_log_reg.fit(X_train, y_train)
    print(poly_log_reg.score(X_train, y_train))
    print(poly_log_reg.score(X_test, y_test))
    # plot_decision_boundary(poly_log_reg, axis=[-4, 4, -4, 4])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()

    poly_log_reg2 = PolynomialLogisticRegression(degree=20)
    poly_log_reg2.fit(X_train, y_train)
    print(poly_log_reg2.score(X_train, y_train))
    print(poly_log_reg2.score(X_test, y_test))
    plot_decision_boundary(poly_log_reg2, axis=[-4, 4, -4, 4])
    plt.scatter(X[y == 0, 0], X[y == 0, 1])
    plt.scatter(X[y == 1, 0], X[y == 1, 1])
    plt.show()

运行结果

0.7933333333333333
0.86
0.9066666666666666
0.94
0.94
0.92

通过图像,我们发现此时决策边界变的不规则,很有可能发生了过拟合。现在我们加入模型的正则化

import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from matplotlib.colors import ListedColormap
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures
from sklearn.preprocessing import StandardScaler

if __name__ == "__main__":

    np.random.seed(666)
    X = np.random.normal(0, 1, size=(200, 2))
    # 生成抛物线的数据
    y = np.array(X[:, 0]**2 + X[:, 1] < 1.5, dtype='int')
    # 添加噪音
    for _ in range(20):
        y[np.random.randint(200)] = 1
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()

    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=666)
    log_reg = LogisticRegression()
    log_reg.fit(X_train, y_train)
    print(log_reg.score(X_train, y_train))
    print(log_reg.score(X_test, y_test))

    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)
        custom_cmap = ListedColormap(['#EF9A9A', '#FFF59D', '#90CAF9'])
        plt.contourf(x0, x1, zz, linewidth=5, cmap=custom_cmap)

    # plot_decision_boundary(log_reg, axis=[-4, 4, -4, 4])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()

    # 使用逻辑回归,添加多项式项
    def PolynomialLogisticRegression(degree):
        return Pipeline([
            ('poly', PolynomialFeatures(degree)),
            ('std_scaler', StandardScaler()),
            ('log_reg', LogisticRegression())
        ])

    poly_log_reg = PolynomialLogisticRegression(degree=2)
    poly_log_reg.fit(X_train, y_train)
    print(poly_log_reg.score(X_train, y_train))
    print(poly_log_reg.score(X_test, y_test))
    # plot_decision_boundary(poly_log_reg, axis=[-4, 4, -4, 4])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()

    poly_log_reg2 = PolynomialLogisticRegression(degree=20)
    poly_log_reg2.fit(X_train, y_train)
    print(poly_log_reg2.score(X_train, y_train))
    print(poly_log_reg2.score(X_test, y_test))
    # plot_decision_boundary(poly_log_reg2, axis=[-4, 4, -4, 4])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()

    def PolynomialLogisticRegression(degree, C):
        return Pipeline([
            ('poly', PolynomialFeatures(degree)),
            ('std_scaler', StandardScaler()),
            ('log_reg', LogisticRegression(C=C))
        ])

    poly_log_reg3 = PolynomialLogisticRegression(degree=20, C=0.1)
    poly_log_reg3.fit(X_train, y_train)
    print(poly_log_reg3.score(X_train, y_train))
    print(poly_log_reg3.score(X_test, y_test))
    plot_decision_boundary(poly_log_reg3, axis=[-4, 4, -4, 4])
    plt.scatter(X[y == 0, 0], X[y == 0, 1])
    plt.scatter(X[y == 1, 0], X[y == 1, 1])
    plt.show()

运行结果

0.7933333333333333
0.86
0.9066666666666666
0.94
0.94
0.92
0.84
0.92

在这里,我们需要知道scikit-learn中的逻辑回归是自带了L2正则项的,只不过默认的C超参数为1,我们调低了这个C为0.1。现在我们来使用L1正则项

import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from matplotlib.colors import ListedColormap
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures
from sklearn.preprocessing import StandardScaler

if __name__ == "__main__":

    np.random.seed(666)
    X = np.random.normal(0, 1, size=(200, 2))
    # 生成抛物线的数据
    y = np.array(X[:, 0]**2 + X[:, 1] < 1.5, dtype='int')
    # 添加噪音
    for _ in range(20):
        y[np.random.randint(200)] = 1
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()

    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=666)
    log_reg = LogisticRegression()
    log_reg.fit(X_train, y_train)
    print(log_reg.score(X_train, y_train))
    print(log_reg.score(X_test, y_test))

    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)
        custom_cmap = ListedColormap(['#EF9A9A', '#FFF59D', '#90CAF9'])
        plt.contourf(x0, x1, zz, linewidth=5, cmap=custom_cmap)

    # plot_decision_boundary(log_reg, axis=[-4, 4, -4, 4])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()

    # 使用逻辑回归,添加多项式项
    def PolynomialLogisticRegression(degree):
        return Pipeline([
            ('poly', PolynomialFeatures(degree)),
            ('std_scaler', StandardScaler()),
            ('log_reg', LogisticRegression())
        ])

    poly_log_reg = PolynomialLogisticRegression(degree=2)
    poly_log_reg.fit(X_train, y_train)
    print(poly_log_reg.score(X_train, y_train))
    print(poly_log_reg.score(X_test, y_test))
    # plot_decision_boundary(poly_log_reg, axis=[-4, 4, -4, 4])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()

    poly_log_reg2 = PolynomialLogisticRegression(degree=20)
    poly_log_reg2.fit(X_train, y_train)
    print(poly_log_reg2.score(X_train, y_train))
    print(poly_log_reg2.score(X_test, y_test))
    # plot_decision_boundary(poly_log_reg2, axis=[-4, 4, -4, 4])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()

    def PolynomialLogisticRegression(degree, C):
        return Pipeline([
            ('poly', PolynomialFeatures(degree)),
            ('std_scaler', StandardScaler()),
            ('log_reg', LogisticRegression(C=C))
        ])

    poly_log_reg3 = PolynomialLogisticRegression(degree=20, C=0.1)
    poly_log_reg3.fit(X_train, y_train)
    print(poly_log_reg3.score(X_train, y_train))
    print(poly_log_reg3.score(X_test, y_test))
    # plot_decision_boundary(poly_log_reg3, axis=[-4, 4, -4, 4])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()

    def PolynomialLogisticRegression(degree, C, penalty='l2'):
        return Pipeline([
            ('poly', PolynomialFeatures(degree)),
            ('std_scaler', StandardScaler()),
            ('log_reg', LogisticRegression(C=C, penalty=penalty, solver='liblinear'))
        ])


    poly_log_reg4 = PolynomialLogisticRegression(degree=20, C=0.1, penalty='l1')
    poly_log_reg4.fit(X_train, y_train)
    print(poly_log_reg4.score(X_train, y_train))
    print(poly_log_reg4.score(X_test, y_test))
    plot_decision_boundary(poly_log_reg4, axis=[-4, 4, -4, 4])
    plt.scatter(X[y == 0, 0], X[y == 0, 1])
    plt.scatter(X[y == 1, 0], X[y == 1, 1])
    plt.show()

运行结果

0.7933333333333333
0.86
0.9066666666666666
0.94
0.94
0.92
0.84
0.92
0.8266666666666667
0.9

通过使用L1正则,我们发现现在的决策边界跟我们的数据集本身的抛物线形式非常接近了,这是因为多项式阶数degree为20,它会有非常多的多项式项,而L1正则可以很好的消除过多的多项式项,让这些多项式项前面的系数为0,使得我们整个决策边界不会弯弯曲曲。不会像L2正则一样,产生这么奇怪的形状。由于我们这个数据本身就是二阶的多项式项生成的数据,所以在degree=20的时候,我们再怎么使用这种正则化的方式都无法达到让degree=2这种程度。我们面对实际的数据的时候,是不知道对于这个数据来说,我们的degree是多少是最合适的。为此,对于degree是多少,C是多少,我们的正则项使用L1还是L2,它们都属于超参数,我们都需要使用网格搜索这样的方式来寻找到最适合进行我们的这个数据相应的预测对应的这些参数,来形成对应的逻辑回归算法

OvR与OvO

我们知道逻辑回归只能解决二分类问题,那么对于多分类问题该怎么解决呢?解决多分类问题,有两种方式:

  1. OvR
  2. OvO

这两种算法不是只针对于逻辑回归的改造方式,而是对近乎所有的二分类算法的改造方式,让其都可以作用在多分类这样的问题上。

OvR(One vs Rest)

这个英文翻译过来就是一对剩余的所有的意思。

在这个图中,我们假设有4个类别。显然对于这一类分类问题,我们显然不能直接使用逻辑回归来解决。相应的我们选取其中的某一个类别(比如说我们选的是红色的类别),对于剩下的类别,我们将它们融合在一起,称之为其他类别,就是这个Rest的意思。

这样就把这个四分类问题转变成了二分类问题。现在我们就可以使用逻辑回归算法来估计一下我们给定一个新的样本点,它是红色的类别相应的概率是多少,它是非红色的类别相应的概率是多少。这样一个过程,我们同样可以在其他颜色的类别上进行,如果我们有四个类别的话,我们就可以形成4个不同的分类问题。

每一次都选择其中的一个类别和剩余的类别,进行这种二分类的任务。对于每一个二分类的任务,我们就都能得到新来的样本点它在其中的一个类别上对应的概率是多少。我们进一步拓展,如果我们有n个类别的话,就进行n次分类,选择分类得分最高的那个类别,我们就判断新来的样本点是在这个类别上。对于逻辑回归来说,这个得分就是指判断那个概率是多少。这样的方法也就被称为OvR。只不过我们处理问题的复杂度提升了,如果我们处理一个二分类的问题用的时间是T,那么现在有n个类别的话,我们就需要有n*T这么多时间来完成这个多分类任务。

OvO(One vs One)

这里其实就是一对一的进行比较。

在OvO中,我们直接就挑出其中的两个类别(比如我们挑出了红、蓝两个类别)

然后进行二分类任务,看对于这个二分类任务来说,我们判断我们新来的这个样本点,它是属于红色的类别还是蓝色的类别。这个过程可以重复进行,如果我们有四个类别的话,那么我们就能形成六个这种两两的对。这个计算方法就是排列组合的=4*3/2=6

它形成了6个二分类问题,对于每一个二分类问题,我们都可以估计一下我们要预测的这个样本点,它处于某一个两个类别中对应的哪一个类别,然后对于这些分类结果,我们就可以投票,看对于这6个分类结果来说,判定它在哪个类别中数量最大,我们最终就判定这个样本在哪个类别中。拓展下来,如果是一个n个类别的话,就进行次分类,选择赢数最高的分类作为我们最终的分类结果。比如说我们是手写数字识别的话,我们一共有10个类别,使用OvO的方式,就要进行=45次二分类任务,最终才可以得到对于一个新的数据来说,它是哪一种类别。很显然这种OvO的方式需要消耗的时间比OvR这种方式要多,这是因为是一个n方级别的数字,而对于OvR来说,我们只需要再进行n次二分类任务就好了。不过实际上对于OvO这个算法来说,虽然它耗时更多,但其实它的分类结果是更加准确的,这是因为我们每次只用真实的两个类别进行比较,而没有混淆其他的类别信息,所以它更倾向于是得到我们新来的那个样本,它真实的属于哪一个类别。现在我们依然使用鸢尾花数据集来进行代码实现。

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from matplotlib.colors import ListedColormap

if __name__ == "__main__":

    iris = datasets.load_iris()
    X = iris.data[:, :2]
    y = iris.target
    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=666)
    log_reg = LogisticRegression()
    log_reg.fit(X_train, y_train)
    print(log_reg.score(X_test, y_test))

    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)
        custom_cmap = ListedColormap(['#EF9A9A', '#FFF59D', '#90CAF9'])
        plt.contourf(x0, x1, zz, linewidth=5, cmap=custom_cmap)

    plot_decision_boundary(log_reg, axis=[4, 8.5, 1.5, 4.5])
    plt.scatter(X[y == 0, 0], X[y == 0, 1])
    plt.scatter(X[y == 1, 0], X[y == 1, 1])
    plt.scatter(X[y == 2, 0], X[y == 2, 1])
    plt.show()

运行结果

0.7894736842105263

这里需要说明的是,在scikit-learn中的逻辑回归,它自带了默认OvR这种多分类的模式,所以上图中对鸢尾花的数据集就是使用了OvR的方式进行了三分类。现在我们使用OvO的方式来进行这个多分类任务。

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from matplotlib.colors import ListedColormap

if __name__ == "__main__":

    iris = datasets.load_iris()
    X = iris.data[:, :2]
    y = iris.target
    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=666)
    # 使用OvR的方式进行多分类
    log_reg = LogisticRegression()
    log_reg.fit(X_train, y_train)
    print(log_reg.score(X_test, y_test))

    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)
        custom_cmap = ListedColormap(['#EF9A9A', '#FFF59D', '#90CAF9'])
        plt.contourf(x0, x1, zz, linewidth=5, cmap=custom_cmap)

    # plot_decision_boundary(log_reg, axis=[4, 8.5, 1.5, 4.5])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.scatter(X[y == 2, 0], X[y == 2, 1])
    # plt.show()
    # 使用OvO的方式进行多分类
    log_reg2 = LogisticRegression(multi_class='multinomial', solver='newton-cg')
    log_reg2.fit(X_train, y_train)
    print(log_reg2.score(X_test, y_test))
    plot_decision_boundary(log_reg2, axis=[4, 8.5, 1.5, 4.5])
    plt.scatter(X[y == 0, 0], X[y == 0, 1])
    plt.scatter(X[y == 1, 0], X[y == 1, 1])
    plt.scatter(X[y == 2, 0], X[y == 2, 1])
    plt.show()

运行结果

0.7894736842105263
0.7894736842105263

由于之前为了可视化,我们只使用了鸢尾花数据集的两个特征(维度)X = iris.data[:, :2],实际上这个数据集有4个特征,现在我们使用所有的特征来看一下

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from matplotlib.colors import ListedColormap

if __name__ == "__main__":

    iris = datasets.load_iris()
    X = iris.data[:, :2]
    y = iris.target
    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=666)
    # 使用OvR的方式进行多分类
    log_reg = LogisticRegression()
    log_reg.fit(X_train, y_train)
    print(log_reg.score(X_test, y_test))

    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)
        custom_cmap = ListedColormap(['#EF9A9A', '#FFF59D', '#90CAF9'])
        plt.contourf(x0, x1, zz, linewidth=5, cmap=custom_cmap)

    # plot_decision_boundary(log_reg, axis=[4, 8.5, 1.5, 4.5])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.scatter(X[y == 2, 0], X[y == 2, 1])
    # plt.show()
    # 使用OvO的方式进行多分类
    log_reg2 = LogisticRegression(multi_class='multinomial', solver='newton-cg')
    log_reg2.fit(X_train, y_train)
    print(log_reg2.score(X_test, y_test))
    # plot_decision_boundary(log_reg2, axis=[4, 8.5, 1.5, 4.5])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.scatter(X[y == 2, 0], X[y == 2, 1])
    # plt.show()

    X = iris.data
    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=666)
    log_reg = LogisticRegression()
    log_reg.fit(X_train, y_train)
    print(log_reg.score(X_test, y_test))

    log_reg2 = LogisticRegression(multi_class='multinomial', solver='newton-cg')
    log_reg2.fit(X_train, y_train)
    print(log_reg2.score(X_test, y_test))

运行结果

0.7894736842105263
0.7894736842105263
1.0
1.0

通过结果,我们可以看出,当使用了所有的特征的时候,逻辑回归对测试数据集的分类准确度无论是采用OvR还是OvO,都达到了100%。

之前我们说了,对于任意的二分类机器学习算法都可以使用OvR和OvO这两个算法来进行多分类,实际上scikit-learn也封装了专门的OvR和OvO两个类供我们使用。

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from matplotlib.colors import ListedColormap
from sklearn.multiclass import OneVsRestClassifier
from sklearn.multiclass import OneVsOneClassifier

if __name__ == "__main__":

    iris = datasets.load_iris()
    X = iris.data[:, :2]
    y = iris.target
    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=666)
    # 使用OvR的方式进行多分类
    log_reg = LogisticRegression()
    log_reg.fit(X_train, y_train)
    print(log_reg.score(X_test, y_test))

    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)
        custom_cmap = ListedColormap(['#EF9A9A', '#FFF59D', '#90CAF9'])
        plt.contourf(x0, x1, zz, linewidth=5, cmap=custom_cmap)

    # plot_decision_boundary(log_reg, axis=[4, 8.5, 1.5, 4.5])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.scatter(X[y == 2, 0], X[y == 2, 1])
    # plt.show()
    # 使用OvO的方式进行多分类
    log_reg2 = LogisticRegression(multi_class='multinomial', solver='newton-cg')
    log_reg2.fit(X_train, y_train)
    print(log_reg2.score(X_test, y_test))
    # plot_decision_boundary(log_reg2, axis=[4, 8.5, 1.5, 4.5])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.scatter(X[y == 2, 0], X[y == 2, 1])
    # plt.show()

    X = iris.data
    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=666)
    log_reg = LogisticRegression()
    log_reg.fit(X_train, y_train)
    print(log_reg.score(X_test, y_test))

    log_reg2 = LogisticRegression(multi_class='multinomial', solver='newton-cg')
    log_reg2.fit(X_train, y_train)
    print(log_reg2.score(X_test, y_test))

    ovr = OneVsRestClassifier(log_reg)
    ovr.fit(X_train, y_train)
    print(ovr.score(X_test, y_test))

    ovo = OneVsOneClassifier(log_reg)
    ovo.fit(X_train, y_train)
    print(ovo.score(X_test, y_test))

运行结果

0.7894736842105263
0.7894736842105263
1.0
1.0
0.9736842105263158
1.0

通过结果,我们可以看出当使用外部的OvR多分类器的时候,它的识别准确率下降了,只有97.36%,可见scikit-learn对于逻辑回归自带的OvR多分类器进行了优化。

评价分类结果

准确度的陷阱和混淆矩阵

我们至今为止都一直在使用分类准确度这一个指标来评判分类算法。那分类准确度这个指标来评判分类算法是不是就足够好了呢?其实分类准确度评判分类算法是有一个很重要的问题的。其实对于分类算法的评价比对回归算法的评价要复杂的多,相应的指标也多很多。具体选取什么样的指标来评价分类算法的时候还是要根据实际的情况,根据实际的数据,根据应用场景来进行恰当的选择。

分类准确度的问题

比如说有一个机器学习系统是用来进行癌症预测的,输入体检信息,可以判断是否有癌症。这个系统的预测准确度达到了99.9%。那这个系统是好的系统还是一个不好的系统呢?很多人可能会觉得这个系统的识别准确率这么高,肯定是一个好的系统,实际不然,这个系统实际是不够好的。因为如果这个癌症的发病率只有0.1%,那意味着这个系统只要预测所有人都是健康的,即可达到99.9%的准确率。这样一来这个系统近乎什么都没做,但是使用准确度这样一个指标来看的话,达到了99.9%。如果这种癌症产生的概率是0.01%,那意味着如果系统预测所有的人都是健康的,此时预测准确率可以达到99.99%,这比99.9%还要高,如果我们真的进行了机器学习训练,发现它的预测准确率是99.9%的话,这意味着这个机器学习系统是失败的。因为它比我纯粹的预测每个人都是健康的得到的准确度还要更低。这就是使用分类准确度来衡量一个分类系统的问题所在。这样的问题发生在我们的数据是极度偏斜(Skewed Data)的情况,我们在之前使用的数据集并没有这种极度偏斜的情况,无论是鸢尾花数据集还是手写识别数据集,这些数据集中各种分类的数量是差不多多的。但是对于癌症患者的数据发生的概率是非常小的,这是一种典型的数据发生了偏斜的情况,这样的数据集就是Skewed Data。面对这样的数据,如果我们只是使用分类准确度的话,可以得到一个非常高的分类准确度,但是其实我们的算法是不够好的,甚至在有一些情况下可能是非常烂的模型也能够得到一个很高的分类准确度。正因为这种情况,我们需要引入其他的分类指标来判断分类算法到底是好是坏。

使用混淆矩阵(Confusion Matrix)做进一步的分析

对于二分类问题,混淆矩阵是一个2*2的矩阵。

这里的行代表真实值,列代表预测值,0代表Negative,也就是负的,或者是阴性的;1代表Positive,也就是正的,或者是阳性的。比如我们在医院进行疾病的监测的时候,会用这种阴性、阳性来表示这种检测结果。如果我们的算法预测为0并且是正确的,放在(0,0)的位置,用TN表示,T表示是true,N表示Negative;如果我们的算法预测为1但是是错误的,放在(0,1)的位置,用FP表示,F表示是false,P表示Positive。如果我们的算法预测为0但是是错误的,放在(1,0)的位置,用FN表示,F表示是false,N表示Negative;如果我们的算法预测为1并且是正确的,放在(1,1)的位置,用TP表示,T表示是true,P表示Positive。

现在假设有10000个人,有一个算法对这10000个人的信息进行预测,预测他们是否患有癌症,这个结果是怎样的。这里用1代表患病,用0代表不患病,那这个混淆矩阵有可能就是下面这个样子

这个混淆矩阵表示的就是这10000个人中有9978个人他们本身并没有患癌症,同时我们的算法预测他们也没有患癌症,这就是TN;同时有12个人本身没有患癌症,可是我们的算法却错误的预测他们患了癌症,这就是FP;有2个人确实患了癌症,我们的算法却预测他们没有患癌症,这就是FN;有8个人确实患了癌症,同时我们的算法也正确的识别出了他们患了癌症,这就是TP。

精准率和召回率

精准率

在我们之前说的这个预测癌症病人中,这个精准率就为8/(8+12)=40%。对于精准率来说,它的分母就是橙色的这一列的数据和,也就是预测值为1的这一列的数据和。而它的分子就是预测为1,实际也是1,预测对了的那个数量。通常在有偏的数据中,我们通常将分类1作为我们关注的那个对象。比如说在医疗中,我们将1作为真正的患病。如果说这个混淆矩阵就是指预测癌症的混淆矩阵的话,那么精准率指的是预测癌症预测的成功率。对于这个40%来讲就是说我们每做100次预测,那么有40次是对的。如果这是一个申请信用卡的风险混淆矩阵,那么我们就把类别1表示有风险的那类人群,这是因为有风险的这类人群是我们的关注对象。精准率就是预测我们关注的那类对象,它预测的有多准

召回率

在我们之前说的这个预测癌症病人中,这个召回率就为8/(8+2)=80%。对于召回率来说,分母就是真实值为1的那一行的数据和。分子就是真实值为1,而预测值也为1,也就是说我们关注的那个事件发生了也被我们正确预测的数量。召回率的分子跟精准率的分子是一样的,不同的在分母。召回率的意思就是说我们关注的那个事件真实的发生了的这些数据中,我们成功的预测了多少。在这个癌症预测的人群中,有10个癌症患者,在这10个癌症患者中,我们成功的预测了8个。也就是说每当有100个癌症患者的话,这个机器学习癌症预测系统能够成功的找到其中的80个癌症患者。

在这个图中,所有的点代表了所有的样本,左边的部分代表我们关注的那一部分的样本,即真实分类为1的那一部分;右边的部分代表真实分类为0的那一部分。这个圈代表我们预测为1的部分,当然,我们的预测会出错,所以这个圈的左半部分代表我们预测成功的这部分,这个圈的右半部分代表我们预测失败的部分。

那么精准率就代表我们预测正确的左边的半圆/我们预测的所有的内容(整个圆)。

召回率就是我们预测正确的左边的半圆/我们真实的所有的关注对象(真实分类为1的左半部分)

为什么精准率和召回率比分类准确度是更好的指标呢?对于这种极度有偏的数据,比如说有10000人,癌症的发病率只有0.1%,我们预测所有人都是健康的,那么分类准确度就是99.9%。那么现在的混淆矩阵就是下面的样子

那么精准率=0/(0+0)无意义。此时我们定义精准率就为0。召回率=0/(10+0)=0。此时我们知道,虽然系统的分类准确度能达到99.9%,但是它的精准率和召回率都是最低值,我们通过这两个指标判断出这个预测算法它是完全没有用的,这就是精准率和召回率的意义。在极其有偏的数据中,我们不看分类准确率,而看精准率和召回率才能更加好的评价我们的分类系统它的好坏。

实现混淆矩阵,精准率和召回率

import numpy as np
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression

if __name__ == "__main__":

    digits = datasets.load_digits()
    X = digits.data
    y = digits.target.copy()
    # 产生偏斜
    y[digits.target == 9] = 1
    y[digits.target != 9] = 0
    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=666)
    log_reg = LogisticRegression()
    log_reg.fit(X_train, y_train)
    print(log_reg.score(X_test, y_test))
    y_log_predict = log_reg.predict(X_test)

    def TN(y_true, y_predict):
        # 预测值为0,真值为0
        assert len(y_true) == len(y_predict)
        return np.sum((y_true == 0) & (y_predict == 0))

    print(TN(y_test, y_log_predict))

    def FP(y_true, y_predict):
        # 预测值为1,真值为0
        assert len(y_true) == len(y_predict)
        return np.sum((y_true == 0) & (y_predict == 1))

    print(FP(y_test, y_log_predict))

    def FN(y_true, y_predict):
        # 预测值为0,真值为1
        assert len(y_true) == len(y_predict)
        return np.sum((y_true == 1) & (y_predict == 0))

    print(FN(y_test, y_log_predict))

    def TP(y_true, y_predict):
        # 预测值为1,真值为1
        assert len(y_true) == len(y_predict)
        return np.sum((y_true == 1) & (y_predict == 1))

    print(TP(y_test, y_log_predict))

    def confusion_matrix(y_true, y_predict):
        # 混淆矩阵
        return np.array([
            [TN(y_true, y_predict), FP(y_true, y_predict)],
            [FN(y_true, y_predict), TP(y_true, y_predict)]
        ])

    print(confusion_matrix(y_test, y_log_predict))

    def precision_score(y_true, y_predict):
        # 精准率
        tp = TP(y_true, y_predict)
        fp = FP(y_true, y_predict)
        try:
            return tp / (tp + fp)
        except:
            return 0

    print(precision_score(y_test, y_log_predict))

    def recall_score(y_true, y_predict):
        # 召回率
        tp = TP(y_true, y_predict)
        fn = FN(y_true, y_predict)
        try:
            return tp / (tp + fn)
        except:
            return 0

    print(recall_score(y_test, y_log_predict))

运行结果

0.9755555555555555
403
2
9
36
[[403   2]
 [  9  36]]
0.9473684210526315
0.8

在这里我们可以看到这些手写数据集经过偏斜后,它的TN为403,FP为2,FN为9,TP为36,精准率为94.7%,召回率为80%。现在来看一下scikit-learn中的混淆矩阵,精准率和召回率

import numpy as np
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import confusion_matrix
from sklearn.metrics import precision_score
from sklearn.metrics import recall_score

if __name__ == "__main__":

    digits = datasets.load_digits()
    X = digits.data
    y = digits.target.copy()
    # 产生偏斜
    y[digits.target == 9] = 1
    y[digits.target != 9] = 0
    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=666)
    log_reg = LogisticRegression()
    log_reg.fit(X_train, y_train)
    print(log_reg.score(X_test, y_test))
    y_log_predict = log_reg.predict(X_test)

    # def TN(y_true, y_predict):
    #     # 预测值为0,真值为0
    #     assert len(y_true) == len(y_predict)
    #     return np.sum((y_true == 0) & (y_predict == 0))
    #
    # print(TN(y_test, y_log_predict))
    #
    # def FP(y_true, y_predict):
    #     # 预测值为1,真值为0
    #     assert len(y_true) == len(y_predict)
    #     return np.sum((y_true == 0) & (y_predict == 1))
    #
    # print(FP(y_test, y_log_predict))
    #
    # def FN(y_true, y_predict):
    #     # 预测值为0,真值为1
    #     assert len(y_true) == len(y_predict)
    #     return np.sum((y_true == 1) & (y_predict == 0))
    #
    # print(FN(y_test, y_log_predict))
    #
    # def TP(y_true, y_predict):
    #     # 预测值为1,真值为1
    #     assert len(y_true) == len(y_predict)
    #     return np.sum((y_true == 1) & (y_predict == 1))
    #
    # print(TP(y_test, y_log_predict))
    #
    # def confusion_matrix(y_true, y_predict):
    #     # 混淆矩阵
    #     return np.array([
    #         [TN(y_true, y_predict), FP(y_true, y_predict)],
    #         [FN(y_true, y_predict), TP(y_true, y_predict)]
    #     ])
    #
    # print(confusion_matrix(y_test, y_log_predict))
    #
    # def precision_score(y_true, y_predict):
    #     # 精准率
    #     tp = TP(y_true, y_predict)
    #     fp = FP(y_true, y_predict)
    #     try:
    #         return tp / (tp + fp)
    #     except:
    #         return 0
    #
    # print(precision_score(y_test, y_log_predict))
    #
    # def recall_score(y_true, y_predict):
    #     # 召回率
    #     tp = TP(y_true, y_predict)
    #     fn = FN(y_true, y_predict)
    #     try:
    #         return tp / (tp + fn)
    #     except:
    #         return 0
    #
    # print(recall_score(y_test, y_log_predict))
    print(confusion_matrix(y_test, y_log_predict))
    print(precision_score(y_test, y_log_predict))
    print(recall_score(y_test, y_log_predict))

运行结果

0.9755555555555555
[[403   2]
 [  9  36]]
0.9473684210526315
0.8

F1 Score

对于精准率和召回率在我们使用的时候该怎么解读它们呢?其实对于这个问题和大多数机器学习领域关于取舍的问题是一样的,我们应该视具体使用场景而定。对于有一些场景来说,我们更加注重精准率,比如说我们做了一个机器学习的算法,这个算法是进行股票的预测。我们有可能将股票预测的问题做成一个二分类的问题,我们预测的是未来这个股票是升还是降。比如说我们关注的是升这种情况,预测到股票升这种情况我们给它分类为1,预测到股票降这种情况我们给它分类为0。在这种情况下,我们更加的关注精准率,精准率的意思就是我们做的所有分类为1的预测中有多少是正确的,换句话说我们做的所有的股票会升的预测有多少是对的,我们希望这个比例越高越好。如果我们预测股票要升的话,那么我就要相应的购买这些股票。但是如果我预测错了的话,犯了FP的错误,在这种情况下我们就会亏钱。在股票预测这个应用中,我们可能对召回率不是特别的关注,比如说有很多上升的股票,或者一支股票它有很多上升周期,在这里我们落掉了一些上升周期,也就是说本来这支股票上升为1,可是我们却错误的给它判断为0,也就是我们漏掉了这些股票上升的周期,那么其实并没有关系,因为我们漏掉了它们,我们也不会把钱投进去,所以我们并不会有损失。但是我们一旦判断它上升,但是这个判断结果却错误的话,会为我们带来实实在在的损失。在这种情况下,精准率就比召回率更加的重要。

而对于一些另外的应用场景,召回率就可能比精准率更加重要,比如在医疗领域,病人的诊断。此时如果召回率低,本来有一个病人他已经得病了,但是我们却没有正确的把他预测出来,这样的结果让这个病人的病情继续恶化下去就很有可能造成更加严重的后果,所以召回率非常重要。我们期望将所有有病的患者都能够正确的预测出来。但是此时精准率低一些并没有关系,因为在这种情况下,精准率低意味着我们犯了FP的错误,也就是有可能有一些人群他们没有得病,但是我们却错误的预测他有病。在这种情况下,对于这些人群,让他们再去做进一步的检查,进行确诊就好了。在进一步更深一步的检查中就会发现其实他们没有病。在此情况下,我们犯了FP的错误只是让一些人群多做了一些检查而已,但是我们关注召回率,也就是不想漏掉真正患病的人群。在这样的应用场景下,召回率就比精准率重要。在不同的应用场景,我们会偏好不同的指标,有些情况下,精准率更重要,有些情况下,召回率更重要。

还有一些情况可能并不像这两个例子这么极端,我们希望获得这两个指标的一个平衡。我们希望同时关注精准率和召回率。在这种情况下,我们使用一种新的指标,这个指标叫做——F1 Score

F1 Socre就是要兼顾精准率和召回率这两个指标。

这个式子本质是描述的精准率precision和召回率recall的调和平均值

所谓调和平均值就是precision和recall它们取倒数的和的结果的1/2,就是F1的倒数,这样求出来的F1就是这二者的调和平均值。调和平均值一个非常重要的特点就是如果这二者极度的不平衡的话,比如说一个值特别低,另外一个值特别高的话,那么我们最终得到的F1 Score的值也将特别的低。只有这二者都非常高的话,我们最终得到的这个F1 Score才会非常高。这个调和平均值的式子经过变型就变成了F1 Score的表达式。现在我们就用代码来看一下这个F1 Score。

import numpy as np

if __name__ == "__main__":

    def f1_score(precision, recall):
        try:
            return 2 * precision * recall / (precision + recall)
        except:
            return 0

    precision = 0.5
    recall = 0.5
    print(f1_score(precision, recall))

    precision = 0.1
    recall = 0.9
    print(f1_score(precision, recall))

    precision = 0
    recall = 1
    print(f1_score(precision, recall))

运行结果

0.5
0.18000000000000002
0.0

通过结果我们可以看出,只有当两个值都相似的时候,F1 Score的值才能达到一定量的值。而当其中一个很小,而另一个很大的时候,F1 Score依然还是很小。当其中一个为0的时候,结果就为0。现在我们依然使用手写数据集做偏斜处理,来看一下F1 Score的值,在这里我们使用scikit-learn为我们提供的F1 Score

import numpy as np
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import confusion_matrix
from sklearn.metrics import precision_score
from sklearn.metrics import recall_score
from sklearn.metrics import f1_score

if __name__ == "__main__":

    # def f1_score(precision, recall):
    #     try:
    #         return 2 * precision * recall / (precision + recall)
    #     except:
    #         return 0

    # precision = 0.5
    # recall = 0.5
    # print(f1_score(precision, recall))
    #
    # precision = 0.1
    # recall = 0.9
    # print(f1_score(precision, recall))
    #
    # precision = 0
    # recall = 1
    # print(f1_score(precision, recall))

    digits = datasets.load_digits()
    X = digits.data
    y = digits.target.copy()
    # 产生偏斜
    y[digits.target == 9] = 1
    y[digits.target != 9] = 0
    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=666)
    log_reg = LogisticRegression()
    log_reg.fit(X_train, y_train)
    print(log_reg.score(X_test, y_test))
    y_log_predict = log_reg.predict(X_test)
    print(confusion_matrix(y_test, y_log_predict))
    print(precision_score(y_test, y_log_predict))
    print(recall_score(y_test, y_log_predict))
    print(f1_score(y_test, y_log_predict))

运行结果

0.9755555555555555
[[403   2]
 [  9  36]]
0.9473684210526315
0.8
0.8674698795180723

最后得到的F1 Score的结果为86.74%,它比分类准确度97.55%要低的多,这是因为我们这个数据是有偏的,对于F1 Score这个指标可以更好衡量我们这个算法它的水平是怎么样的。

精准率和召回率的平衡

通过F1 Score,我们知道了,只有当精准率和召回率都比较大的时候,F1 Score才会比较大。那我们当然是希望我们的系统中精准率和召回率都是比较大的。但对于这样一个目标其实是实现不了的,因为精准率和召回率这二者之间是互相矛盾的两个指标。如果想让精准率提高,召回率就会不可避免的下降;如果想让召回率提高,精准率就会不可避免的下降。我们要做到的是找到精准率和召回率这两个指标之间的一个平衡。

为了理解这两个指标是互相矛盾的,我们从之前的逻辑回归的决策边界来看一下这个问题。我们知道逻辑回归的决策边界为,它是一根直线。如果我把这个决策边界变为

这个threshold就是一个阈值,这个阈值就是一个决策边界。如果≥threshold,我们就给它分类为1;<threshold,我们就给它分类为0。这样相当于为我们的逻辑回归算法引入了一个新的超参数threshold。我们指定这个threshold可以平移这个决策边界对应的这根直线,从而影响我们的分类结果。

在上图中,这根横轴表示决策边界的位置,而此时这个位置就在的地方。上面的圆形和五角星代表着样本数据,圆形代表真实值为0的数据,五角星代表着真实值为1的数据。从图上可以看出,大于0这边有5个样本,我们都把它们预测成了1;小于0这边有7个样本,我们都把它们预测成了0。所以我们此时的精准率和召回率就如下图所示

因为我们总共预测了5个数据,有4个是对的,所以精准率就为80%;总共有6个真实值为1的样本,而我们只预测对了4个,所以召回率为67%。现在我们来挪动阈值的位置,改变决策边界。

我们将阈值向右移动了,决策边界变成了右边这根直线。此时我们会发现

它的精准率提高了,因为我们只预测了2个数据,而这两个数据都预测正确了,精准率变成了100%;但是召回率下降了,因为有6个真实值为1的数据,而我们只预测出了2个,所以召回率就变成了33%。现在我们来向左移动阈值

我们将阈值向左移动了,决策边界变成了左边这根直线。此时我们会发现

它的召回率提升了,而精准率下降了,因为我们预测了8个数据,而只有6个数据预测正确,所以精准率为75%;总共有6个真实值为1的数据,而我们都预测出来了,所以召回率为100%。通过这样的一个过程,我们看到了当阈值取0的时候,精准率和召回率是多少;取比0大的时候,精准率和召回率是多少;取比0小的时候,精准率和召回率是多少。随着阈值的变大,精准率再逐渐的提高,相应的召回率在逐渐的减小;所以精准率和召回率是互相牵制,互相矛盾的一对指标。由此也可以看到,如果要提高精准率,我们只对那些特别有把握的数据,预测概率达到90%的时候,才把它们分类为1,那很多真实值为1的样本被排除了在判断值为1的外面,那么召回率相应的就会降低了。反过来,召回率要升高的话,我们就要降低预测概率,比如在癌症预测系统中,一个人有10%患癌的可能,我们都判定他是一个癌症患者,去做进一步的确诊,进一步的检查,所以我们就会拉低阈值,召回率确实得到了提高,不过此时精准率也就会下降。这就是精准率和召回率两个指标之间的平衡。现在我们通过代码来看一下如何调整这个阈值来改变精准率和召回率

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import confusion_matrix
from sklearn.metrics import f1_score
from sklearn.metrics import precision_score
from sklearn.metrics import recall_score

if __name__ == "__main__":

    digits = datasets.load_digits()
    X = digits.data
    y = digits.target.copy()
    # 产生偏斜
    y[digits.target == 9] = 1
    y[digits.target != 9] = 0
    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=666)
    log_reg = LogisticRegression()
    log_reg.fit(X_train, y_train)
    y_predict = log_reg.predict(X_test)
    print(f1_score(y_test, y_predict))
    print(confusion_matrix(y_test, y_predict))
    print(precision_score(y_test, y_predict))
    print(recall_score(y_test, y_predict))
    print(log_reg.decision_function(X_test))

运行结果

0.8674698795180723
[[403   2]
 [  9  36]]
0.9473684210526315
0.8
[-21.39111629 -32.89968323 -16.44016163 -79.83561566 -48.02243137
 -24.17617018 -44.63362641 -24.22571118  -1.14596638 -19.01171386
 -65.83248217 -50.98042089 -30.89891353 -45.95104769 -37.35822175
 -29.53882118 -36.91822879 -82.84001664 -37.66571977  -9.86437625
  -9.27829403 -85.28059442 -16.75452611 -45.33628474  -5.02567423
 -48.32625942 -11.65431236 -37.33148679 -25.08336919 -13.5901564
 -16.58231541 -28.78799912 -34.37955118 -28.52012209  -8.13542264
  -4.63313032 -21.87928168 -21.88869519 -31.09892577 -23.3984867
 -26.91006829 -62.24346371 -37.67365812 -66.36651609 -20.11619791
 -16.65433589 -18.17972943 -21.54203464 -28.96458855 -19.60879607
   2.43887956   7.71912841 -34.849044   -42.71549783 -25.64715678
 -34.7680629   -7.59448578 -49.53016396 -51.52838815  19.67482305
 -10.13108736 -32.02192086 -11.49856852  -1.4699478  -48.70848257
 -43.83625893 -24.84311106 -19.61286116 -36.66586896  -3.53065043
  -4.46961169 -19.22715107 -20.35477842 -40.90104406 -11.85320764
 -32.75099324 -35.76726754 -28.56663926 -55.41823901 -18.8405687
   4.57583101 -16.45881461 -76.79761689 -58.23275352 -30.22585503
 -29.42590775 -33.40990489  -8.41240435 -47.92542597 -65.5025185
 -16.93172073 -22.17307034 -11.27872628 -18.68046457 -69.24012583
 -46.38513493 -39.45623248 -35.95182192 -17.7478779  -62.96997249
 -16.87368524 -55.16843836 -28.78237436 -68.50484033 -68.89439645
  -6.49868293 -25.51499672 -38.34875938 -27.45932835 -15.5478568
 -27.48248247 -20.36181876  12.07838868 -23.10599106 -35.98724084
 -29.8672301  -68.97216398 -27.27941662 -54.28218303 -24.62341497
 -11.83237722 -47.37099878  -2.76558372 -59.69792222 -31.00151261
  -9.01540621 -70.85805094 -56.99696801 -20.05707198 -21.5152789
 -68.29244741 -18.92014998 -38.58907364 -57.36570694  -0.91281913
 -22.54992466 -22.66467292 -29.0235113  -32.72820688 -20.45664617
 -11.34951143   4.65601946   6.26977861   1.49141768  -7.6322439
 -39.25685273  12.16623393 -74.55460816 -75.09624528 -49.97601234
 -11.65031862 -47.60536987 -75.43101252 -29.91347448 -63.95185681
  -7.2689773   -6.63843399 -18.20822617 -32.46673697 -17.9546501
 -43.29361712 -32.68994198 -34.3180716  -72.74551409 -15.2019442
  11.46964243 -56.42500557  -6.03279066 -48.37246419 -16.43765927
  -2.13270161 -11.86214884 -33.25295195 -51.38609337 -10.38441353
 -17.19114075  -5.24520507 -25.21953049 -15.69850493   3.54659252
 -45.00033016 -12.56806987 -25.36295033 -16.5912305  -22.13313929
 -82.49022097  -5.87113613 -20.29637087 -20.47640811 -26.83351282
 -25.97467095 -40.49162548 -37.9859323  -26.96640644 -23.75413381
 -20.13250209  -9.67658136 -19.68893877 -42.53474809 -44.17265609
 -15.66069524 -64.04688059 -24.55090181 -56.30937433 -13.01786914
 -29.65196663   3.88268612 -44.34077042  -7.89931516   1.1386725
  -2.83023993 -11.93127768   7.50903966  -7.17392227 -46.37055337
 -48.655327    -4.5917329  -19.0386492  -24.07764547 -48.76584476
 -15.0310363  -24.94392526 -16.69726519 -18.67621386 -15.69759961
 -16.87847809 -38.54539917 -31.09627629  -9.38955202 -71.46839701
 -22.79405742 -14.42232531 -23.07557033 -34.32759766  -0.88099854
 -32.76513112 -11.2433954  -18.68720525  -8.20466874 -45.46311901
 -22.30356666 -62.43491449 -46.77168427 -65.15925742 -33.23520093
 -23.46541576 -28.48344277 -64.79554203   1.45558119  -4.09259583
 -25.67552639 -22.3174507  -54.70729372 -16.34598362 -12.09697892
 -35.30133487  -5.75958739 -13.47308585 -72.33098862  -6.16740739
  -1.17230738 -35.54773249 -24.18715193 -68.33044284  14.76753318
 -63.08785126   9.91043525 -24.14305198 -32.46286967 -14.40715008
 -85.75262185 -12.79528441   9.00228282 -16.503391   -36.69845123
 -16.512096   -19.35133411 -32.60588814  -5.63915201   7.69497867
   9.40915283   5.87506061 -35.64726625 -12.99147964 -54.43160989
 -41.10194425   5.62194842 -79.52182102 -15.81768172 -19.23006978
 -10.86823812 -42.52899396 -19.83047039 -15.69814303 -17.9881368
 -18.03782243  -6.75389093 -20.79051566 -16.58471449 -70.44670178
  -9.20322073 -31.68326159 -19.70349907 -21.98390417 -24.76579065
 -16.374491   -13.37493066 -22.92694274  11.05274035 -15.40526762
 -32.93952208 -13.75616482 -50.36690474 -20.48555219 -56.27962062
 -28.68557026 -21.86287716 -30.40747782 -69.26566596 -59.34914687
  14.35788248   8.58428239 -25.6877792    2.74470278   4.93805248
 -19.68361559 -58.84898426 -10.01585482 -28.80478682 -27.20954018
   6.30088159 -80.50052526 -34.46318844 -50.31508122 -35.96203418
 -48.65003927 -17.95728207 -62.35151791  -3.08443829 -25.258016
 -64.10656011  -9.63906144 -21.72007978  19.92910802 -18.75421777
  -4.47831729 -13.15140197 -21.63965789 -43.10672748 -52.13053585
 -28.51712477 -14.5823744   -2.46609074  -6.12681966   3.71455999
 -14.9977847  -40.84905871 -26.67149074  14.11187575 -17.69770862
  15.2076739  -33.08394198   5.28155289 -14.27159115 -53.5969583
 -50.03282267 -30.66725832 -38.04409354 -23.29943335 -24.69438718
 -13.55230615 -22.59913745 -27.21119584 -19.64039532 -28.17022763
 -19.93440541 -29.78707338 -11.30148933 -17.24723071 -24.0270052
 -24.36164647  10.39497979 -17.24456276 -38.0307681  -16.09768011
 -37.60773045 -16.35556132 -69.1377758  -33.70524439 -43.63203695
 -26.54795828 -10.30511187 -66.3722633  -31.88204434 -45.55916781
 -14.58978562 -36.09888421 -14.97052471 -70.02824527 -11.3627415
 -40.87243347 -32.6886302  -19.74769031 -27.56948378 -15.73011542
 -31.59281444  -8.52806892 -21.36544185 -34.08460531 -11.67038657
 -36.4478906  -34.76021702 -22.23118301   4.78055769 -21.33805183
  -4.46928482 -20.85352719 -32.24890674 -41.15145747 -25.07561195
 -19.76267733 -47.88375454 -30.95384768 -45.58272734 -71.52627785
  -6.26013362 -32.5437332    2.29251349  11.95767159   7.11098246
 -31.39894635 -63.96650655 -23.79370668  -5.74111066 -32.41736889
 -24.74623259 -67.7157729  -32.81326364 -33.60404409 -31.55108505
 -51.97532241 -22.53962738  -7.74759499 -17.29216614 -25.76348152
 -32.39491728 -29.51688896 -66.44155839 -45.69256747 -16.05903712]

由于scikit-learn自带的逻辑回归并不能通过调整预测概率值来改变分类的结果,但是我们可以通过log_reg.decision_function(X_test)来看到的决策边界值,它是一个一维数组。现在我们来取出它的前十个数据来看看预测值是多少。

print(log_reg.decision_function(X_test)[: 10])
print(log_reg.predict(X_test)[: 10])

运行结果

[-21.39111629 -32.89968323 -16.44016163 -79.83561566 -48.02243137
 -24.17617018 -44.63362641 -24.22571118  -1.14596638 -19.01171386]
[0 0 0 0 0 0 0 0 0 0]

通过结果,我们可以看到前十个决策边界值都是负的小于0,所以它们都被预测为0。因为逻辑回归默认是以0作为基准作为分类的。

现在我们要调大这个决策边界的基准为5,该怎么办呢?此时我们可以将这些决策边界值给储存起来,然后根据这些值的大小来重新生成预测值

decision_scores = log_reg.decision_function(X_test)
print(np.min(decision_scores))
print(np.max(decision_scores))
y_predict_2 = np.array(decision_scores >= 5, dtype='int')
print(confusion_matrix(y_test, y_predict_2))
print(precision_score(y_test, y_predict_2)) 
print(recall_score(y_test, y_predict_2))

运行结果

-85.75262185435525
19.929108017559418
[[404   1]
 [ 21  24]]
0.96
0.5333333333333333

此时我们会发现,当基准为5的时候,精准率变成了96%,比基准为0的时候的94.73%提升了,但是召回率变成了53.33%,比基准为0的时候的80%下降了。现在我们来调低基准为-5

y_predict_3 = np.array(decision_scores >= -5, dtype='int')
print(confusion_matrix(y_test, y_predict_3))
print(precision_score(y_test, y_predict_3))
print(recall_score(y_test, y_predict_3))

运行结果

[[390  15]
 [  5  40]]
0.7272727272727273
0.8888888888888888

我们会发现基准下降后,精准率变成了72.72%,比基准为0的时候的94.73%下降了,召回率变成了88.88%,比基准为0的时候的80%上升了。

精准率-召回率曲线

现在我们使用可视化的方式来看看精准率和召回率的关系。

precisions = []
recalls = []
# 定义从最小到最大的概率值,步长为0.1的队列
thresholds = np.arange(np.min(decision_scores), np.max(decision_scores), 0.1)
for threshold in thresholds:
    # 以队列中的每一个值为决策边界,获取预测的分类值
    y_predict = np.array(decision_scores >= threshold, dtype='int')
    # 添加该次的精准率
    precisions.append(precision_score(y_test, y_predict))
    # 添加该次的召回率
    recalls.append(recall_score(y_test, y_predict))
# 绘制精确率曲线
plt.plot(thresholds, precisions)
# 绘制召回率曲线
plt.plot(thresholds, recalls)
plt.show()

运行结果

图中蓝色的曲线为精准率,它随着阈值的增大在逐渐的增大;黄色的曲线为召回率,它随着阈值的增大在逐渐的降低。如果我们希望精准率在90%以上,我们就可以看到纵轴0.9对应的蓝色的线的位置对应的阈值(横轴)是多少。现在我们将精准率作为横轴,召回率作为纵轴,来看这样一个Precision-Recall曲线。

plt.plot(precisions, recalls)
plt.show()

运行结果

这个趋势也非常的明显,这根曲线是逐渐下降的,换句话说,随着精准率precisions逐渐增大,召回率recalls逐渐减小。再一次印证了这两个指标,它们是互相制约,互相平衡的。而在图中的红圈,我们可以看到这个位置这根曲线突然下降最快的位置,而这个位置也就是精准率和召回率一个最好的平衡位置。在这之后召回率将极度的下降,在这之前召回率下降的幅度并不高。这就是精准率-召回率曲线。现在我们来看看在scikit-learn中的调用精准率-召回率曲线。

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import precision_recall_curve

if __name__ == "__main__":

    digits = datasets.load_digits()
    X = digits.data
    y = digits.target.copy()
    # 产生偏斜
    y[digits.target == 9] = 1
    y[digits.target != 9] = 0
    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=666)
    log_reg = LogisticRegression()
    log_reg.fit(X_train, y_train)
    y_predict = log_reg.predict(X_test)

    decision_scores = log_reg.decision_function(X_test)
    precisions, recalls, thresholds = precision_recall_curve(y_test, decision_scores)
    print(precisions.shape)
    print(recalls.shape)
    print(thresholds.shape)

运行结果

(151,)
(151,)
(150,)

这里我们可以看到阈值列表比精准率和召回率列表要少一个值,这是因为precision_recall_curve返回的最后一个精准率为1召回率为0的那个值没有对应的阈值。我们用这三个列表就可以绘制出之前的两个图。

plt.plot(thresholds, precisions[: -1])
plt.plot(thresholds, recalls[: -1])
plt.show()

运行结果

通过这个图我们可以看到它跟我们之前绘制的图不一样,这是因为scikit-learn封装的precision_recall_curve中,自动寻找了它认为最重要的那一部分数据,但依然可以看到精准率和召回率是互相制约的关系。

plt.plot(precisions, recalls)
plt.show()

运行结果

现在有这样一个问题

虽然P-R曲线都是这个样子的,但假如有两个不同的算法或者对于同一个算法我们用两组超参数训练出了两个不同的P-R曲线,那么外面的这一根曲线对应的那个模型是优于里面的这根曲线对应的模型的。因为在外面的这根曲线上,它的每一个点对应的精准率值和召回率值都比里面的这根曲线要大。所以整体如果我们有一个模型,它的P-R曲线更靠外的话,那么通常这个模型就更加好,以此我们可以看到P-R模型可以作为选择模型,选择算法,选择超参数的指标。一般我们会用这个曲线和x,y轴包的面积来看这个模型的好坏,因为面积越大,这个模型就越好。不过一般我们不用P-R曲线的面积来衡量模型的优劣,相应的我们用另外一根曲线和x,y轴包围的面积来衡量我们对应模型的优劣,这个曲线就是非常著名的ROC曲线

ROC曲线

ROC的全称是Receiver Operation Characteristic Curve,它是描述TPR和FPR之间的关系。

TPR

TPR其实就是之前说的召回率

TPR就是预测为1,并且预测对了的数量占真实为1的百分比是多少。

FPR

FPR就是用FP去除以真实值为0的所有的数字和。FPR是指预测为1,可惜我们预测错了,这个数量占真实值为0的百分比是多少。

我们依然来看不同的阈值,对应的TPR和FPR是多少。在上图中最右边的这条阈值线,此时我们预测对了2个样本数据,真实值为1的样本数据有6个,所以TPR=33%,我们预测错的样本数量一个也没有,真实值为0的样本数量有6个,所以FPR为0%。再来看中间这条阈值线,此时我们预测对了4个样本数据,真实值为1的样本数据有6个,所以TPR=67%,我们预测错的样本数量有1个,真实值为0的样本数量有6个,所以FPR为16%。最后看一下最左边但阈值线,我们预测对的样本数量有6个,真实值为1的样本数量也是6个,所以TPR=100%,而我们预测错的样本数量有2个,真实值为0的样本数量有6个,所以FPR为33%。在这里我们可以看到,随着阈值逐渐降低的过程中,FPR是在逐渐升高的,TPR也是在逐渐升高的。FPR和TPR是呈现着一个相一致的趋势。TPR越高,FPR也越高;TPR越低,FPR也越低。这和精准率和召回率之间的关系是正好相反的。这是因为要提升TPR,我们需要拉低阈值,同时我们犯FP的错误也会增高,所以FPR也会增高。ROC曲线就是刻画TPR和FPR之间的关系。我们先把之前用到的代码全部放入我们自己的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之间的MAE"""

    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)

def TN(y_true, y_predict):
    # 预测值为0,真值为0
    assert len(y_true) == len(y_predict)
    return np.sum((y_true == 0) & (y_predict == 0))

def FP(y_true, y_predict):
    # 预测值为1,真值为0
    assert len(y_true) == len(y_predict)
    return np.sum((y_true == 0) & (y_predict == 1))

def FN(y_true, y_predict):
    # 预测值为0,真值为1
    assert len(y_true) == len(y_predict)
    return np.sum((y_true == 1) & (y_predict == 0))

def TP(y_true, y_predict):
    # 预测值为1,真值为1
    assert len(y_true) == len(y_predict)
    return np.sum((y_true == 1) & (y_predict == 1))

def confusion_matrix(y_true, y_predict):
    # 混淆矩阵
    return np.array([
        [TN(y_true, y_predict), FP(y_true, y_predict)],
        [FN(y_true, y_predict), TP(y_true, y_predict)]
    ])

def precision_score(y_true, y_predict):
    # 精准率
    tp = TP(y_true, y_predict)
    fp = FP(y_true, y_predict)
    try:
        return tp / (tp + fp)
    except:
        return 0

def recall_score(y_true, y_predict):
    # 召回率
    tp = TP(y_true, y_predict)
    fn = FN(y_true, y_predict)
    try:
        return tp / (tp + fn)
    except:
        return 0

def f1_score(y_true, y_predict):
    precision = precision_score(y_true, y_predict)
    recall = recall_score(y_true, y_predict)
    try:
        return 2 * precision * recall / (precision + recall)
    except:
        return 0

def TPR(y_true, y_predict):
    tp = TP(y_true, y_predict)
    fn = FN(y_true, y_predict)
    try:
        return tp / (tp + fn)
    except:
        return 0

def FPR(y_true, y_predict):
    fp = FP(y_true, y_predict)
    tn = TN(y_true, y_predict)
    try:
        return fp / (fp + tn)
    except:
        return 0
import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from playLA.metrics import FPR, TPR

if __name__ == "__main__":

    digits = datasets.load_digits()
    X = digits.data
    y = digits.target.copy()
    # 产生偏斜
    y[digits.target == 9] = 1
    y[digits.target != 9] = 0
    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=666)
    log_reg = LogisticRegression()
    log_reg.fit(X_train, y_train)
    # 获取所有测试数据集的概率值
    decision_scores = log_reg.decision_function(X_test)
    fprs = []
    tprs = []
    # 定义从最小到最大的概率值,步长为0.1的队列
    thresholds = np.arange(np.min(decision_scores), np.max(decision_scores), 0.1)
    for threshold in thresholds:
        # 以队列中的每一个值为决策边界,获取预测的分类值
        y_predict = np.array(decision_scores >= threshold, dtype='int')
        fprs.append(FPR(y_test, y_predict))
        tprs.append(TPR(y_test, y_predict))
    # 绘制ROC曲线
    plt.plot(fprs, tprs)
    plt.show()

运行结果

根据结果,我们可以看出随着FPR值的逐渐增高,TPR的值也在逐渐增高。我们来看一下scikit-learn中的ROC曲线。

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from playLA.metrics import FPR, TPR
from sklearn.metrics import roc_curve

if __name__ == "__main__":

    digits = datasets.load_digits()
    X = digits.data
    y = digits.target.copy()
    # 产生偏斜
    y[digits.target == 9] = 1
    y[digits.target != 9] = 0
    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=666)
    log_reg = LogisticRegression()
    log_reg.fit(X_train, y_train)
    # 获取所有测试数据集的概率值
    decision_scores = log_reg.decision_function(X_test)
    fprs = []
    tprs = []
    # 定义从最小到最大的概率值,步长为0.1的队列
    thresholds = np.arange(np.min(decision_scores), np.max(decision_scores), 0.1)
    for threshold in thresholds:
        # 以队列中的每一个值为决策边界,获取预测的分类值
        y_predict = np.array(decision_scores >= threshold, dtype='int')
        fprs.append(FPR(y_test, y_predict))
        tprs.append(TPR(y_test, y_predict))
    # 绘制ROC曲线
    # plt.plot(fprs, tprs)
    # plt.show()

    fprs, tprs, thresholds = roc_curve(y_test, decision_scores)
    plt.plot(fprs, tprs)
    plt.show()

运行结果

对于ROC曲线来说,我们关注的是曲线下面的面积,这个面积越大,可以说我们分类的模型越好。在横轴坐标越小的时候,也就是FPR越低的时候,我们犯FP的错误越少的时候,相应的这个曲线的值越高,也就是TPR值越大,也就是TP这种正确结果越多的时候,那么这根曲线整体就会被抬的越高,它底下的面积相应也就会越大。在这种情况下,我们的分类算法就会更好。所以ROC曲线下的面积可以作为我们分类算法的一个指标。这个面积,scikit-learn中也为我们提供了方法。

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from playLA.metrics import FPR, TPR
from sklearn.metrics import roc_curve
from sklearn.metrics import roc_auc_score

if __name__ == "__main__":

    digits = datasets.load_digits()
    X = digits.data
    y = digits.target.copy()
    # 产生偏斜
    y[digits.target == 9] = 1
    y[digits.target != 9] = 0
    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=666)
    log_reg = LogisticRegression()
    log_reg.fit(X_train, y_train)
    # 获取所有测试数据集的概率值
    decision_scores = log_reg.decision_function(X_test)
    fprs = []
    tprs = []
    # 定义从最小到最大的概率值,步长为0.1的队列
    thresholds = np.arange(np.min(decision_scores), np.max(decision_scores), 0.1)
    for threshold in thresholds:
        # 以队列中的每一个值为决策边界,获取预测的分类值
        y_predict = np.array(decision_scores >= threshold, dtype='int')
        fprs.append(FPR(y_test, y_predict))
        tprs.append(TPR(y_test, y_predict))
    # 绘制ROC曲线
    # plt.plot(fprs, tprs)
    # plt.show()

    fprs, tprs, thresholds = roc_curve(y_test, decision_scores)
    # plt.plot(fprs, tprs)
    # plt.show()
    # 获取ROC曲线下的面积
    print(roc_auc_score(y_test, decision_scores))

运行结果

0.9823868312757201

由于ROC曲线的定义域和值域都是在[0,1]之间的,所以它的总面积就为1。我们计算出来的ROC曲线下的面积为0.98,这是一个非常高的值。通过这个值,我们可以看出ROC曲线下的面积对非常有偏的数据并不敏感。不像精准率和召回率这两个指标那样对有偏的数据分类结果那么的敏感。所以针对有偏的数据,我们来看看它的精准率和召回率是非常有必要的。那么ROC曲线下的面积应用场合在哪里呢?

ROC的应用场合主要就在比较两个模型的优劣。上图中的两根曲线就代表了两个模型或者是同一个算法,它们使用超参数不同使用的模型对应的结果。在这种情况下,我们应该选择ROC曲线下面的面积更大的那个模型被认为是更好的模型。

多分类问题中的混淆矩阵

之前我们说的都是指二分类的,那么多分类来说,我们也可以看看对应的数据

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import precision_score
from sklearn.metrics import recall_score

if __name__ == "__main__":

    digits = datasets.load_digits()
    X = digits.data
    y = digits.target
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.8, random_state=666)
    log_reg = LogisticRegression()
    # 使用OvR的方式解决多分类的问题
    log_reg.fit(X_train, y_train)
    # 获取多分类准确度
    print(log_reg.score(X_test, y_test))
    y_predict = log_reg.predict(X_test)
    # 获取多分类的精准率
    print(precision_score(y_test, y_predict, average='micro'))
    # 获取多分类的召回率
    print(recall_score(y_test, y_predict, average='micro'))

运行结果

0.9408901251738526
0.9408901251738526
0.9408901251738526

现在我们重点看一下多分类的混淆矩阵,而混淆矩阵天然支持多分类

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import precision_score
from sklearn.metrics import recall_score
from sklearn.metrics import confusion_matrix

if __name__ == "__main__":

    digits = datasets.load_digits()
    X = digits.data
    y = digits.target
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.8, random_state=666)
    log_reg = LogisticRegression()
    # 使用OvR的方式解决多分类的问题
    log_reg.fit(X_train, y_train)
    # 获取多分类准确度
    print(log_reg.score(X_test, y_test))
    y_predict = log_reg.predict(X_test)
    # 获取多分类的精准率
    print(precision_score(y_test, y_predict, average='micro'))
    # 获取多分类的召回率
    print(recall_score(y_test, y_predict, average='micro'))
    # 获取这个十分类的混淆矩阵
    print(confusion_matrix(y_test, y_predict))

运行结果

0.9408901251738526
0.9408901251738526
0.9408901251738526
[[148   0   1   0   0   0   0   0   0   0]
 [  0 125   2   0   0   0   0   3   2  11]
 [  0   1 134   0   0   0   0   0   1   0]
 [  0   0   1 138   0   5   0   1   4   0]
 [  2   4   0   0 138   0   1   3   0   2]
 [  1   2   1   0   0 146   1   0   0   1]
 [  0   2   0   0   0   1 132   0   1   0]
 [  0   0   0   0   0   0   0 135   0   1]
 [  0   8   2   1   3   3   0   1 120   2]
 [  0   1   0   6   0   1   0   1   1 137]]

这是一个10*10的矩阵,就是我们这个十分类问题的混淆矩阵。这个混淆矩阵跟我们之前说的二分类的混淆矩阵是一样的,在这个矩阵中第i行第j列的位置表示的就是行方向依然是真值,也就是真值为i,列方向表示的是预测值,也就是将它预测为j的样本数量有多少。我们可以看到这个矩阵数字最大的是在对角线的位置,也就是第i行第i列的位置,也就是这个数字本身真值为i,同时我们也正确的预测为i相应的数量。不过在这里,我们依然可以看到在一些位置,会犯一些错误。我们现在来直观的把犯错的地方给绘制出来。

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import precision_score
from sklearn.metrics import recall_score
from sklearn.metrics import confusion_matrix

if __name__ == "__main__":

    digits = datasets.load_digits()
    X = digits.data
    y = digits.target
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.8, random_state=666)
    log_reg = LogisticRegression()
    # 使用OvR的方式解决多分类的问题
    log_reg.fit(X_train, y_train)
    # 获取多分类准确度
    print(log_reg.score(X_test, y_test))
    y_predict = log_reg.predict(X_test)
    # 获取多分类的精准率
    print(precision_score(y_test, y_predict, average='micro'))
    # 获取多分类的召回率
    print(recall_score(y_test, y_predict, average='micro'))
    # 获取这个十分类的混淆矩阵
    print(confusion_matrix(y_test, y_predict))
    cfm = confusion_matrix(y_test, y_predict)
    plt.matshow(cfm, cmap=plt.cm.gray)
    plt.show()

运行结果

这样我们就将数字矩阵给转化成了一个图像,越亮的地方数字越大,越暗的地方数字越小。但是我们只关注正确的部分是没有意义的,我们真正想找到的是我们预测犯错误的地方在哪。所以我们需要对绘制的这个矩阵做相应的处理。

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import precision_score
from sklearn.metrics import recall_score
from sklearn.metrics import confusion_matrix

if __name__ == "__main__":

    digits = datasets.load_digits()
    X = digits.data
    y = digits.target
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.8, random_state=666)
    log_reg = LogisticRegression()
    # 使用OvR的方式解决多分类的问题
    log_reg.fit(X_train, y_train)
    # 获取多分类准确度
    print(log_reg.score(X_test, y_test))
    y_predict = log_reg.predict(X_test)
    # 获取多分类的精准率
    print(precision_score(y_test, y_predict, average='micro'))
    # 获取多分类的召回率
    print(recall_score(y_test, y_predict, average='micro'))
    # 获取这个十分类的混淆矩阵
    print(confusion_matrix(y_test, y_predict))
    cfm = confusion_matrix(y_test, y_predict)
    # plt.matshow(cfm, cmap=plt.cm.gray)
    # plt.show()
    # 获取每行的样本
    row_sums = np.sum(cfm, axis=1)
    err_matrix = cfm / row_sums
    np.fill_diagonal(err_matrix, 0)
    print(err_matrix)
    plt.matshow(err_matrix, cmap=plt.cm.gray)
    plt.show()

运行结果

0.9408901251738526
0.9408901251738526
0.9408901251738526
[[148   0   1   0   0   0   0   0   0   0]
 [  0 125   2   0   0   0   0   3   2  11]
 [  0   1 134   0   0   0   0   0   1   0]
 [  0   0   1 138   0   5   0   1   4   0]
 [  2   4   0   0 138   0   1   3   0   2]
 [  1   2   1   0   0 146   1   0   0   1]
 [  0   2   0   0   0   1 132   0   1   0]
 [  0   0   0   0   0   0   0 135   0   1]
 [  0   8   2   1   3   3   0   1 120   2]
 [  0   1   0   6   0   1   0   1   1 137]]
[[0.         0.         0.00735294 0.         0.         0.
  0.         0.         0.         0.        ]
 [0.         0.         0.01470588 0.         0.         0.
  0.         0.02205882 0.01428571 0.07482993]
 [0.         0.00699301 0.         0.         0.         0.
  0.         0.         0.00714286 0.        ]
 [0.         0.         0.00735294 0.         0.         0.03289474
  0.         0.00735294 0.02857143 0.        ]
 [0.01342282 0.02797203 0.         0.         0.         0.
  0.00735294 0.02205882 0.         0.01360544]
 [0.00671141 0.01398601 0.00735294 0.         0.         0.
  0.00735294 0.         0.         0.00680272]
 [0.         0.01398601 0.         0.         0.         0.00657895
  0.         0.         0.00714286 0.        ]
 [0.         0.         0.         0.         0.         0.
  0.         0.         0.         0.00680272]
 [0.         0.05594406 0.01470588 0.00671141 0.02       0.01973684
  0.         0.00735294 0.         0.01360544]
 [0.         0.00699301 0.         0.04026846 0.         0.00657895
  0.         0.00735294 0.00714286 0.        ]]

在这个图像中,越亮的部分就是我们犯错最多的地方。我们将很多1预测成了9,将很多8预测成了1。有了这样的提示,其实我们就可以进一步改进我们的算法。我们得到的结果又把它规约成了一个二分类问题,我们现在的这个多分类很容易混淆1和9,也很容易混淆1和8。我们可以相应的微调在1和9、1和8这两个二分类问题中相应的阈值,来提高我们整体多分类问题的准确度。我们也应该把样本为1的图片,样本为9的图片和样本为8的图片都拿出来看看,感性的理解一下为什么会有那么多的分类错误的问题。在我们具体的要解决机器学习要解决的任务的时候,有些时候我们并不能通过算法来很好的解决这个问题。我们要回到数据中来看看数据有没有什么问题。我们能不能更好的整理我们的数据,清理我们的数据,提取我们数据中的特征等等,这些都是在机器学习领域非常重要的事情。

支撑向量机SVM

什么是SVM

SVM是Support Vector Machine的缩写。支撑向量机既可以解决分类问题,也可以解决回归问题。

在之前说逻辑回归的时候,我们提到了决策边界。

如果在决策边界的一侧,我们就认为这个样本点属于某一类,在另一侧就属于另外一类。但是对于一些样本数据,它们的决策边界并不唯一。

比方说在这个图中,这两条线都可以作为决策边界。对于这种决策边界不唯一的问题,通常有一个术语叫做不适定的问题。对于逻辑回归来说,定义了一个概率函数函数。根据这个概率函数进行建模,形成了一个损失函数,我们最小化这个损失函数,从而求出一条决策边界。在这里这个损失函数是完全由训练数据集所决定的。支撑向量机跟逻辑回归有一些不同,对于算法的泛化能力,对于完全未知的数据,能不能很好的分类是机器学习分类算法真正的目的。

现在假设有一个蓝色陌生的点,它离决策边界非常的近,虽然是蓝色,但它离决策边界太近了,很明显它离红色的点比其他蓝色的点更加的近,那么它极有可能是一个红色的点,只不过被决策边界分到了蓝色这边。这就是这个决策边界泛化能力较差的原因。根据这种思想,什么样的决策边界比较好呢?

那么这个决策边界就应该是这样的一根直线。这根直线的特别就是离这根直线最近的那些点离这根直线的距离都尽可能的远。

换句话说,我期望这根直线离红色的点尽可能的远,又离蓝色的点尽可能的远。同时它又能很好的分别红色和蓝色的数据点。所以就得到了这样的一个决策边界。它相当于是我要找到一个决策边界,这个决策边界不仅要将现在的训练数据集的样本进行一个很好的划分,它同时还考虑到了未来,它期望我们找到的这个决策边界相应的泛化能力尽可能的好。SVM的这种思想,对未来的泛化能力尽可能的好这种考量,没有寄望在数据的预处理阶段,或者是我们找到这个模型后再对模型进行正则化的方式。而是对泛化能力的考量直接放在了算法内部。也就是我们找到的决策边界离我们的分类样本都要尽可能的远。我们直观的看,这样的决策边界相应的泛化能力就是好的。实际上它背后是有数学理论的,在数学上可以严格的证明出在一个不适定的问题中,使用SVM的思路找到的这样的一个决策边界,它的泛化能力相应的是好的。SVM本身也是统计学中一种非常重要的方法,它的背后有极强的统计理论的支撑。

在这两个类别中,离决策边界最近的那些点,到这个决策边界的距离尽可能的远。在上面的图像中,红色类别有两个点,蓝色类别中有一个点,这三个点到这个决策边界的距离是一样的。同时这三个点也是所有的点中到这个决策边界距离最近的那些点。这些点它们又定义出了两根直线。

这两根直线跟决策边界是平行的。这两根直线相当于定义了一个区域,在这两根直线之间将不再有任何的数据点,而SVM最终得到的决策边界就是这个区域中间的那根线。SVM尝试寻找一个最优的决策边界,距离两个类别的最近的样本最远。而最近的这些点就被称为支撑向量

这也是支撑向量机这个名称的由来。这些支撑向量是SVM算法中非常重要的元素。而支撑向量到决策边界的距离被称为d,无论从哪一边到决策边界的d都是相等的。而两边的支撑向量的距离被称为marginSVM要做的事情就是最大化margin。目前我们说的都是线性可分问题。也就是说要存在一根直线或者在高维空间中存在一个超平面,可以将这些点划分。在这种情况下,我们才定义出了这个margin。这样的算法通常又叫做Hard Margin SVM。但是在真实情况下,我们的很多数据有可能是线性不可分的,在这种情况下,SVM可以进一步的改进成Soft Margin SVM。SVM要最大化margin是有数学语言来表达的,我们要求出这个margin的数学表达式,找到这个未知数的一组取值,能够最大化这个margin。

SVM背后的最优化问题

SVM要最大化margin,而margin=2d,所以只要我们找到这个d,就可以得到margin。

在二维平面上,一条直线可以写成Ax+By+C=0。如果有一个点(x,y)到这根直线的距离,就为

如果将二维平面扩展到n维空间,那么这条直线就为,这是包含了截距的写法,如果我们把截取给单独拿出来就为,其中w就是系数向量,x是特征向量,b就是截距。那么n维空间中有一个点到这条直线的距离就为

其中

这里w1、w2...wn表示向量w在每一个维度的分量。从这个距离公式,我们可以看一下

我们把决策边界的直线就定义成,那么决策边界两边的数据点到决策边界的距离都是大于等于d的,那么我们可以推导出这样一组式子

对于这两类数据来讲,我管一类叫做1,另一类叫做-1(这里跟之前的逻辑回归不同,逻辑回归是以1和0来分类的)。对于1这一类数据点来说,表示任取一个真实值为1的数据点到决策边界这条直线的距离,代入距离公式。由于这些点代入直线方程>0,所以距离公式去掉绝对值,就有了这组式子的第一条。对于-1这一类数据点来说,表示任取一个真实值为-1的数据点到决策边界这条直线的距离,代入距离公式。由于这些点代入直线方程<0,所以距离公式去掉绝对值,就有了这组式子的第二条。将上面的这组式子同时除以d,就有

由于||w||是一个数,d也是一个数,我们将和b分别除以这两个数的乘积,将该组式子化为

由此我们知道了决策边界两边的代表支撑向量的两根直线的方程就为

我们也可以将代表决策边界的这根直线的方程也化成跟代表支撑向量的两根直线方程保持一致

为了方便书写,现在我们将重新命名为和b,只不过它们已经不是之前的和b,它们之间相差一个系数

现在我们将上图中右边的这组式子合并成一个式子

这是因为当=1的时候,该式子与等价,当=-1的时候,该式子与等价。我们的支撑向量机就变成了对于我们所有的数据点都要满足这个不等式关系。在这种情况下,我们的目标是最大化d。对于任意支撑向量x,我们需要最大化

由于x是支撑向量,代表支撑向量的直线要么为1,要么为-1。所以=1,所以我们最大化的目标变成了

这相当于最小化w的模

通常我们求解这个问题,是

而这么写的原因也是因为方便求导数的操作。因此我们整个支撑向量机就变成了这样一个最优化的问题

我们最优化的目标函数是,但是它有一个限定条件,对于我们所有的数据要满足不等式,通常在数学表达上,这个条件写为,s.t.表示在这样的条件下,我们来最小化。这就是支撑向量机背后的最优化的问题。这跟我们之前讲的最优化问题(无论是线性回归,逻辑回归)是不一样的,之前的最优化问题是一个全局最优化问题,是没有限定条件的。而对于支撑向量机,它是一个有条件的最优化问题。加不加这个条件,在最优化领域方法是大不相同的。如果对于一个全局最优化问题,它没有条件,那只需要为我们目标函数求导,让它等于0相应的极值点就应该是取最大值或者是最小值的位置,对于每个点我们再进一步判断一下它是最大值还是最小值。但当是有条件的最优化问题的时候,这个求解方法变的复杂了很多,需要使用拉普拉斯算子来进行求解。把约束条件融合到优化目标函数中,将不等式约束转化为等式约束,从而将问题转化为拉格朗日求极值的问题。建立拉格朗日乘数法辅助函数:(关于拉格朗日乘数法求条件极值可以参考高等数学整理(二) 中的条件极值)

优化目标变为

求对偶问题

对w,b求极小值,即对w和b求偏导,并令其值为0,则得

把w和b回代入L(w,b,a),则得

求导,即可求出的极大值,采用SMO方法。

Soft Margin SVM

上一小节我们了解了对Hard Margin SVM的最优化问题

但是在这里有一个问题,如果我们训练数据中其中的一类中的支撑向量离另一类的支撑向量非常的近,那么决策边界就有可能是这个样子的

虽然这条直线将图中的两种分类给分开了,但是我们对这个决策边界的泛化能力产生怀疑。因为很显然对于大多数蓝色的点都集中在一个区域,只有一个蓝色的点在靠近红色点的位置,而决策边界非常强的受到这一个蓝色点的影响,从而形成了这样一根直线。但很有可能这个蓝色的点它根本就是一个错误的点或即使它是一个正确的点,但它是一个极度特殊的点,并不能代表一般情况。在这种情况下很有可能这样一个决策边界(与之前的决策边界交叉的那条直线)

是更好的一个决策边界。这个决策边界虽然将其中的一个蓝色的点进行了错误的分类,但很有可能我们将这样一个决策边界放到真实环境中,生产环境中进行实际的预测的时候,它的预测能力比之前那个将所有的训练数据集都进行正确分类的决策边界还要好。也就是所谓的泛化能力更加的高。所以我们要思考一个机制,这个机制对于我们的SVM算法来说,得到的这个决策边界,要能够有一定的容错能力。在一些情况下,它可以考虑到把一些点进行错误的分类。最终要达到的结果还是希望泛化能力尽可能的高。更一般的情况下,如果这个蓝色的点在这个位置

此时我们的数据根本是线性不可分的,没有任何一条直线能够将我们现在这个数据分成两类,在这种情况下,Hard Margin SVM的算法就已经不是泛化能力强不强的问题,而是根本无法应用,无法得到结果的问题了。不管从那个角度分析,我们都必须做出一个拥有容错能力的SVM。这种SVM就叫做Soft Margin SVM

之前我们在Hard Margin SVM中的限制条件就是,它表示在margin区域内不能再有任务点,所有的点都必须在margin区域外。

现在我们对这个条件加以宽松,让这些点不一定在margin外,这个宽松量为,那么这个条件就变成了

意思就是说我们允许一部分点在支撑向量的直线和这根虚线之间,而这根虚线和支撑向量的直线的距离就是

而这根虚线的方程就是这个不是一个固定的值,而是对于每一个样本数据i,都有一个相应的,并且≥0。换句话说如果我们有m个数据点的话,相应的也有m个。也就是说对于每一个数据点,我们都求出它对应的相应的容错空间。当然光有这个条件是不够的,如果让=+∞,也就是这根虚线在坐标空间下面无限远的地方,那么很显然对于我们所有的数据点都将满足这样的条件,那么此时我们这个容错的范围就太大了。我们希望这个能做的事情就是希望它有一定的容错空间,但是又不希望这个容错空间太大。我们怎么表征这个容错空间不能太大呢?

就是求出所有的和,放入目标函数中,求这个目标函数的最小化。这样一来,我们就顾及了Hard Margin SVM对应的最优化的内容,也就是做到了SVM的思想所要做的事情,同时表示我们设定了一堆,使得我们的SVM算法可以容忍一定的分类错误,但是这个容忍的程度要尽量的小。这二者之间应该取得一个平衡。这就是Soft Margin SVM对应的那个最优化问题相应的数学表达式。按照这个式子,这两部分所占的比例是同样重要的,但是按照我们之前在正则化的内容中,我们通常会加入一个超参数来调节它们的比例

如果C=1,代表这两部分一样重要;如果0<C<1的话,我们整个优化的式子将以优化前半部分为主;如果C特别大的话,这个式子主要优化的目标就是后半部分。我们可以使用网格搜索的方法来解决我们自己所要解决的问题和我们自己采集到的数据集而言,最好的C。

这里其实也相当于是给Soft Margin SVM中加入了L1正则。换句话说,我们加入的这个项本身是一个正则化项,它在避免我们训练出来的模型向一个极端方向发展,把这个模型往回拉一拉。对于我们的线性回归算法也好,逻辑回归算法也好,我们加入正则化项,其实它的本质也依然是让我们的模型相应的针对我们的训练数据集有更高的容错能力,当我们应用了这样的容错能力之后,使得我们的模型本身对训练数据集中的那些极端的数据点不那么敏感,用这种方式期望它的泛化能力得到提升,面对真实的未知的数据的时候,它的预测能力能够有所提升。当然有L1正则就有L2正则。L2正则就为

对比与线性回归的L1正则,为什么SVM的L1正则不是对的绝对值求和呢?因为我们已经限定了≥0,所以我们不需要加入绝对值。我们在逻辑回归中的正则超参数C是放在损失函数前面的,在SVM中,这个C放的位置虽然变了,但是它的表意和逻辑回归中那个模型中C的表意是一样的。它们的表意其实都是C越大,相应的容错空间越小。在SVM这个式子中,如果C取+∞,那么就意味着我们逼迫着这些都必须等于0,那么此时Soft Margin SVM就变成了Hard Margin SVM,也就是说我们给它的容错空间更小。在逻辑回归中,C在损失函数的前面,C越大就逼着我们要尽量多的去顾及损失函数那部分,而不去顾及后面的正则化,所以也是容错空间越小,而C越小,意味着我们可以有更大的容错空间

scikit-learn中的SVM

使用SVM算法和KNN一样,要做数据标准化处理。SVM算法寻找的是margin中间的那根线。而我们衡量margin的方式是数据点之间的距离。所以它是涉及距离的。如果我们的数据点在不同的维度上的量纲不同的话,就会使得我们对距离的估计是有问题的。

比如有4个数据点,在横轴上是在0到1之间,而纵轴是在0到10000之间,那么我们产生的决策边界就是一条横线。这是因为中间的两个点是支撑向量,这两个点离决策边界的距离是最大的。虽然我们现在看上去觉得这个距离很短,但是由于在纵向的尺度上是从0到10000之间的,所以在纵向范围里一个很短的距离都代表一个很大的数,此时我们的margin是这个样子的。

如果我们的纵轴也是在0到1之间的话。那我们的决策边界就是这个样子的。

在这种情况下,这4个点都是这个决策边界的支撑向量,我们的margin就是这个样子的

对于SVM算法来说,如果我们的特征在不同的维度上,它们数据尺度不同的话,就会非常严重的影响我们的SVM算法得到的最终的结果,那个决策边界是什么样子的。为了避免出现这种情况,在具体使用SVM之前,应该对数据进行标准化处理。现在依然使用鸢尾花数据集,先获取鸢尾花数据集。

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets

if __name__ == "__main__":

    iris = datasets.load_iris()
    X = iris.data
    y = iris.target
    # 只进行二分类,并且为了可视化只取前两个特征
    X = X[y < 2, :2]
    y = y[y < 2]
    plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
    plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    plt.show()

运行结果

首先,我们对数据集进行标准化处理,然后使用Hard Margin SVM来进行分类

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.preprocessing import StandardScaler
from sklearn.svm import LinearSVC
from matplotlib.colors import ListedColormap

if __name__ == "__main__":

    iris = datasets.load_iris()
    X = iris.data
    y = iris.target
    # 只进行二分类,并且为了可视化只取前两个特征
    X = X[y < 2, :2]
    y = y[y < 2]
    # plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
    # plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    # plt.show()
    standardScaler = StandardScaler()
    standardScaler.fit(X)
    X_standard = standardScaler.transform(X)
    # 生成一个Hard Margin SVM对象,调大超参数C
    svc = LinearSVC(C=1e9)
    svc.fit(X_standard, y)

    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)
        custom_cmap = ListedColormap(['#EF9A9A', '#FFF59D', '#90CAF9'])
        plt.contourf(x0, x1, zz, linewidth=5, cmap=custom_cmap)

    plot_decision_boundary(svc, axis=[-3, 3, -3, 3])
    plt.scatter(X_standard[y == 0, 0], X_standard[y == 0, 1], color='red')
    plt.scatter(X_standard[y == 1, 0], X_standard[y == 1, 1], color='blue')
    plt.show()

运行结果

现在我们调低超参数C,使用Soft Margin SVM来看看结果

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.preprocessing import StandardScaler
from sklearn.svm import LinearSVC
from matplotlib.colors import ListedColormap

if __name__ == "__main__":

    iris = datasets.load_iris()
    X = iris.data
    y = iris.target
    # 只进行二分类,并且为了可视化只取前两个特征
    X = X[y < 2, :2]
    y = y[y < 2]
    # plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
    # plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    # plt.show()
    standardScaler = StandardScaler()
    standardScaler.fit(X)
    X_standard = standardScaler.transform(X)
    # 生成一个Hard Margin SVM对象,调大超参数C
    svc = LinearSVC(C=1e9)
    svc.fit(X_standard, y)

    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)
        custom_cmap = ListedColormap(['#EF9A9A', '#FFF59D', '#90CAF9'])
        plt.contourf(x0, x1, zz, linewidth=5, cmap=custom_cmap)

    # plot_decision_boundary(svc, axis=[-3, 3, -3, 3])
    # plt.scatter(X_standard[y == 0, 0], X_standard[y == 0, 1], color='red')
    # plt.scatter(X_standard[y == 1, 0], X_standard[y == 1, 1], color='blue')
    # plt.show()
    svc2 = LinearSVC(C=0.01)
    svc2.fit(X_standard, y)
    plot_decision_boundary(svc2, axis=[-3, 3, -3, 3])
    plt.scatter(X_standard[y == 0, 0], X_standard[y == 0, 1], color='red')
    plt.scatter(X_standard[y == 1, 0], X_standard[y == 1, 1], color='blue')
    plt.show()

运行结果

在这里,我们会发现有一个红色的点被错误分类了,这就是我们将C给放小之后的结果。C越大容错空间越小,C越小容错空间越大。现在我们将margin的上下两条直线也绘制出来。先绘制Hard Margin SVM

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.preprocessing import StandardScaler
from sklearn.svm import LinearSVC
from matplotlib.colors import ListedColormap

if __name__ == "__main__":

    iris = datasets.load_iris()
    X = iris.data
    y = iris.target
    # 只进行二分类,并且为了可视化只取前两个特征
    X = X[y < 2, :2]
    y = y[y < 2]
    # plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
    # plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    # plt.show()
    standardScaler = StandardScaler()
    standardScaler.fit(X)
    X_standard = standardScaler.transform(X)
    # 生成一个Hard Margin SVM对象,调大超参数C
    svc = LinearSVC(C=1e9)
    svc.fit(X_standard, y)

    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)
        custom_cmap = ListedColormap(['#EF9A9A', '#FFF59D', '#90CAF9'])
        plt.contourf(x0, x1, zz, linewidth=5, cmap=custom_cmap)

    # plot_decision_boundary(svc, axis=[-3, 3, -3, 3])
    # plt.scatter(X_standard[y == 0, 0], X_standard[y == 0, 1], color='red')
    # plt.scatter(X_standard[y == 1, 0], X_standard[y == 1, 1], color='blue')
    # plt.show()
    svc2 = LinearSVC(C=0.01)
    svc2.fit(X_standard, y)
    # plot_decision_boundary(svc2, axis=[-3, 3, -3, 3])
    # plt.scatter(X_standard[y == 0, 0], X_standard[y == 0, 1], color='red')
    # plt.scatter(X_standard[y == 1, 0], X_standard[y == 1, 1], color='blue')
    # plt.show()

    def plot_svc_decision_boundary(model, axis):
        # 绘制不规则决策边界以及margin上下的两根直线
        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)
        custom_cmap = ListedColormap(['#EF9A9A', '#FFF59D', '#90CAF9'])
        plt.contourf(x0, x1, zz, linewidth=5, cmap=custom_cmap)
        w = model.coef_[0]
        b = model.intercept_[0]
        plot_x = np.linspace(axis[0], axis[1], 200)
        up_y = - w[0] / w[1] * plot_x - b / w[1] + 1 / w[1]
        down_y = - w[0] / w[1] * plot_x - b / w[1] - 1 / w[1]
        up_index = (up_y >= axis[2]) & (up_y <= axis[3])
        down_index = (down_y >= axis[2]) & (down_y <= axis[3])
        plt.plot(plot_x[up_index], up_y[up_index], color='black')
        plt.plot(plot_x[down_index], down_y[down_index], color='black')

    # 打印系数
    print(svc.coef_)
    # 打印截距
    print(svc.intercept_)
    plot_svc_decision_boundary(svc, axis=[-3, 3, -3, 3])
    plt.scatter(X_standard[y == 0, 0], X_standard[y == 0, 1], color='red')
    plt.scatter(X_standard[y == 1, 0], X_standard[y == 1, 1], color='blue')
    plt.show()

运行结果

[[ 4.03240941 -2.5070088 ]]
[0.92734009]

通过结果,我们可以看到,对于SVM算法来说相应的系数值,这个系数值有2个,这是因为对于我们的这个数据集来说它的特征有2个,这两个特征每个特征都对应着一个系数。在这里我们可以看到,返回的系数是一个二维的数组,这是因为scikit-learn封装的SVM算法直接可以处理多分类问题。如果有多个类别的话,算法就会相应的有多条直线来分割我们的这个特征平面,每一根直线相应的都有它的系数,所以这里是一个二维的数组。而返回的截距是一维数组,也是由于上面的原因产生的。从图像来看,由于我们现在用的是一个Hard Margin SVM算法,所以在margin中间是没有任何点的。现在我们对Soft Margin SVM也这样绘制一下。

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.preprocessing import StandardScaler
from sklearn.svm import LinearSVC
from matplotlib.colors import ListedColormap

if __name__ == "__main__":

    iris = datasets.load_iris()
    X = iris.data
    y = iris.target
    # 只进行二分类,并且为了可视化只取前两个特征
    X = X[y < 2, :2]
    y = y[y < 2]
    # plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
    # plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    # plt.show()
    standardScaler = StandardScaler()
    standardScaler.fit(X)
    X_standard = standardScaler.transform(X)
    # 生成一个Hard Margin SVM对象,调大超参数C
    svc = LinearSVC(C=1e9)
    svc.fit(X_standard, y)

    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)
        custom_cmap = ListedColormap(['#EF9A9A', '#FFF59D', '#90CAF9'])
        plt.contourf(x0, x1, zz, linewidth=5, cmap=custom_cmap)

    # plot_decision_boundary(svc, axis=[-3, 3, -3, 3])
    # plt.scatter(X_standard[y == 0, 0], X_standard[y == 0, 1], color='red')
    # plt.scatter(X_standard[y == 1, 0], X_standard[y == 1, 1], color='blue')
    # plt.show()
    svc2 = LinearSVC(C=0.01)
    svc2.fit(X_standard, y)
    # plot_decision_boundary(svc2, axis=[-3, 3, -3, 3])
    # plt.scatter(X_standard[y == 0, 0], X_standard[y == 0, 1], color='red')
    # plt.scatter(X_standard[y == 1, 0], X_standard[y == 1, 1], color='blue')
    # plt.show()

    def plot_svc_decision_boundary(model, axis):
        # 绘制不规则决策边界以及margin上下的两根直线
        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)
        custom_cmap = ListedColormap(['#EF9A9A', '#FFF59D', '#90CAF9'])
        plt.contourf(x0, x1, zz, linewidth=5, cmap=custom_cmap)
        w = model.coef_[0]
        b = model.intercept_[0]
        plot_x = np.linspace(axis[0], axis[1], 200)
        up_y = - w[0] / w[1] * plot_x - b / w[1] + 1 / w[1]
        down_y = - w[0] / w[1] * plot_x - b / w[1] - 1 / w[1]
        up_index = (up_y >= axis[2]) & (up_y <= axis[3])
        down_index = (down_y >= axis[2]) & (down_y <= axis[3])
        plt.plot(plot_x[up_index], up_y[up_index], color='black')
        plt.plot(plot_x[down_index], down_y[down_index], color='black')

    # 打印系数
    # print(svc.coef_)
    # 打印截距
    # print(svc.intercept_)
    # plot_svc_decision_boundary(svc, axis=[-3, 3, -3, 3])
    # plt.scatter(X_standard[y == 0, 0], X_standard[y == 0, 1], color='red')
    # plt.scatter(X_standard[y == 1, 0], X_standard[y == 1, 1], color='blue')
    # plt.show()
    print(svc2.coef_)
    print(svc2.intercept_)
    plot_svc_decision_boundary(svc2, axis=[-3, 3, -3, 3])
    plt.scatter(X_standard[y == 0, 0], X_standard[y == 0, 1], color='red')
    plt.scatter(X_standard[y == 1, 0], X_standard[y == 1, 1], color='blue')
    plt.show()

运行结果

[[ 0.43789691 -0.41091867]]
[0.00592624]

由于C=0.01,所以它的容错空间是比较大的,在margin中存在了很多点。

SVM中使用多项式特征和核函数

现在我们使用SVM来处理非线性数据的问题,首先我们先生成数据

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets

if __name__ == "__main__":

    # 制造一个月亮形状的数据集
    X, y = datasets.make_moons()
    print(X.shape)
    print(y.shape)
    plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
    plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    plt.show()

运行结果

(100, 2)
(100,)

通过结果,我们可以看到我们产生的这个月亮形数据集有100个样本,每个样本有2个特征。由于这个形状太规整了,我们给这些数据集加入噪音。

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets

if __name__ == "__main__":

    # 制造一个月亮形状的数据集
    # X, y = datasets.make_moons()
    # print(X.shape)
    # print(y.shape)
    # plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
    # plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    # plt.show()
    X, y = datasets.make_moons(noise=0.15, random_state=666)
    plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
    plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    plt.show()

运行结果

现在我们使用多项式特征的SVM算法来进行分类

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.preprocessing import PolynomialFeatures
from sklearn.preprocessing import StandardScaler
from sklearn.svm import LinearSVC
from sklearn.pipeline import Pipeline
from matplotlib.colors import ListedColormap

if __name__ == "__main__":

    # 制造一个月亮形状的数据集
    # X, y = datasets.make_moons()
    # print(X.shape)
    # print(y.shape)
    # plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
    # plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    # plt.show()
    X, y = datasets.make_moons(noise=0.15, random_state=666)
    # plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
    # plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    # plt.show()
    # 使用多项式特征的SVM
    def PolynomialSVC(degree, C=1):
        return Pipeline([
            ("poly", PolynomialFeatures(degree=degree)),
            ("std_scaler", StandardScaler()),
            ("linearSVC", LinearSVC(C=C))
        ])

    poly_svc = PolynomialSVC(degree=3)
    poly_svc.fit(X, y)

    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)
        custom_cmap = ListedColormap(['#EF9A9A', '#FFF59D', '#90CAF9'])
        plt.contourf(x0, x1, zz, linewidth=5, cmap=custom_cmap)

    plot_decision_boundary(poly_svc, axis=[-1.5, 2.5, -1, 1.5])
    plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
    plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    plt.show()

运行结果

通过图像我们可以看到决策边界不再是一条直线了,而是一条曲线。对于scikit-learn封装的SVM算法其实可以不使用PolynomialFeatures的方式,先对我们的数据升维,再放入LinearSVC中。SVM算法有一个特殊的方式可以直接使用多项式特征,这个特殊的方式就称为多项式核函数

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.preprocessing import PolynomialFeatures
from sklearn.preprocessing import StandardScaler
from sklearn.svm import LinearSVC
from sklearn.pipeline import Pipeline
from matplotlib.colors import ListedColormap
from sklearn.svm import SVC

if __name__ == "__main__":

    # 制造一个月亮形状的数据集
    # X, y = datasets.make_moons()
    # print(X.shape)
    # print(y.shape)
    # plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
    # plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    # plt.show()
    X, y = datasets.make_moons(noise=0.15, random_state=666)
    # plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
    # plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    # plt.show()
    # 使用多项式特征的SVM
    def PolynomialSVC(degree, C=1):
        return Pipeline([
            ("poly", PolynomialFeatures(degree=degree)),
            ("std_scaler", StandardScaler()),
            ("linearSVC", LinearSVC(C=C))
        ])

    poly_svc = PolynomialSVC(degree=3)
    poly_svc.fit(X, y)

    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)
        custom_cmap = ListedColormap(['#EF9A9A', '#FFF59D', '#90CAF9'])
        plt.contourf(x0, x1, zz, linewidth=5, cmap=custom_cmap)

    # plot_decision_boundary(poly_svc, axis=[-1.5, 2.5, -1, 1.5])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
    # plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    # plt.show()

    def PolynomialKernelSVC(degree, C=1):
        return Pipeline([
            ("std_scaler", StandardScaler()),
            ("kernelSVC", SVC(kernel='poly', degree=degree, C=C))
        ])

    poly_kernel_svc = PolynomialKernelSVC(degree=3)
    poly_kernel_svc.fit(X, y)
    plot_decision_boundary(poly_kernel_svc, axis=[-1.5, 2.5, -1, 1.5])
    plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red')
    plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue')
    plt.show()

运行结果

这里生成的决策边界虽然跟之前的决策边界不一样,但整体它依然是一个非线性的决策边界。

什么是核函数

我们之前在说SVM的目标函数进行最优化的时候

使用了拉格朗日乘数法来求条件极值。在对w、b求偏导,代入辅助函数L后,得到这样一个式子

通过这个式子,我们可以看到,它对数据集内任意两个特征向量都进行了点乘。之前我们在多项式回归中,给线性回归的特征向量添加了多项式特征进行升维。而升维的规则是类似于这个样子的

这相当于给原始数据增加了很多新的数据。相当于对进行了转化。

对于SVM的最优化问题来说,也就是说在进行多项式升维后的点乘也相应变成了的点乘。而核函数是不将特征向量进行升维扩大数据量,而是设置一个函数

直接转化成。也就是对原来的两个样本直接进行运算,直接计算出 ,如果我们能设计出这样一个函数K的话,那么之前求最优化的式子就变成了

这个K函数取不同的函数,其实就相当于在我们原始的样本上进行了一个变化,之后把变成了,把变成了,再对这两个新的样本进行点乘。这个K函数就叫做核函数。对于一些比较复杂的变型,通常当我们使用核函数的话,一方面计算量会有所减少,另外一方面节省了存储的空间,因为我们在对原始数据变型的时候,通常都是把低维的样本变成高维的数据,我们存储这个高维的数据相应就会花很多存储的空间。如果使用核函数的技巧的话,完全不用管把原始数据变成什么样子,不需要存储这个变化后的结果。核函数本身不是SVM特有的一种技巧,只要我们的算法转化成了一个最优化的问题,在求解这个最优化的问题中,存在这样的式子或者类似这样的式子,我们都可以使用核函数技巧。

多项式核函数

这是一个二次多项式的核函数。这里只是为了展示方便。x、y是两个向量,它们做点乘,我们将其展开为各个维度相乘的和的形式。

那么这个的转化就是这个样子

而y'也是这个样子,x'•y'其实就是

这个过程,我们可以不把这个变换后的数据x'、y'表示出来,然后再计算那么复杂的一个式子,我们可以直接用原始数据来计算就好了。这个结果和我们先变换x'、y'点乘的结果是一致的。这就是核函数的优势,它大大的降低了计算的复杂度。所以对于之前的SVM求最优化的式子就变成了

这样我们就完成了对我们的样本给它添加二次方这样的特征。之后再使用SVM算法进行运算。

我们的多项式核函数不仅仅能制造二次方这样的多项式特征,我们可以计算任意的阶数degree这样的特征。

这个幂d就是程序中的degree的值。同时这个1也可以是c。这个c和d就是两个超参数,这就是多项式核函数。而线性核函数其实就是

高斯RBF核函数

对于核函数是表示K(x,y)表示x和y的点乘。

对于高斯核函数,它的定义是这样的

这里的γ也是一个超参数。这个函数为什么叫做高斯核函数呢?首先我们来看一下正态分布(高斯分布)的式子(有关正态分布的介绍请参考概率论整理(二) )

在这里我们可以看到高斯核函数和正态分布是拥有一样的形态的,在高斯核函数中的的系数γ相当于正态分布中1/(2)。其实γ与1/(2)是有着共同的关系的。在正态分布中是,而在高斯核函数中是两个向量的差的模的平方。高斯核函数也被称为RBF(Radial Basis Function Kernel)核。在中文翻译中也称为径向基函数

之前我们在说多项式核函数的本质就是将所有的数据点首先添加了多项式项,再将这些有了多项式项的数据特征进行点乘,就形成了多项式核函数。相应的高斯核函数本质也应该是将原来的数据点先映射成一种新的特征向量,然后是这种新特征向量点乘的结果。高斯核函数表达出来这种数据的映射是非常复杂的。高斯核函数的本质是将每一个样本点映射到一个无穷维的特征空间。这里我们不做类似于多项式核函数映射后的推导,但是经过这样的变型之后进行的点乘却是非常简单的,就是,这就是核函数的威力,它可以让我们不需要具体的计算出来对于每一个样本点x和y到底变成了一个怎样的新的样本点,直接关注映射后的点乘的结果即可。我们用简单的例子来模拟一下高斯核函数在做什么样的事情。

我们先来看一下多项式特征为什么可以处理非线性的数据问题。它的基本原理是依靠升维使得原本线性不可分的数据线性可分。假设原本的数据只有一个特征x的话,现在我们为每一个x同时添加上一个x方,这样就可能使得原本线性不可分的数据线性可分。比如说我原本的数据就是一维数据。

这里一维的样本点有两个类,被分成了红色和蓝色的样本点。现在很显然这些数据是线性不可分的。我们没法画一根直线就把红色和蓝色完全分开。但是我们添加多项式特征的话,相当于我们是在把数据升维。也就是让每一个数据点不但有一个横轴维度的值,还有一个纵轴维度的值。在多项式特征中,我们让纵轴维度的值就是x的平方。

这些数据点在x轴上的位置是没有变的,只不过在y轴上它相应的也有一个取值了。这个取值就是x的平方。原来的数据点在一维上是线性不可分的,可是此时我们很容易找到一根直线把现在的两个类别区分开。

这就是升维的意义。我们依靠对原始数据进行升维可以将原本线性不可分的数据线性可分。其实高斯核做的事情也是这样一个事情。高斯核函数本来是这个样子的,现在我们将y固定,不取样本点而取固定的点,取两个固定点分别为l1和l2,这个l1、l2被称为地标

 

高斯核函数的升维的过程是,原本每一个x值有2个地标的话,这样就将一维的x映射到了二维空间。我们用代码来模拟一下这个过程。我们先生成线性不可分的一维数据。

import numpy as np
import matplotlib.pyplot as plt

if __name__ == "__main__":

    # 生成一维样本数据
    x = np.arange(-4, 5, 1)
    print(x)
    # 生成线性不可分的分类
    y = np.array((x >= -2) & (x <= 2), dtype='int')
    print(y)
    plt.scatter(x[y == 0], [0] * len(x[y == 0]))
    plt.scatter(x[y == 1], [0] * len(x[y == 1]))
    plt.show()

运行结果

[-4 -3 -2 -1  0  1  2  3  4]
[0 0 1 1 1 1 1 0 0]

现在我们使用高斯核函数的思路将一维的数据映射到二维空间。

import numpy as np
import matplotlib.pyplot as plt

if __name__ == "__main__":

    # 生成一维样本数据
    x = np.arange(-4, 5, 1)
    # print(x)
    # 生成线性不可分的分类
    y = np.array((x >= -2) & (x <= 2), dtype='int')
    # print(y)
    # plt.scatter(x[y == 0], [0] * len(x[y == 0]))
    # plt.scatter(x[y == 1], [0] * len(x[y == 1]))
    # plt.show()

    def gaussian(x, l):
        # 高斯核函数
        gamma = 1
        return np.exp(- gamma * (x - l)**2)

    # 设定地标
    l1, l2 = -1, 1
    X_new = np.empty((len(x), 2))
    for i, data in enumerate(x):
        X_new[i, 0] = gaussian(data, l1)
        X_new[i, 1] = gaussian(data, l2)
    plt.scatter(X_new[y == 0, 0], X_new[y == 0, 1])
    plt.scatter(X_new[y == 1, 0], X_new[y == 1, 1])
    plt.show()

运行结果

对于这样的一个二维数据,显然它是线性可分的。我们可以画出这样一条直线,将这些数据分成两部分,一侧是橙色的点,一侧是蓝色的点。

我们使用的是这样的映射将一维数据变成了二维数据,只不过我们将特征向量y那部分固定成了地标l1和l2。但是实际高斯核函数是,y其实是每一个数据点,高斯核函数干的事情其实和这里干的事情一样,只不过那个地标点取的比这里取的要多的多,我们的样本数据里有多少样本,它就有多少个地标点。换句话说,每一个样本x都要和样本y进行核函数的计算,成为高维数据相对应的某一个维度的元素。它是将一个m*n的数据映射成了m*m的数据。也就是说原本只有n个维度的数据的话,那么经过高斯核函数之后就将有m个维度的数据。如果m非常大的话,我们映射成了一个非常非常高维空间的样本数据。之前我们说高斯核函数可以将样本数据映射成无穷维的空间,是因为我们样本数据可以有无穷多个。但是具体我们拿到数据的时候,m是有限个,就映射成了一个m*m的数据。当我们使用高斯核函数的时候,这个计算开销是非常大的。SVM使用高斯核函数运行的时间会比较长。但是对于一些场景,运用高斯核函数是比较合适的,比如初始的样本数据维度就非常高,但是样本数据量可能并不多,如果我们的数据有m<n这样的特点的话,使用SVM这种高斯核函数就比较划算。最典型的应用的领域就是NPL(自然语言处理),在NPL领域通常会构建一个非常高维的特征空间,但是很有可能样本数量是并不多的。

高斯RBF核函数中的γ

我们依然来看一下高斯核函数和正态分布之间的关系

我们知道在正态分布中,μ代表均值,也就是数学期望,σ代表标准差,它的平方就是方差。

在正态分布中,μ决定了中心轴的位置,σ是用来描述样本数据分布的情况,σ越小,那么正态分布的图案就会越窄越集中;σ越大,正态分布图案就会越平缓越胖。在高斯核函数中,1/(2)这一项变成了γ,所以γ的趋势和σ的趋势刚好是相反的。在高斯核函数中,γ越大,正态分布越窄;γ越小,正态分布越宽。

scikit-learn中的高斯RBF核

首先我们依然先创建加入噪音的月亮型的数据集

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets

if __name__ == "__main__":

    X, y = datasets.make_moons(noise=0.15, random_state=666)
    plt.scatter(X[y == 0, 0], X[y == 0, 1])
    plt.scatter(X[y == 1, 0], X[y == 1, 1])
    plt.show()

运行结果

下面我们来使用SVM的高斯核

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.svm import SVC
from matplotlib.colors import ListedColormap

if __name__ == "__main__":

    X, y = datasets.make_moons(noise=0.15, random_state=666)
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()

    def RBFKernelSVC(gamma=1):
        # 生成一个高斯核函数的SVM模型
        return Pipeline([
            ("std_scaler", StandardScaler()),
            ("svc", SVC(kernel='rbf', gamma=gamma))
        ])

    svc = RBFKernelSVC()
    svc.fit(X, y)

    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)
        custom_cmap = ListedColormap(['#EF9A9A', '#FFF59D', '#90CAF9'])
        plt.contourf(x0, x1, zz, linewidth=5, cmap=custom_cmap)

    plot_decision_boundary(svc, axis=[-1.5, 2.5, -1, 1.5])
    plt.scatter(X[y == 0, 0], X[y == 0, 1])
    plt.scatter(X[y == 1, 0], X[y == 1, 1])
    plt.show()

运行结果

这就是使用高斯核函数,当γ=1的时候的决策边界。这个图像似乎跟之前加入了多项式项的图像没什么区别,现在我们来看改变γ后的图像,我们先调大γ的值

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.svm import SVC
from matplotlib.colors import ListedColormap

if __name__ == "__main__":

    X, y = datasets.make_moons(noise=0.15, random_state=666)
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()

    def RBFKernelSVC(gamma=1):
        # 生成一个高斯核函数的SVM模型
        return Pipeline([
            ("std_scaler", StandardScaler()),
            ("svc", SVC(kernel='rbf', gamma=gamma))
        ])

    svc = RBFKernelSVC()
    svc.fit(X, y)

    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)
        custom_cmap = ListedColormap(['#EF9A9A', '#FFF59D', '#90CAF9'])
        plt.contourf(x0, x1, zz, linewidth=5, cmap=custom_cmap)

    # plot_decision_boundary(svc, axis=[-1.5, 2.5, -1, 1.5])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()
    svc_gamma100 = RBFKernelSVC(gamma=100)
    svc_gamma100.fit(X, y)
    plot_decision_boundary(svc_gamma100, axis=[-1.5, 2.5, -1, 1.5])
    plt.scatter(X[y == 0, 0], X[y == 0, 1])
    plt.scatter(X[y == 1, 0], X[y == 1, 1])
    plt.show()

运行结果

之前我们说γ越大,正态分布越窄。上图中决策边界对于某一类的样本点,每一个样本点在它周围都形成了一个正态分布图案。我们可以想象成在俯视这个正态分布图案,每一个蓝色的点都是正态分布图中的一个尖。在每一个蓝色点的周围都围绕了一个决策区域,只在这样的区域中,才判定蓝色的分类。否则的话,都会被判定为橙色的分类。这也是高斯核函数直观的几何意义。我们这样得到的决策边界显然是过拟合了,它过于跟训练数据集相关了。现在我们稍微降低γ。

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.svm import SVC
from matplotlib.colors import ListedColormap

if __name__ == "__main__":

    X, y = datasets.make_moons(noise=0.15, random_state=666)
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()

    def RBFKernelSVC(gamma=1):
        # 生成一个高斯核函数的SVM模型
        return Pipeline([
            ("std_scaler", StandardScaler()),
            ("svc", SVC(kernel='rbf', gamma=gamma))
        ])

    svc = RBFKernelSVC()
    svc.fit(X, y)

    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)
        custom_cmap = ListedColormap(['#EF9A9A', '#FFF59D', '#90CAF9'])
        plt.contourf(x0, x1, zz, linewidth=5, cmap=custom_cmap)

    # plot_decision_boundary(svc, axis=[-1.5, 2.5, -1, 1.5])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()
    svc_gamma100 = RBFKernelSVC(gamma=100)
    svc_gamma100.fit(X, y)
    # plot_decision_boundary(svc_gamma100, axis=[-1.5, 2.5, -1, 1.5])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()
    svc_gamma10 = RBFKernelSVC(gamma=10)
    svc_gamma10.fit(X, y)
    plot_decision_boundary(svc_gamma10, axis=[-1.5, 2.5, -1, 1.5])
    plt.scatter(X[y == 0, 0], X[y == 0, 1])
    plt.scatter(X[y == 1, 0], X[y == 1, 1])
    plt.show()

运行结果

对于这个图,我们可以理解成对于每一个蓝色的点来说,它周围的正态分布变的宽了一些。所有蓝色点的正态分布区域融合在了一起就有了这样一个决策边界。现在我们来看一下γ<1的情况。

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.svm import SVC
from matplotlib.colors import ListedColormap

if __name__ == "__main__":

    X, y = datasets.make_moons(noise=0.15, random_state=666)
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()

    def RBFKernelSVC(gamma=1):
        # 生成一个高斯核函数的SVM模型
        return Pipeline([
            ("std_scaler", StandardScaler()),
            ("svc", SVC(kernel='rbf', gamma=gamma))
        ])

    svc = RBFKernelSVC()
    svc.fit(X, y)

    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)
        custom_cmap = ListedColormap(['#EF9A9A', '#FFF59D', '#90CAF9'])
        plt.contourf(x0, x1, zz, linewidth=5, cmap=custom_cmap)

    # plot_decision_boundary(svc, axis=[-1.5, 2.5, -1, 1.5])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()
    svc_gamma100 = RBFKernelSVC(gamma=100)
    svc_gamma100.fit(X, y)
    # plot_decision_boundary(svc_gamma100, axis=[-1.5, 2.5, -1, 1.5])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()
    svc_gamma10 = RBFKernelSVC(gamma=10)
    svc_gamma10.fit(X, y)
    # plot_decision_boundary(svc_gamma10, axis=[-1.5, 2.5, -1, 1.5])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()
    svc_gamma05 = RBFKernelSVC(gamma=0.5)
    svc_gamma05.fit(X, y)
    plot_decision_boundary(svc_gamma05, axis=[-1.5, 2.5, -1, 1.5])
    plt.scatter(X[y == 0, 0], X[y == 0, 1])
    plt.scatter(X[y == 1, 0], X[y == 1, 1])
    plt.show()

运行结果

现在我们继续调低γ值

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.svm import SVC
from matplotlib.colors import ListedColormap

if __name__ == "__main__":

    X, y = datasets.make_moons(noise=0.15, random_state=666)
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()

    def RBFKernelSVC(gamma=1):
        # 生成一个高斯核函数的SVM模型
        return Pipeline([
            ("std_scaler", StandardScaler()),
            ("svc", SVC(kernel='rbf', gamma=gamma))
        ])

    svc = RBFKernelSVC()
    svc.fit(X, y)

    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)
        custom_cmap = ListedColormap(['#EF9A9A', '#FFF59D', '#90CAF9'])
        plt.contourf(x0, x1, zz, linewidth=5, cmap=custom_cmap)

    # plot_decision_boundary(svc, axis=[-1.5, 2.5, -1, 1.5])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()
    svc_gamma100 = RBFKernelSVC(gamma=100)
    svc_gamma100.fit(X, y)
    # plot_decision_boundary(svc_gamma100, axis=[-1.5, 2.5, -1, 1.5])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()
    svc_gamma10 = RBFKernelSVC(gamma=10)
    svc_gamma10.fit(X, y)
    # plot_decision_boundary(svc_gamma10, axis=[-1.5, 2.5, -1, 1.5])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()
    svc_gamma05 = RBFKernelSVC(gamma=0.5)
    svc_gamma05.fit(X, y)
    # plot_decision_boundary(svc_gamma05, axis=[-1.5, 2.5, -1, 1.5])
    # plt.scatter(X[y == 0, 0], X[y == 0, 1])
    # plt.scatter(X[y == 1, 0], X[y == 1, 1])
    # plt.show()
    svc_gamma01 = RBFKernelSVC(gamma=0.1)
    svc_gamma01.fit(X, y)
    plot_decision_boundary(svc_gamma01, axis=[-1.5, 2.5, -1, 1.5])
    plt.scatter(X[y == 0, 0], X[y == 0, 1])
    plt.scatter(X[y == 1, 0], X[y == 1, 1])
    plt.show()

运行结果

现在我们可以看到,这个决策边界跟一个线性的差不多了。这显然是欠拟合了,不能非常好的反映出数据相应的样子。而γ值非常大的时候,又显然是过拟合了,对训练数据太敏感了,我们可以遇见它的泛化能力肯定是非常低。γ值相当于是在调整我们模型的复杂度,γ值越小,模型复杂度越低;γ值越高,模型复杂度就越高。γ值越低,模型更倾向于欠拟合,γ值越高,模型更倾向于过拟合。

SVM思想解决回归问题

回归问题的本质其实就是找到一根直线也好,曲线也好,能够最佳程度的拟合数据点。在这里怎么定义拟合其实就是不同回归算法的关键。对于线性回归来说,这个拟合就是定义数据点到这根直线的MSE均方误差值最小。而对于SVM这个算法来说,拟合的定义为在margin范围里,我们可以包含进的样本数据点越多越好。如果在margin范围里的数据点越多,就代表在margin范围里能够比较好的来表达我们的样本数据。在这种情况下,中间的直线作为真正的回归结果,用它来预测其他未知点相应的y的值。用SVM解决回归问题和解决分类问题是一个相反的思路,用SVM解决分类问题,我们期望的是在margin范围里的点越少越好。在极端的情况下Hard Margin SVM要求在margin区域里是一个点都不能有的。但是对于解决回归问题来说,恰恰相反,我们期望的是margin范围里,点越多越好。在我们具体训练SVM回归问题的结果的时候,我们是对margin的范围进行一个指定的,这里就引入了一个超参数ε,这个ε来指定margin有上下两根直线,任意的一根直线到中间的这根直线的距离。现在我们使用scikit-learn中的SVM来看一下回归问题。

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.svm import LinearSVR
# 可以传入不同的核函数
from sklearn.svm import SVR

if __name__ == "__main__":

    boston = datasets.load_boston()
    X = boston.data
    y = boston.target
    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=666)

    def StandardLinearSVR(epsilon=0.1):
        return Pipeline([
            ("std_scaler", StandardScaler()),
            ("linearSVR", LinearSVR(epsilon=epsilon))
        ])

    svr = StandardLinearSVR()
    svr.fit(X_train, y_train)
    print(svr.score(X_test, y_test))

运行结果

0.6361902835717916

当然这个R Squared的准确度是比较低的,我们可以通过调整其他超参数来提升。

展开阅读全文
打赏
0
0 收藏
分享
加载中
更多评论
打赏
0 评论
0 收藏
0
分享
返回顶部
顶部