逻辑回归

逻辑回归是一种解决分类问题的机器学习算法。

逻辑回归的思想是将样本特征和样本发生的概率联系起来,概率是一个0到1之间的数。线性回归中,通过回归方程可以得到一个预测值

而在逻辑回归中,通过回归函数得到的应该是一个概率值

p是一个概率,可以通过p去对应正事件或负事件,例如

可以表示当p>0.5时,得到的样本预测值对应正事件,反之对应负事件。

逻辑回归可以视为回归算法也可以视为分类算法,但通常用于分类,######但只适用于二分类问题!
对于线性回归



θ是θ0到θn的系数集合,Xb是在原有数据特征的基础上增加X0=1的一列之后得到的矩阵。通过f(X)得到的结果是一个在-∞ -> +∞之间的一个值,但是在逻辑回归中,我们需要的p是一个介于0和1之间的值,因此,我们需要对原有的回归函数做做手脚。

这个函数称为sigmoid函数。sigmoid函数表达式

t就是通过原有的回归函数得到的值,经过sigmoid函数的转化变成一个介于0到1之间的值。反映到代码中:
def sigmoid(t):
    return 1 / (1 + np.exp(-t))

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

绘制出sigmoid函数对应的曲线。



当sigmoid函数值t为0时,函数值为0.5,t越大,函数值越趋近于1,t越小越趋近于0。

逻辑回归的损失函数

通过之前的讲解,当得到的预测概率p大于0.5时,得到的预测值为1,否则为0。损失函数反应的是预测值和实际值之间的偏差情况,假设实际值为1,因为预测值根据预测概率得到,因此预测概率p越小,损失函数越大,因为p越小,低于0.5的概率就越高,以此得到的实际值为0的可能性月就越高,离实际的偏差也就越大。同样的,如果实际值为0,那么p越大,损失函数值就越大,道理同上。



给出逻辑回归的损失函数



下图是y=-log(x)的图像

因为p是一个0到1之间的值,因此横轴一下的部分可以不要,得到样本实际值y=1时的损失函数图像



对于y=0的情况,绘出相应的曲线

也是只取横轴以上的部分,加上之前y=1的曲线,得到的最终曲线

因为y只能去0或者1两个值,因此可以将损失函数归并为

对于样本集而言,含有m个样本的样本集,每个样本可以通过上面的cost函数求得单个样本的代价,求模型的代价函数则是通过样本集中每个样本的代价之和的均值得到。

将p代入代价函数J



这个函数没有公式解,只能通过梯度下降的方式求得最优解。

逻辑回归损失函数的梯度推导

根据梯度下降公式

需要求出代价函数对每个θ的偏导。首先从sigmoid函数入手。

求出sigmoid函数的导数后,返回看代价函数。把代价函数拆开来看。

先求出log(σ(t))和log(1-σ(t))的导数,接下来将整个代价函数拆分成两部分,依次求yi * log(σ(t))和(1-yi) * log(1-σ(t))对于θj的偏导。



求出log(σ(t))和log(1-σ(t))的导数后代入代价函数,就可以求出每个θj的导数公式。再将每个θj的导数代入梯度。

X0是添加的值为1的一列,因此对θ0的导数项可以省略X0,最终可以得到矩阵相乘的表达式结果。这也就是将在代码中实现的方式。
根据上述推导,模拟一次逻辑回归的批量梯度下降。
import numpy as np
from sklearn.metrics import mean_squared_error
import sklearn.datasets as dataset
from sklearn.model_selection import train_test_split

