天天看點
傳回

邏輯回歸

邏輯回歸

目錄

    • 一、認識Logistic Regression
    • 二、Logistic Regression的損失函數
    • 三、Logistic Regression的損失函數的梯度
    • 四、程式設計實作Logistic Regression
    • 五、決策邊界
      • 1、Logistic Regression的決策邊界
      • 2、KNN的決策邊界
    • 七、邏輯回歸中使用多項式特征
    • 八、邏輯回歸中使用正則化
      • 1、使用Logistic Regression L2正則
      • 2、使用Logistic Regression L1正則
    • 九、邏輯回歸解決多分類問題
      • 1、OvR
      • 2、OvO
      • 3、Logistic Regression的OvR和OvO的程式設計實作
      • 4、sklearn中的OvR和OvO
    • 我是尾巴

邏輯回歸（Logistics Regression）,用來解決分類問題，那回歸怎麼解決分類問題？

将樣本特征和樣本發生的機率聯系起來，機率是一個數。

這是一個統計資料，Logistic Regression是最廣泛使用的一種算法。

邏輯回歸通常既可以看做回歸算法，又可以看做分類算法，通常作為分類算法，隻可以解決二分類問題。

通常我們在直線回歸過程中使用第一行的公式，但是他的值域是從（-infineity， +infinity）而所需的機率的值域為[0,1],是以做一下改進，

接下裡繪制一下這個函數的圖形：

import numpy as np
import matplotlib.pyplot as plt

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

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

那麼問題來了：

對于給定的樣本資料集x，y如何找到參數\theta，使得使用這樣的方式，可以最大程度地獲得樣本資料集x，對應的分類輸出y？

這個函數是一個凸函數，他沒有公式解，隻能使用梯度下降法求解。

對于這種複雜的損失函數，我們首先來看一個sigmoid函數的求導。

import numpy as np
from play_ML.multi_linear_regression.metrics import accuracy_score

class LogisticRegression(object):

    def __init__(self):
        "初始化logistic回歸模型"
        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):
        assert x_train.shape[0] == y_train.shape[0], "the size of x_train must be equal to the size of y_train"

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

        def dj(theta, X_b, y):
            return X_b.T.dot(self._sigmoid(X_b.dot(theta)) - y) / len(X_b)

        def gradient_descent(X_b, y, init_theta, eta, n_iters=1e4, epsilon=1e-8):

            theta = init_theta
            i_iter = 0

            while i_iter < n_iters:
                gradient = dj(theta, X_b, y)
                last_theta = theta
                theta = last_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])
        init_theta = np.zeros(X_b.shape[1])
        self._theta = gradient_descent(X_b, y_train, init_theta, eta, n_iters=1e4)
        self.interception_ = self._theta[0]
        self.coef_ = self._theta[1:]

        return self

    def predict_prob(self, x_predict):
        assert self.interception_ is not None and self.coef_ is not None, "must fit before predict"
        assert x_predict.shape[1] == len(self.coef_), "the feature number must be equal to x_train"

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

    def predict(self, x_predict):
        assert self.interception_ is not None and self.coef_ is not None, "must fit before predict"
        assert x_predict.shape[1] == len(self.coef_), "the feature number must be equal to x_train"

        prob = self.predict_prob(x_predict)
        return np.array(prob >= 0.5, dtype='int')

    def score(self, x_test, y_test):
        y_preict = self.predict(x_test)
        return accuracy_score(y_test, y_preict)

    def __repr__(self):
        return "Logistic Regression"

測試一下：

if __name__ == '__main__':
    import numpy as np
    from sklearn import datasets
    import matplotlib.pyplot as plt
    from sklearn.model_selection import train_test_split

    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()

    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))
    print(log_reg.predict_prob(x_test))
    print(y_test)
    print(log_reg.coef_)
    print(log_reg.interception_)

輸出結果：

