第5章 Logistic迴歸

假設現在有一些資料點，我們用一條直線對這些點進行擬合（該線稱為最佳擬合直線），這個擬合過程就稱作迴歸。利用Logistic迴歸進行分類的主要思想是：根據現有資料對分類邊界線建立迴歸公式，以此進行分類。

Logistic迴歸的一般流程

1. 收集資料：可以使用任何方法。
2. 準備資料：由於需要進行距離計算，因此要求資料型別為數值型。另外，結構化資料格式則最佳。
3. 分析資料：採用任意方法對資料進行分析。
4. 訓練演算法：大部分時間將用於訓練，訓練的目的是為了找到最佳的分類迴歸係數。
5. 測試演算法：一旦訓練步驟完成，分類將會很快。
6. 使用演算法：首先，我們需要輸入一些資料，並將其轉換成對應的結構化數值；接著，基於訓練好的迴歸係數就可以對這些數值進行簡單的迴歸計算，判定它們屬於哪個類別；在這之後，我們就可以在輸出的類別上做一些其他分析工作。

5.1 基於 Logistic 迴歸和 Sigmoid 函式的分類

Logistic迴歸：

• 優點：計算代價不高，易於理解和實現。
• 缺點：容易欠擬合，分類精度可能不高。
• 適用資料型別：數值型和標稱型資料。

我們想要的函式應該是，能接受所有的輸入然後預測出類別。例如，在兩個類的情況下，上述函式輸出0或1。Sigmoid函式具體的計算公式如下：
$\sigma(z)=\frac{1}{1+e^{-z}}$

當x為0時，Sigmoid函式值為0.5。隨著x的增大，對應的Sigmoid值將逼近於1；而隨著x的減小，Sigmoid值將逼近於0。

import matplotlib.
 
pyplot as plt
import numpy as np
import pandas as pd

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

x = np.linspace(-10, 10, 100)
plt.title("Logical Function")
plt.xlabel('x')
plt.ylabel('y')
plt.xlim(-12, 12)
plt.ylim(0.0, 1.1)
plt.plot(x, sigmoid(x))
plt.plot([-12, 12], [0.5, 0.5])
plt.plot([0, 0], [0, 1])
 

plt.show()

為了實現Logistic迴歸分類器，我們可以在每個特徵上都乘以一個迴歸係數，然後把所有的結果值相加，將這個總和代入Sigmoid函式中，進而得到一個範圍在0~1之間的數值。任何大於0.5的資料被分入1類，小於0.5即被歸入0類。所以，Logistic迴歸也可以被看成是一種概率估計。

5.2 基於最優化方法的最佳迴歸係數確定

Sigmoid函式的輸入記為z，由下面公式得出：
$z=w_{0}x_{0}+w_{1}x_{1}+w_{2}x_{2}+\cdots+w_{n}x_{n}$

如果採用向量的寫法，上述公式可以寫成$z=w^{T}x$，它表示將這兩個數值向量對應元素相乘然後全部加起來即得到z值。其中的向量x是分類器的輸入資料，向量w也就是我們要找到的最佳引數（係數），從而使得分類器儘可能地精確。為了尋找該最佳引數，需要用到最優化理論的一些知識。

5.2.1 梯度上升法

梯度上升法基於的思想是：要找到某函式的最大值，最好的方法是沿著該函式的梯度方向探尋。如果梯度記為∇，則函式f(x,y)的梯度由下式表示：
$\nabla f(x,y)= \left(\begin{matrix} \frac{df(x,y)}{dx} \\ \frac{df(x,y)}{dy} \end{matrix} \right)$

這是機器學習中最易造成混淆的一個地方，但在數學上並不難，需要做的只是牢記這些符號的意義。這個梯度意味著要沿x的方向移動$\frac{df(x,y)}{dx}$，沿y的方向移動$\frac{df(x,y)}{dy}$。其中，函式f(x,y)必須要在待計算的點上有定義並且可微。