class Logistic_Regression:

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

    def _sigmoid(self, t):
        '''
        sigmoid函数
        :param t:
        :return:
        '''
        return 1 / (1 + np.exp(-t))

    def fit(self, X_train, y_train, eta=0.01, n_iters=1e4, epsilon=0.0001):
        '''
        梯度下降
        :param X_train:
        :param y_train:
        :param eta:
        :param n_iters:
        :return:
        '''
        assert X_train.shape[0] == y_train.shape[0], 'the size of X_train must equals the size of y_train'

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

        def dJ(theta, X_b, y):
            '''
            计算梯度
            :param theta:
            :param X_b:
            :param y:
            :return:
            '''
            return X_b.T.dot(self._sigmoid(X_b.dot(theta)) - y) / len(X_b)

        def gradient_descent(X_b, y, theta_init, eta=eta, n_iters=n_iters, epsilon=epsilon):
            '''
            梯度下降
            :param X_b:
            :param y:
            :param theta_init:
            :param eta:
            :param n_iters:
            :param epsilon:
            :return:
            '''
            theta = theta_init
            cur_iters = 0

            while cur_iters < n_iters:
                gradient = dJ(theta, X_b, y)
                last_theta = theta
                theta = theta - gradient * eta

                if abs(J(theta, X_b, y) - J(last_theta, X_b, y)) < epsilon:
                    break
                cur_iters = cur_iters + 1

            return theta

        X_b = np.hstack([np.ones((len(X_train), 1)), X_train])
        theta_init = np.zeros(X_b.shape[1])
        self._theta = gradient_descent(X_b, y_train, theta_init, eta, n_iters, epsilon)
        self.interept_ = self._theta[0]
        self.coef_ = self._theta[1:]
        return self

    def predict_probability(self, X_test):
        '''
        预测概率函数
        :param X_test:
        :return:
        '''
        assert self.coef_ is not None, 'coef can not be None'
        assert X_test.shape[1] == len(self.coef_), 'the size of X_test must equals the size of coef'

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


    def predict(self, X_test):
        '''
        预测函数
        :param X_test:
        :return:
        '''
        assert self.coef_ is not None, 'coef can not be None'
        assert X_test.shape[1] == len(self.coef_), 'the size of X_test must equals the size of coef'

        prob = self.predict_probability(X_test)
        return np.array(prob >= 0.5, dtype='int')

    def mse(self, X_test, y_test):
        '''
        测试预测准确度
        :param X_test:
        :param y_test:
        :return:
        '''
        y_predict = self.predict(X_test)
        return mean_squared_error(y_predict, y_test)

data = dataset.load_iris()
X = data.data
y = data.target
X = X[y < 2, :2]
y = y[y<2]
X_train, X_test, y_train, y_test = train_test_split(X, y)
logistics = Logistic_Regression()
logistics.fit(X_train, y_train)
mse = logistics.mse(X_test, y_test)
print(mse)

定义一个Logistic_Regression类,通过coef_和interept_记录特征的系数和回归曲线的截距。fit()函数封装了梯度下降相关的所有函数,包括代价函数J、求梯度的函数dJ以及梯度下降函数gradient_descent(),gradient_descent()函数的执行过程和总结线性回归时的批量梯度下降基本一样,不过是代价函数和求梯度的方式有所变化而已。之后可以调用predict()方法进行测试预测并使用mse函数求预测值与实际值之间的均方误差检验模型预测效果。得到输出

mse is  0.0
the coef of 2 features is  [ 1.19406112 -2.03850496]
the intercept is  -0.21296952322024035

可见预测准确度100%正确,相应的特征系数θ和截距θ0也能打印得到。

决策边界

逻辑回归本质上是一种二分类问题,决策边界就是将样本划分为不同类别的一条边界。



sigmoid函数将回归函数得到的值求得一个介于0和1之间的概率,当p>0.5时,预测值为1,p<0.5时,预测值为0.那么p=0.5就是一个临界值,此时e的系数就是0,也就是说

因此将上面的式子称为决策边界。
以二维特征矩阵为例,含有两个特征的决策边界为:

以此推得

于是我们定义求解X2的函数
def X2(logistic, X1):
    return (-logistic.intercept_ - logistic.coef_[0] * X1) / logistic.coef_[1]

logistic是一个逻辑回归的实例,intercept_为回归曲线的截距也就是θ0,coef_记录从θ1起的所有特征系数。

