機器學習——從線性回歸到邏輯回歸【附詳細推導和代碼】

本文始發於個人公眾號：TechFlow，原創不易，求個關注

在之前的文章當中，我們推導了線性回歸的公式，線性回歸本質是線性函數，模型的原理不難，核心是求解模型參數的過程。通過對線性回歸的推導和學習，我們基本上了解了機器學習模型學習的過程，這是機器學習的精髓，要比單個模型的原理重要得多。

新關注和有所遺忘的同學可以點擊下方的鏈接回顧一下之前的線性回歸和梯度下降的內容。

一文講透梯度下降法

詳細推導線性回歸模型

回歸與分類

在機器學習當中，模型根據預測結果的不同分為兩類，如果我們希望模型預測一個或者多個連續值，這類問題被稱為是回歸問題。像是常見的未來股票價格的估計、未來溫度估計等等都算是回歸問題。還有一類呢是分類問題，模型的預測結果是一個離散值，也就是說只有固定的種類的結果。常見的有垃圾郵件識別、網頁圖片鑒黃等等。

我們之前介紹的邏輯回歸顧名思義是一個回歸問題，今天的文章講的呢是如何將這個回歸模型轉化成分類模型，這個由線性回歸推導得到的分類模型稱為邏輯回歸。

邏輯回歸

邏輯回歸這個模型很神奇，雖然它的本質也是回歸，但是它是一個分類模型，並且它的名字當中又包含”回歸“兩個字，未免讓人覺得莫名其妙。

如果是初學者，覺得頭暈是正常的，沒關系，讓我們一點點捋清楚。

讓我們先回到線性回歸，我們都知道，線性回歸當中\(y=WX+b\)。我們通過W和b可以求出X對應的y，這里的y是一個連續值，是回歸模型對吧。但如果我們希望這個模型來做分類呢，應該怎么辦？很容易想到，我們可以人為地設置閾值對吧，比如我們規定y > 0最后的分類是1，y < 0最后的分類是0。從表面上來看，這當然是可以的，但實際上這樣操作會有很多問題。

最大的問題在於如果我們簡單地設計一個閾值來做判斷，那么會導致最后的y是一個分段函數，而分段函數不連續，使得我們沒有辦法對它求梯度，為了解決這個問題，我們得找到一個平滑的函數使得既可以用來做分類，又可以解決梯度的問題。

很快，信息學家們找到了這樣一個函數，它就是Sigmoid函數，它的表達式是：

\[\sigma(x)= \frac{1}{1+e^{-x}} \]

它的函數圖像如下：

可以看到，sigmoid函數在x=0處取值0.5，在正無窮處極限是1，在負無窮處極限是0，並且函數連續，處處可導。sigmoid的函數值的取值范圍是0-1，非常適合用來反映一個事物發生的概率。我們認為\(\sigma(x)\)表示x發生的概率，那么x不發生的概率就是\(1-\sigma(x)\)。我們把發生和不發生看成是兩個類別，那么sigmoid函數就轉化成了分類函數，如果\(\sigma(x)>0.5\)表示類別1，否則表示類別0.

到這里就很簡單了，通過線性回歸我們可以得到\(y = WX + b, P(1)=\sigma(y), P(0)=1-\sigma(y)\)。也就是說我們在線性回歸模型的外面套了一層sigmoid函數，我們通過計算出不同的y，從而獲得不同的概率，最后得到不同的分類結果。

損失函數

下面的推導全程高能，我相信你們看完會三連的（點贊、轉發、關注）。

讓我們開始吧，我們先來確定一下符號，為了區分，我們把訓練樣本當中的真實分類命名為\(y\)，\(y\)的矩陣寫成\(Y\)。同樣，單條樣本寫成\(x\),\(x\)的矩陣寫成\(X\)。單條預測的結果寫成\(\hat{y}\)，所有的預測結果寫成\(\hat{Y}\)。

對於單條樣本來說，y有兩個取值，可能是1，也可能是0，1和0代表兩個不同的分類。我們希望\(y = 1\)的時候，\(\hat{y}\)盡量大，\(y = 0\)時，\(1 - \hat{y}\)盡量大，也就是\(\hat{y}\)盡量小，因為它取值在0-1之間。我們用一個式子來統一這兩種情況：