1.0
# 由此可以看出Logistic的簡單的鸢尾花資料集具有良好的表現。後續會用複雜的資料集盡心測試。
# 下面分别是Logistic預測出來的結果，我們以0.5将資料分為兩類，通過對比真實标簽可以看出預測與實際的差距
[0.93292947 0.98717455 0.15541379 0.01786837 0.03909442 0.01972689
 0.05214631 0.99683149 0.98092348 0.75469962 0.0473811  0.00362352
 0.27122595 0.03909442 0.84902103 0.80627393 0.83574223 0.33477608
 0.06921637 0.21582553 0.0240109  0.1836441  0.98092348 0.98947619
 0.08342411]
[1 1 0 0 0 0 0 1 1 1 0 0 0 0 1 1 1 0 0 0 0 0 1 1 0]
# 其實Logistic回歸算法是從線性回歸改進而來的，同樣我們得到類似相性回歸的系數和截距，那它有什麼幾何意義呢？其實就是決策邊界。cc
[ 3.01749692 -5.03046934]
-0.6827383698993108

首先回顧一下Logistic Regression進行分類的原理：

接下來再看一下Sigmoid這個函數

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()

第四節預測明明是100%的準确率為什麼還會有個紅點在決策邊界下方呢？因為這是訓練集，接下裡試一下測試集：

那就把這個繪制決策邊界封裝成一個函數，

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)
    
    from matplotlib.colors import ListedColormap
    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作用在兩個類别上： 
from sklearn.neighbors import KNeighborsClassifier

knn_clf = KNeighborsClassifier()
knn_clf.fit(x_train, y_train)

knn_clf.score(x_test, y_test)
# 1.0
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作用在三個類别上： 
from sklearn.neighbors import KNeighborsClassifier

knn_clf_all = KNeighborsClassifier()
knn_clf_all.fit(iris.data[:,:2], iris.target)
#KNeighborsClassifier(algorithm='auto', leaf_size=30, metric='minkowski',
#          metric_params=None, n_jobs=None, n_neighbors=5, p=2,
#          weights='uniform')

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='blue')
plt.show()

通過上面的圖可以發現KNN其實已經過拟合了，這是因為在執行個體化KNN的時候預設使用的k=5，k越小模型就越複雜,k越大模型越簡單。複雜的表現就是決策邊界不規整，那麼把k調大試一下。

from sklearn.neighbors import KNeighborsClassifier

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='blue')
plt.show()

​ 邏輯回歸中如果使用直線分類方式就隻能針對二分類了，如果像下圖中不可能使用一根直線完成分割，但是很顯然可以使用圓形或者橢圓形完整這個分類任務。其實線上性回歸到多項式回歸我們思想就是給訓練資料集添加多項式項。同理我們把這個東西用到邏輯回歸中。

首先生成需要的資料：

import numpy as np
import matplotlib.pyplot as plt

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)
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()

顯然有很多錯誤的分類，那接下來給邏輯回歸添加多項式項。

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures
from sklearn.preprocessing import StandardScaler

def PolynomiaLogisticRegression(degree):
    return Pipeline([
        ('poly', PolynomialFeatures(degree=degree)),
        ('std_scale', StandardScaler()),
        ('log_reg', LogisticRegression())
    ])

poly_log_reg = PolynomiaLogisticRegression(degree=2)
poly_log_reg.fit(x, y)
poly_log_reg.score(x, y)
# 0.94999999999999996
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()

那接下裡把degree調大，試一下