梯度上升演算法到達每個點後都會重新估計移動的方向。從P0開始，計算完該點的梯度，函式就根據梯度移動到下一點P1。在P1點，梯度再次被重新計算，並沿新的梯度方向移動到P2。如此迴圈迭代，直到滿足停止條件。迭代的過程中，梯度運算元總是保證我們能選取到最佳的移動方向。

梯度上升演算法沿梯度方向移動了一步。可以看到，梯度運算元總是指向函式值增長最快的方向。這裡所說的是移動方向，而未提到移動量的大小。該量值稱為步長，記做$\alpha$。用向量來表示的話，梯度演算法的迭代公式如下：
$w:=w+\alpha\nabla_{w} f(w)$
該公式將一直被迭代執行，直至達到某個停止條件為止，比如迭代次數達到某個指定值或演算法達到某個可以允許的誤差範圍。

梯度下降演算法:
你最經常聽到的應該是梯度下降演算法，它與這裡的梯度上升演算法是一樣的，只是公式中的加法需要變成減法。因此，對應的公式可以寫成:
$w:=w-\alpha\nabla_{w} f(w)$
梯度上升演算法用來求函式的最大值，而梯度下降演算法用來求函式的最小值。

5.2.2 訓練演算法：使用梯度上升找到最佳引數

def load_data():
    data = pd.read_csv('./data/05testSet.txt', header=None, sep='\t', names=['x1', 'x2', 'class'])
    return data

def showData(myData):
    fig = plt.figure()
    ax = fig.add_subplot(111)
    plt.title('Scatter plot of training data')
    plt.xlabel('x1')
    plt.ylabel('x2')
    positive = myData[myData['class'].isin([1])]
    negative = myData[myData['class'].isin([0])]
    ax.scatter(positive['x1'], positive['x2'], s=30, c='purple', marker='+', label='Class 0')
    ax.scatter(negative['x1'], negative['x2'], s=30, c='b', marker='o', label='Class 1')
    plt.legend(loc='best')
    plt.show()

mydata = load_data()
showData(mydata)

圖中有100個樣本點，每個點包含兩個數值型特徵：X1和X2。在此資料集上，我們將通過使用梯度上升法找到最佳迴歸係數，也就是擬合出Logistic迴歸模型的最佳引數。

梯度上升法的虛擬碼如下：

每個迴歸係數初始化為1
重複R次：
   計算整個資料集的梯度
   使用alpha × gradient更新迴歸係數的向量
   返回迴歸係數

def loadDataSet():
    dataMat = []
    labelMat = []
    fr = open('./data/05testSet.txt')
    for line in fr.readlines():
        lineArr = line.strip().split()
        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
        labelMat.append(int(lineArr[2]))
    return dataMat, labelMat


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


def gradAscent(dataMat, classLabels, alpha, iteration):
    """
    梯度上升方法找最佳引數
    :param dataMat: 訓練集
    :param classLabels: 標籤
    :param alpha: 學習率
    :param iteration: 迭代次數
    :return:
    """
    dataMatrix = np.matrix(dataMat)
    labelMat = np.matrix(classLabels).transpose()
    m, n = dataMatrix.shape
    weights = np.ones((n, 1))
    for k in range(iteration):
        h = sigmoid(dataMatrix*weights)
        error = (labelMat - h)
        weights = weights + alpha * dataMatrix.T * error
    return weights

data, lables = loadDataSet()

alpha = 0.001
iteration = 500
my_weights = gradAscent(data, lables, alpha, iteration)
print(my_weights)

[[ 4.12414349]
 [ 0.48007329]
 [-0.6168482 ]]

5.2.3 分析資料：畫出決策邊界

上面已經解出了一組迴歸係數，它確定了不同類別資料之間的分隔線。那麼怎樣畫出該分隔線，從而使得優化的過程便於理解呢？