\[P(y|x)=\hat{y}^y(1-\hat{y})^{1-y} \]

我們代入一下，\(y=0\)時前項為1，表達式就只剩下后項，同理，\(y=1\)時，后項為1，只剩下前項。所以這個式子就可以表示預測准確的概率，我們希望這個概率盡量大。顯然，\(P(y|x)>0\)，所以我們可以對它求對數，因為log函數是單調的。所以\(P(y|x)\)取最值時的取值，就是\(logP(y|x)\)取最值的取值。

\[\log P(y|x)=y\log \hat{y} + (1-y)\log (1-\hat{y}) \]

我們期望這個值最大，也就是期望它的相反數最小，我們令\(J=-\log P(y|x)\)，這樣就得到了它的損失函數：

\[J(\theta)=-\frac{1}{m}Y\log \hat{Y} + (1-Y)\log (1-\hat{Y}) \]

如果知道交叉熵這個概念的同學，會發現這個損失函數的表達式其實就是交叉熵。交叉熵是用來衡量兩個概率分布之間的”距離“，交叉熵越小說明兩個概率分布越接近，所以經常被用來當做分類模型的損失函數。關於交叉熵的概念我們這里不多贅述，會在之后文章當中詳細介紹。我們隨手推導的損失函數剛好就是交叉熵，這並不是巧合，其實底層是有一套信息論的數學邏輯支撐的，我們不多做延伸，感興趣的同學可以了解一下。

硬核推導

損失函數有了，接下來就是求梯度來實現梯度下降了。

這個函數看起來非常復雜，要對它直接求偏導過於硬核（危），如果是許久不碰高數的同學直接肝不亞於硬抗葦名一心。

為了簡化難度，我們先來做一些准備工作。首先，我們先來看下\(\sigma\)函數，它本身的形式很復雜，我們先把它的導數搞定。