poly_log_reg2 = PolynomiaLogisticRegression(degree=20)
poly_log_reg2.fit(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()

出現這樣的邊界形狀是因為發生了過拟合現象，degree越來模型越複雜，也就越容易發生過拟合。接下來我們就來解決這個過拟合問題。這裡使用模型正則化。

​ 在實際的應用過程中，很少有問題直接用直線就能完成分類或者回歸任務，是以正則化必不可少。之前學過模型泛化的時候提到的L1正則、L2正則化的方式：左邊

但是在sklearn中對邏輯回歸中的正則化：右邊

  • 先使用直線Logistic Regression：

    接下來看看sklearn中的邏輯回歸是如何加入正則化的。還是先生成樣本資料：

import numpy as np
import matplotlib.pyplot as plt

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 i 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()

           
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression

x_train, x_test, y_train, y_test = train_test_split(x, y)
log_reg = LogisticRegression()
log_reg.fit(x, y)
# LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,
#          intercept_scaling=1, max_iter=100, multi_class='warn',
#         n_jobs=None, penalty='l2', random_state=None, solver='warn',
#          tol=0.0001, verbose=0, warm_start=False)

​ 通過輸出結果我們可以發現預設C=1.0，這個C就是最開始提到的邏輯回歸中正則中的C，penalty='l2'說明sklearn預設使用L2正則來進行模型正則化。

log_reg.score(x_train, y_train)
# 0.7933333333333333
log_reg.score(x_test, y_test)
# 0.7933333333333333
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()
  • 使用多項式Logistic Regression 
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures
from sklearn.preprocessing import StandardScaler

def PolynomiaLogisticRegression(degree):
    return Pipeline([
        ('poly', PolynomialFeatures(degree=degree)),
        ('std_scale', StandardScaler()),
        ('log_reg', LogisticRegression())
    ])

poly_log_reg = PolynomiaLogisticRegression(degree=2)
poly_log_reg.fit(x_train, y_train)
poly_log_reg.score(x_train, y_train)
# 0.9133333333333333
poly_log_reg.score(x_test, y_test)
# 0.94
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 = PolynomiaLogisticRegression(degree=20)
poly_log_reg2.fit(x_train, y_train)

poly_log_reg2.score(x_train, y_train)
# 0.94
poly_log_reg2.score(x_test, y_test)
# 0.92
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()

           

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures
from sklearn.preprocessing import StandardScaler

# 傳入一個新的參數C
def PolynomiaLogisticRegression(degree, C):
    return Pipeline([
        ('poly', PolynomialFeatures(degree=degree)),
        ('std_scale', StandardScaler()),
        ('log_reg', LogisticRegression(C=C))
    ])

poly_log_reg3 = PolynomiaLogisticRegression(degree=20, C=0.1)
poly_log_reg3.fit(x, y)

poly_log_reg3.score(x_train, y_train)
# 0.8533333333333334
poly_log_reg3.score(x_test, y_test)
# 0.92
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()

           

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures
from sklearn.preprocessing import StandardScaler

def PolynomiaLogisticRegression(degree, C, penalty='l2'):
    return Pipeline([
        ('poly', PolynomialFeatures(degree=degree)),
        ('std_scale', StandardScaler()),
        ('log_reg', LogisticRegression(C=C, penalty=penalty))
    ])

poly_log_reg4 = PolynomiaLogisticRegression(degree=20, C=0.1, penalty='l1')
poly_log_reg4.fit(x_train, y_train)

poly_log_reg4.score(x_train, y_train)
# 0.8266666666666667
poly_log_reg4.score(x_test, y_test)
# 0.9
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()

通過以上的例子，可以大緻了解sklearn中已經封裝了正則化的内容，但是在實際問題中我們并不知道degree，C，penalty這些超參數，是以就需要進行網格搜尋進行尋優。

​ 在開始之初，說邏輯回歸隻可以解決二分類問題， 其實可以稍加改造使其能夠解決多分類問題。當然這個改造方式并不是隻針對邏輯回歸這一種算法，這是一種通用的近乎于可以改造所有的二分類。

​ （One vs Rest）一對剩餘，有些說法也叫（One vs All，OVA）。比如下圖中的這個四分類問題，顯然如果使用邏輯回歸一條直線分出4類是不現實的，但是如果我們取出其中任意一種，将剩下的作為另一種，這種就是一個2分類問題，同理将每一個類别分别做一次這樣的2分類，如果有n個類别就進行n次分類，選擇分類得分最高的。就類似于下圖這種，進行C（n，1）中分類。進而完成多分類。

​ （One vs One）一對一，就是在多個類别中，先挑出2個來進行2分類，然後逐個進行，也就是C（n，2）中情況進行2分類，選擇赢數最高的分類。比如一個手寫數字的識别任務來說就要進行C（10,2）=45次分類，才能完成任務。很顯然它相比OvR（n級别）來說，OvO這是一個n[^2]級别的，要比耗時更多，但是它分類的結果更加準确。這是因為每次都是在用兩個真實的類别再進行2分類，而沒有混淆其他的類别資訊。

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

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

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)
# LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,
#          intercept_scaling=1, max_iter=100, multi_class='warn',
#          n_jobs=None, penalty='l2', random_state=None, solver='warn',
#          tol=0.0001, verbose=0, warm_start=False)