descision_boundary = X2(logistic, X_train[:][0])

还是以鸢尾花数据集为例,只取前两个特征,依据数据的特点,在区间4到8之间绘制决策边界。

data = dataset.load_iris()
X = data.data
y = data.target
X = X[y < 2, :2]
y = y[y < 2]

logistic = Logistic_Regression()
logistic.fit(X_train, y_train)

X1_plot = np.linspace(4, 8, 1000)
X2_plot = X2(logistic, 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()

得到决策边界
注意决策边界并非逻辑回归专有,很多其他机器学习模型也有决策边界的概念,如KNN等。

逻辑回归与多项式回归

之前模拟的逻辑回归都是基于一次项的逻辑回归,在实际环境中,一次项并不能很好的适应数据,因此,多项式回归和逻辑回归相结合,可以得到更好的逻辑回归模型。
构造虚拟数据

#模拟测试用例
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()

根据表达式,这个曲线应该是一个圆,将圆内和圆外的点分为不同类别。打印图像。



使用之前定义的逻辑回归模型进行训练

logistic = Logistic_Regression()
X_train, X_test, y_train, y_test = train_test_split(X, y)

logistic.fit(X_train, y_train)
mse = logistic.score(X_test, y_test)
print('theta[0] is %s, theta[1] is %s, interceptor is %s' % (logistic.coef_[0], logistic.coef_[1], logistic.intercept_))
print('regression is %s * X1 + %s * X2 + %s' % (logistic.coef_[0], logistic.coef_[1], logistic.intercept_))
print('the mse is ', mse)

得到输出结果

theta[0] is 0.0004486610126586672, theta[1] is 0.0001887681532517242, interceptor is 0.0007333333333333333
regression is 0.0004486610126586672 * X1 + 0.0001887681532517242 * X2 + 0.0007333333333333333
the mse is  0.42

可以看到截距和两个特征分别对应的系数,这种情况下的均方误差为0.42。
打印决策边界的图像

def X2(logistic, X1):
    return (-logistic.intercept_ - logistic.coef_[0] * X1) / logistic.coef_[1]

descision_boundary = X2(logistic, X_train[:][0])

X1_plot = np.linspace(-4, 4, 1000)
X2_plot = X2(logistic, 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()

得到图像



可以看出,这个决策边界并不能很好的划分不同类别。因此,需要使用多项式逻辑回归进行更细致的划分。借助先前讲过的sklearn中的PipeLine进行多项式逻辑回归。

def polynomialLogisticRegression(degree):
    '''
    多项式逻辑回归
    :param degree:
    :return:
    '''
    return Pipeline([
        ('poly', PolynomialFeatures(degree=degree)),
        ('std', StandardScaler()),
        ('Logistic_Regression', Logistic_Regression())
    ])

poly_log_reg = polynomialLogisticRegression(degree=2)
poly_log_reg.fit(X_train, y_train)
print('the mse of polynomial logistic regression is ', poly_log_reg.score(X_test, y_test))

得到输出结果

the mse of polynomial logistic regression is  0.1

明显的再使用多项式回归后,误差减小。

逻辑回归中使用正则化

在使用多项式逻辑回归的时候容易出现多项式项数过高引起过拟合,因此需要进行正则化。结合之前讲过的L1正则和L2正则,本节给出sklearn中将多项式逻辑回归和模型正则化一起使用的方法。

在讲解多项式回归时,将正则化项乘以正则化系数追加在代价函数后面。

逻辑回归的正则化表达式有所不同,将正则化系数乘以代价函数再加上正则化项。

C越大,与损失函数的乘积也就越大,这样在回归过程中若想达到达到极小值点,就需要将代价函数尽可能小。C越小,乘积越小,若想到达极小值点时模型约精确,就需要注意正则化项的值不能太小,否则求出来的模型准确度不高。

先模拟数据

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

#模拟200个样本,每个样本有2个特征
np.random.seed(666)
X = np.random.normal(0, 1, size=(2000, 2))
y = np.array(X[:, 0]**2 + X[:, 1] < 1.5, dtype='int')
#添加噪音,强制样本20个点值为1
for _ in range(20):
    y[np.random.randint(200)] = 1

X_train, X_test, y_train, y_test = train_test_split(X, y, 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()

模拟2000个样本,每个样本含有两个特征。y = X1^2 + 1.5 * X2,并随机模拟200个噪声数据。

logistic = LogisticRegression()
logistic.fit(X_train, y_train)
score = logistic.score(X_test, y_test)
print(logistic.coef_)
print(logistic.intercept_)
print('C = 1, penalty = l2 score is ', score)

logistic2 = LogisticRegression(C=0.1)
logistic2.fit(X_train, y_train)
score = logistic2.score(X_test, y_test)
print(logistic2.coef_)
print(logistic2.intercept_)
print('C = 0.1, penalty = l2 score is ', score)

logistic3 = LogisticRegression(C=0.1, penalty='l1')
logistic3.fit(X_train, y_train)
score = logistic3.score(X_test, y_test)
print(logistic3.coef_)
print(logistic3.intercept_)
print('C = 1, penalty = l1 score is ', score)

sklearn的逻辑回归中,默认的C=1且使用L2正则,通过修改C和penalty参数进行调参得到不同结果。
得到输出

[[ 0.12528696 -1.18653217]]
[1.06721009]
C = 1, penalty = l2 score is  0.804
[[ 0.11725372 -1.11016952]]
[1.0096957]
C = 0.1, penalty = l2 score is  0.808
[[ 0.08271496 -1.11618329]]
[1.00965503]
C = 1, penalty = l1 score is  0.808

通过不同参数得到的特征系数不同且预测的准确度也不同。但都在0.8左右,因为使用的是一次项表达式,而模型是2次项,因此使用多项式逻辑回归准确度应该会好些。

def polynomialLogisiticRegression(degree, C=1.0, penalty='l2'):
    '''
    多项式逻辑回归
    :param degree: 
    :param C: 
    :param penalty: 
    :return: 
    '''
    return Pipeline([
        ('poly', PolynomialFeatures(degree=degree)),
        ('std', StandardScaler()),
        ('reg', LogisticRegression(C=C, penalty=penalty)),
    ])

poly_reg = polynomialLogisiticRegression(degree=2)
poly_reg.fit(X_train, y_train)
print('the score of polynomial logisitic regression with degree=2, C=1, penalty=l2 is ', poly_reg.score(X, y))

poly_reg = polynomialLogisiticRegression(degree=20)
poly_reg.fit(X_train, y_train)
print('the score of polynomial logisitic regression with degree=20, C=1, penalty=l2 is ', poly_reg.score(X, y))

poly_reg = polynomialLogisiticRegression(degree=20, C=0.1)
poly_reg.fit(X_train, y_train)
print('the score of polynomial logisitic regression with degree=20, C=0.1, penalty=l2 is ', poly_reg.score(X, y))

poly_reg = polynomialLogisiticRegression(degree=20, C=0.1, penalty='l1')
poly_reg.fit(X_train, y_train)
print('the score of polynomial logisitic regression with degree=20, C=0.1, penalty=l1 is ', poly_reg.score(X, y))

使用PipeLine封装多项式回归和逻辑回归,通过将项数为2和20进行测试,很明显,20对于训练数据而言是过拟合的,因此,在项数20的基础上进行L1正则和L2正则。
得到输出结果

the score of polynomial logisitic regression with degree=2, C=1, penalty=l2 is  0.991
the score of polynomial logisitic regression with degree=20, C=1, penalty=l2 is  0.9905
the score of polynomial logisitic regression with degree=20, C=0.1, penalty=l2 is  0.9765
the score of polynomial logisitic regression with degree=20, C=0.1, penalty=l1 is  0.9935

可以对比出当项数过高时,模型的准确度略微降低,当使用L1正则并设置参数C=0.1后模型的准确度提升,但整体上高于之前一次项回归。可以看出正则化和多项式回归的效果。

OvR与OvO

逻辑回归本身只能解决二分类问题,但通过一些手段可以解决多分类问题。常见的方式有OvR和OvO两种。

假设现在有四种类别,分别以不同颜色表示。
①OvR的思路

OvR(One vs Rest)的思路是将所有类别分为两个类,当前类别是一类,其他类别合并视为一个类。

那么有K个类别的数据样本就会被分为K个由两个类组成的新样本集合,这样就将多分类问题转化为二分类问题,然后对每个数据样本进行模型训练,得到模型使用样本进行验证,选择分类得分最高的作为最终的样本类别。这里的得分最高指的就是样本属于某类概率最高。
sklearn提供了OvR的实现,可以在分类器中进行实现,也可以使用OvR类,以分类器作为参数进行实现。

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

#以鸢尾花数据为例
iris = dataset.load_iris()
X = iris.data
y = iris.target
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=666)

#打印LogisticRegression实例信息所知,默认采用OvR方式
log_reg = LogisticRegression()
log_reg.fit(X_train, y_train)
print(log_reg)
score = log_reg.score(X_test, y_test)
print('with multi class OvR, the score is ', score)

#使用sklearn的OvR分类器
ovr = OneVsRestClassifier(log_reg)
ovr.fit(X_train, y_train)
score = ovr.score(X_test, y_test)
print('use sklearn OneVsRestClassifier, the score is ', score)

得到输出

LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,
          intercept_scaling=1, max_iter=100, multi_class='ovr', n_jobs=1,
          penalty='l2', random_state=None, solver='liblinear', tol=0.0001,
          verbose=0, warm_start=False)
with multi class OvR, the score is  0.9473684210526315
use sklearn OneVsRestClassifier, the score is  0.9473684210526315

LogisticRegression默认使用OvR的方式进行分类,可以通过设置multiclass参数调整。在sklearn的multiclass包下也对OvR和OvO进行封装,可以将分类器对象作为参数传入,对象可以是sklearn自带的,也可以是自定义的,只要满足函数规范即含有fit、predict、score等函数即可实现分类。
使用鸢尾花数据集,解决三分类问题,OvR模式下将生成三个分类模型,分类的准确度为94%。

②OvO的思路

OvO即One vs One,不再将出当前类别外的其他类别视为一个大类,而是将类别拆分,每两个训练一个模型。例如有四个类别的样本将会被分为如下6个OvO类别。

多分类问题转化为两个两个类别进行分类的问题,一个含有N个类别的样本将被分为C(N, 2)个分类问题。当对一个未知样本进行分类时,用这C(N, 2)个模型进行分类,最后得票最多的类别即为该未知样本的类别。
因为OvO产生的模型数量高于OvR,因此,一般情况下,OvO的分类准确度好于OvR,但训练更多模型的系统和时间开销要大,因此OvR更加高效。

#使用OvO进行分类
log_reg2 = LogisticRegression(multi_class='multinomial', solver='newton-cg')
log_reg2.fit(X_train, y_train)
score = log_reg2.score(X_test, y_test)
print('with multi class OvO, the score is ', score)

#使用sklearn的OvO分类器
ovo = OneVsOneClassifier(log_reg)
ovo.fit(X_train, y_train)
score = ovo.score(X_test, y_test)
print('use sklearn OneVsOneClassifier, the score is ', score)

如果是使用LogisticRegression进行OvO训练,除了设置multi_class参数为multinomial外,还需要指定solver参数,否则无效。若使用OneVsOneClassifier,使用方式和OneVsRestClassifier没有什么差别。
得到输出结果

with multi class OvO, the score is  1.0
use sklearn OneVsOneClassifier, the score is  1.0

分类结果100%正确,可见在当前样本下,OvO确实好于OvR。

最后说明,OvO和OvR不是逻辑回归专有的,也可以运用与SVM等只能解决二分类问题的模型上。

推荐阅读更多精彩内容