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

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

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越大，损失函数值就越大，道理同上。

给出逻辑回归的损失函数

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

将p代入代价函数J

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

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

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也能打印得到。

决策边界

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

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的思路

那么有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的思路

多分类问题转化为两个两个类别进行分类的问题，一个含有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等只能解决二分类问题的模型上。