從輸出結果來看，這個multi_class='warn',這個multi_class就是多分類的意思，從官方文檔來看預設是OvR，不知道這個warn是什麼意思。先不管，把這個參數傳進去就是了。想具體了解每個參數的意義，請檢視官方文檔：from sklearn.linear_model import LogisticRegression

solver這個參數是因為在sklearn中并不是簡單地使用梯度下降，是以我們需要給不同的方法傳入不同的解決辦法。

  • OvR 
log_reg = LogisticRegression(multi_class='ovr',solver='liblinear')
log_reg.fit(x_train, y_train)
log_reg.score(x_test, y_test)
# 0.6578947368421053
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)
# 0.7894736842105263
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()

通過上面這兩個示例發現準确率并不高，這是因為鸢尾花資料集共有4個特征，而隻用了前兩個,為了友善可視化，下面就是使用所有資料。

iris = datasets.load_iris()
x = iris.data[:, :]
y = iris.target

x_train, x_test, y_train, y_test = train_test_split(x, y, random_state=666)

# OvR
log_reg3 = LogisticRegression(multi_class='ovr', solver='liblinear')
log_reg3.fit(x_train, y_train)
log_reg3.score(x_test, y_test)
# 0.9473684210526315

# OvO
log_reg4 = LogisticRegression(multi_class='multinomial', solver='newton-cg')
log_reg4.fit(x_train, y_train)
log_reg4.score(x_test, y_test)
# 1.0

           

from sklearn.linear_model import LogisticRegression

log_reg = LogisticRegression()

from sklearn.multiclass import OneVsRestClassifier

ovr = OneVsRestClassifier(log_reg)
ovr.fit(x_train, y_train)
# OneVsRestClassifier(estimator=LogisticRegression(C=1.0, class_weight=None, 
# dual=False, fit_intercept=True,
#          intercept_scaling=1, max_iter=100, multi_class='ovr',
#          n_jobs=None, penalty='l2', random_state=None, solver='liblinear',
#          tol=0.0001, verbose=0, warm_start=False),
#          n_jobs=None)
ovr.score(x_test, y_test)
# 0.9473684210526315

from sklearn.multiclass import OneVsOneClassifier

ovo = OneVsOneClassifier(log_reg)
ovo.fit(x_train, y_train)
# OneVsOneClassifier(estimator=LogisticRegression(C=1.0, class_weight=None, 
# dual=False, fit_intercept=True,
#          intercept_scaling=1, max_iter=100, multi_class='multinomial',
#          n_jobs=None, penalty='l2', random_state=None, solver='newton-cg',
#          tol=0.0001, verbose=0, warm_start=False),
#          n_jobs=None)

本節就先到這。

本次推薦：

一個幫助開發者成長的社群

毒雞湯：孩子，窮怎麼了？窮也要挺起胸膛來，讓别人看看，你不僅窮而且還矮。

矮又如何？擡起你的頭來，讓他們知道，你不僅矮，而且還醜。

醜不要緊，用你的言談舉止讓其他人明白，你還是一個沒有内涵的人。

沒有内涵也不要放棄，從現在開始學習。當你讀了足夠多的書的時候，你會發現自己還笨。

繼續加油~

堅持到無能為力，拼搏到感動自己！

邏輯回歸 機器學習
上一篇: 模型泛化
下一篇: Mysql高性能優化規範建議

繼續閱讀