邏輯迴歸

邏輯迴歸是一種解決分類問題的機器學習演算法。

邏輯迴歸的思想是將樣本特徵和樣本發生的機率聯絡起來，機率是一個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 npfrom sklearn.metrics import mean_squared_errorimport sklearn.datasets as datasetfrom sklearn.model_selection import train_test_splitclass 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.0the 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.0007333333333333333regression is 0.0004486610126586672 * X1 + 0.0001887681532517242 * X2 + 0.0007333333333333333the 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 npimport matplotlib.pyplot as pltfrom sklearn.model_selection import train_test_splitfrom sklearn.linear_model import LogisticRegressionfrom sklearn.preprocessing import StandardScalerfrom sklearn.pipeline import Pipelinefrom 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個點值為1for _ in range(20):
    y[np.random.randint(200)] = 1X_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.991the score of polynomial logisitic regression with degree=20, C=1, penalty=l2 is  0.9905the score of polynomial logisitic regression with degree=20, C=0.1, penalty=l2 is  0.9765the 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 npimport sklearn.datasets as datasetimport matplotlib.pyplot as pltfrom sklearn.model_selection import train_test_splitfrom sklearn.linear_model import LogisticRegressionfrom sklearn.multiclass import OneVsRestClassifierfrom 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.9473684210526315use 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.0use sklearn OneVsOneClassifier, the score is  1.0

分類結果100%正確，可見在當前樣本下，OvO確實好於OvR。

最後說明，OvO和OvR不是邏輯迴歸專有的，也可以運用與SVM等只能解決二分類問題的模型上。

作者：Chuck_Hu
連結：