\[\begin{aligned} (\frac{1}{1+e^{-x}})'&=\frac{e^{-x}}{(1+e^{-x})^2} \\ &=\frac{1+e^{-x}-1}{(1+e^{-x})^2} \\ &=\frac{1}{1+e^{-x}}\cdot (1-\frac{1}{1+e^{-x}}) \\ &=\sigma(x) (1-\sigma(x)) \end{aligned} \]

因為\(\hat{y}=\sigma(\theta X)\)，我們將它帶入損失函數，可以得到，其中\(\sigma(\theta)=\sigma(\theta X)\)：

\[J(\theta)=-\frac{1}{m}Y\log \sigma(\theta) + (1 - Y)\log( 1- \sigma(\theta)) \]

接着我們求\(J(\theta)\)對\(\theta\)的偏導，這里要代入上面對\(\sigma(x)\)求導的結論：

\[\begin{aligned} \frac{\partial}{\partial \theta}J(\theta)&=-\frac{1}{m}(Y\cdot \frac{1}{\sigma(\theta)}\cdot \sigma'(\theta)+(1-Y)\cdot \frac{1}{1-\sigma(\theta)}\cdot \sigma'(\theta)\cdot(-1)) \\ &=-\frac{1}{m}(Y\cdot \frac{1}{\sigma(\theta)}\cdot \sigma(\theta)(1-\sigma(\theta))\cdot X-(1-Y)\cdot \frac{1}{1-\sigma(\theta)}\cdot \sigma(\theta)(1-\sigma(\theta))\cdot X) \\ &=-\frac{1}{m}(Y(1-\sigma(\theta)) - \sigma(\theta)+Y\cdot \sigma(\theta))\cdot X\\ &=-\frac{1}{m}(\sigma(\theta) - Y)\cdot X\\ &=-\frac{1}{m}(\hat{Y} - Y)\cdot X \end{aligned} \]

代碼實戰

梯度的公式都推出來了，離寫代碼實現還遠嗎？

不過巧婦難為無米之炊，在我們擼模型之前，我們先試着造一批數據。

我們選擇生活中一個很簡單的場景——考試。假設每個學生需要參加兩門考試，兩門考試的成績相加得到最終成績，我們有一批學生是否合格的數據。希望設計一個邏輯回歸模型，幫助我們直接計算學生是否合格。

為了防止sigmoid函數產生偏差，我們把每門課的成績縮放到(0, 1)的區間內。兩門課成績相加超過140分就認為總體及格。

import numpy as np
def load_data():
    X1 = np.random.rand(50, 1)/2 + 0.5
    X2 = np.random.rand(50, 1)/2 + 0.5
    y = X1 + X2
    Y = y > 1.4
    Y = Y + 0
    x = np.c_[np.ones(len(Y)), X1, X2]
    return x, Y

這樣得到的訓練數據有兩個特征，分別是學生兩門課的成績，還有一個偏移量1，用來記錄常數的偏移量。

接着，根據上文當中的公式，我們不難（真的不難）實現sigmoid以及梯度下降的函數。

def sigmoid(x):
    return 1.0/(1+np.exp(-x))
    
def gradAscent(data, classLabels):
    X = np.mat(data)
    Y = np.mat(classLabels)
    m, n = np.shape(dataMat)
    # 學習率
    alpha = 0.01
    iterations = 10000
    wei = np.ones((n, 1))
    for i in range(iterations):
        y_hat = sigmoid(X * wei)
        error = Y - y_hat
        wei = wei + alpha * X.transpose()*error
    return wei

這段函數實現的是批量梯度下降，對Numpy熟悉的同學可以看得出來，這就是在直接套公式。

最后，我們把數據集以及邏輯回歸的分割線繪制出來。

import matplotlib.pyplot as plt
def plotPic():
    dataMat, labelMat = load_data()
    #獲取梯度
    weis = gradAscent(dataMat, labelMat).getA()
    n = np.shape(dataMat)[0]
    xc1 = []
    yc1 = []
    xc2 = []
    yc2 = []
    #數據按照正負類別划分
    for i in range(n):
        if int(labelMat[i]) == 1:
            xc1.append(dataMat[i][1])
            yc1.append(dataMat[i][2])
        else:
            xc2.append(dataMat[i][1])
            yc2.append(dataMat[i][2])
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(xc1, yc1, s=30, c='red', marker='s')
    ax.scatter(xc2, yc2, s=30, c='green')
    #繪制分割線
    x = np.arange(0.5, 1, 0.1)
    y = (-weis[0]- weis[1]*x)/weis[2]
    ax.plot(x, y)
    plt.xlabel('X1')
    plt.ylabel('X2')
    plt.show()

最后得到的結果如下：

隨機梯度下降版本

可以發現，經過了1萬次的迭代，我們得到的模型已經可以正確識別所有的樣本了。

我們剛剛實現的是全量梯度下降算法，我們還可以利用隨機梯度下降來進行優化。優化也非常簡單，我們計算梯度的時候不再是針對全量的數據，而是從數據集中選擇一條進行梯度計算。

基本上可以復用梯度下降的代碼，只需要對樣本選取的部分加入優化。

def stocGradAscent(data, classLab, iter_cnt):
    X = np.mat(data)
    Y = np.mat(classLab)
    m, n = np.shape(dataMat)
    weis = np.ones((n, 1))
    for j in range(iter_cnt):
        for i in range(m):
            # 使用反比例函數優化學習率
            alpha = 4 / (1.0 + j + i) + 0.01
            randIdx = int(np.random.uniform(0, m))
            y_hat = sigmoid(np.sum(X[randIdx] * weis))
            error = Y[randIdx] - y_hat
            weis = weis + alpha * X[randIdx].transpose()*error
    return weis

我們設置迭代次數為2000，最后得到的分隔圖像結果如下：

當然上面的代碼並不完美，只是一個簡單的demo，還有很多改進和優化的空間。只是作為一個例子，讓大家直觀感受一下：其實自己親手寫模型並不難，公式的推導也很有意思。這也是為什么我會設置高數專題的原因。CS的很多知識也是想通的，在學習的過程當中靈感迸發旁征博引真的是非常有樂趣的事情，希望大家也都能找到自己的樂趣。

今天的文章就是這些，如果覺得有所收獲，請順手掃碼點個關注吧，你們的舉手之勞對我來說很重要。