def show_boundary(data, weights):
    w = weights.getA()
    fig = plt.figure()
    ax = fig.add_subplot(111)
    plt.title('Decision Boundary')
    plt.xlabel('x1')
    plt.ylabel('x2')
    positive = data[data['class'].isin([1])]
    negative = data[data['class'].isin([0])]
    ax.scatter(positive['x1'], positive['x2'], s=30, c='purple', marker='+', label='Class 0')
    ax.scatter(negative['x1'], negative['x2'], s=30, c='b', marker='o', label='Class 1')
    x = np.arange(-3.0, 3.0, 0.1)
    y = (-w[0] - w[1]*x) / w[2]
    ax.plot(x, y, color='g')
    plt.legend(loc='best')
    plt.show()
    
show_boundary(mydata, my_weights)

5.2.4 訓練演算法：隨機梯度上升

梯度上升演算法在每次更新迴歸係數時都需要遍歷整個資料集，該方法在處理100個左右的資料集時尚可，但如果有數十億樣本和成千上萬的特徵，那麼該方法的計算複雜度就太高了。一種改進方法是一次僅用一個樣本點來更新迴歸係數，該方法稱為隨機梯度上升演算法。由於可以在新樣本到來時對分類器進行增量式更新，因而隨機梯度上升演算法是一個線上學習演算法。與“線上學習”相對應，一次處理所有資料被稱作是“批處理”。

隨機梯度上升演算法可以寫成如下的虛擬碼：

所有迴歸係數初始化為1
對資料集中每個樣本
   計算該樣本的梯度
   使用alpha × gradient更新迴歸係數值
返回迴歸係數值

def stocgradAscent(dataMat, classLabels, alpha):
    """
    隨機梯度上升
    :param dataMat: 訓練集
    :param classLabels: 標籤
    :param alpha: 學習率
    :return:
    """
    dataMat = np.array(dataMat)
    m, n = dataMat.shape
    weights = np.ones(n)
    for i in range(m):
        h = sigmoid(sum(dataMat[i] * weights))
        error = classLabels[i] - h
        weights = weights + alpha * error * dataMat[i]
    return weights

alpha = 0.01
my_weights = stocgradAscent(data, lables, alpha)
my_weights = np.matrix(my_weights.reshape([3, 1]))
print(my_weights)
show_boundary(mydata, my_weights)

[[ 1.01702007]
 [ 0.85914348]
 [-0.36579921]]

可以看到，隨機梯度上升演算法與梯度上升演算法在程式碼上很相似，但也有一些區別：第一，後者的變數h和誤差error都是向量，而前者則全是數值；第二，前者沒有矩陣的轉換過程，所有變數的資料型別都是NumPy陣列。

圖5-6展示了隨機梯度上升演算法在200次迭代過程中迴歸係數的變化情況。其中的係數2，也就是圖5-5中的X2只經過了50次迭代就達到了穩定值，但係數1和0則需要更多次的迭代。另外值得注意的是，在大的波動停止後，還有一些小的週期性波動。不難理解，產生這種現象的原因是存在一些不能正確分類的樣本點（資料集並非線性可分），在每次迭代時會引發係數的劇烈改變。我們期望演算法能避免來回波動，從而收斂到某個值。另外，收斂速度也需要加快。對於圖5-6存在的問題，可以通過修改程式的隨機梯度上升演算法來解決。

def stocgradAscent02(dataMat, classLabels, alpha, iteration):
    """
    改進後的隨機梯度上升
    :param dataMat: 訓練集
    :param classLabels: 標籤
    :param alpha: 學習率
    :param iteration: 迭代次數
    :return: 迴歸係數值
    """
    dataMat = np.array(dataMat)
    m, n = dataMat.shape
    weights = np.ones(n)
    for i in range(iteration):
        dataIndex = list(range(m))
        for j in range(m):
            # alpha每次迭代時需要調整
            alpha = 4 / (1.0 + i + j) + 0.01
            # 隨機選取樣本來更新迴歸係數
            randIndex = int(np.random.uniform(0, len(dataIndex)))
            h = sigmoid(sum(dataMat[randIndex] * weights))
            error = classLabels[randIndex] - h
            weights = weights + alpha * error * dataMat[randIndex]
            del(dataIndex[randIndex])
    return weights

alpha = 0.01
iteration = 150
my_weights = stocgradAscent02(data, lables, alpha, iteration)
my_weights = np.matrix(my_weights.reshape([3, 1]))
print(my_weights)
show_boundary(mydata, my_weights)

[[13.77472812]
 [ 0.89055145]
 [-2.17190635]]

比較圖5-7和圖5-6可以看到兩點不同。第一點是，圖5-7中的係數沒有像圖5-6裡那樣出現週期性的波動，這歸功於stocGradAscent02()裡的樣本隨機選擇機制；第二點是，圖5-7的曲線比圖5-6收斂的快很多，這是由於stocGradAscent02()可以收斂得更快。這次我們僅僅對資料集做了20次遍歷，而之前的方法是500次。

5.3 示例：從疝氣病症預測病馬的死亡率

本節將使用Logistic迴歸來預測患有疝病的馬的存活問題。這裡的資料包含368個樣本和28個特徵。我並非育馬專家，從一些文獻中瞭解到，疝病是描述馬胃腸痛的術語。然而，這種病不一定源自馬的胃腸問題，其他問題也可能引發馬疝病。該資料集中包含了醫院檢測馬疝病的一些指標，有的指標比較主觀，有的指標難以測量，例如馬的疼痛級別。

另外需要說明的是，除了部分指標主觀和難以測量外，該資料還存在一個問題，資料集中有30%的值是缺失的。下面將首先介紹如何處理資料集中的資料缺失問題，然後再利用Logistic迴歸和隨機梯度上升演算法來預測病馬的生死。

5.3.1 準備資料：處理資料中的缺失值

資料中的缺失值是個非常棘手的問題，有很多文獻都致力於解決這個問題。那麼，資料缺失究竟帶來了什麼問題？假設有100個樣本和20個特徵，這些資料都是機器收集回來的。若機器上的某個感測器損壞導致一個特徵無效時該怎麼辦？此時是否要扔掉整個資料？這種情況下，另外19個特徵怎麼辦？它們是否還可用？答案是肯定的。因為有時候資料相當昂貴，扔掉和重新獲取都是不可取的，所以必須採用一些方法來解決這個問題。

下面給出了一些可選的做法：

• 使用可用特徵的均值來填補缺失值；
• 使用特殊值來填補缺失值，如-1；
• 忽略有缺失值的樣本；
• 使用相似樣本的均值添補缺失值；
• 使用另外的機器學習演算法預測缺失值。

現在，我們對下一節要用的資料集進行預處理，使其可以順利地使用分類演算法。在預處理階段需要做兩件事：第一，所有的缺失值必須用一個實數值來替換，因為我們使用的NumPy資料型別不允許包含缺失值。這裡選擇實數0來替換所有缺失值，恰好能適用於Logistic迴歸。這樣做的直覺在於，我們需要的是一個在更新時不會影響係數的值。迴歸係數的更新公式如下：
$weights=weights+alpha*error*dataMatrix[randIndex]$
如果dataMatrix的某特徵對應值為0，那麼該特徵的係數將不做更新，即：
$weights=weights$
另外，由於sigmoid(0)=0.5，即它對結果的預測不具有任何傾向性，因此上述做法也不會對誤差項造成任何影響。基於上述原因，將缺失值用0代替既可以保留現有資料，也不需要對優化演算法進行修改。此外，該資料集中的特徵取值一般不為0，因此在某種意義上說它也滿足“特殊值”這個要求。

預處理中做的第二件