第 5 章 Logistic 回歸

[TOC]

本章內(nèi)容

• Sigmoid 函數(shù)和 Logistic 回歸分類器
• 最優(yōu)化理論初步
• 梯度下降最優(yōu)化算法
• 數(shù)據(jù)中的缺失項處理

假設(shè)現(xiàn)在有一些數(shù)據(jù)點胳徽，我們用一條直線對這些點進(jìn)行擬合（該線稱為最佳擬合直線）缤骨，這個擬合過程就稱為回歸亲桥。利用 Logistic 回歸進(jìn)行分類的主要思想是：根據(jù)現(xiàn)有數(shù)據(jù)對分類邊界建立回歸公式，以此進(jìn)行分類憎夷。

Logistic 回歸的一般過程：

1. 收集數(shù)據(jù)：采用任意方法收集數(shù)據(jù)。
2. 準(zhǔn)備數(shù)據(jù)：由于需要進(jìn)行距離計算，因此要求數(shù)據(jù)類型為數(shù)值型携御。另外诱咏，結(jié)構(gòu)化數(shù)據(jù)格式則最佳苔可。
3. 分析數(shù)據(jù)：采用任意方法對數(shù)據(jù)進(jìn)行分析。
4. 訓(xùn)練算法：大部分時間將用于訓(xùn)練袋狞，訓(xùn)練的目的是為了找到最佳的分類回歸系數(shù)硕蛹。
5. 測試算法：一旦訓(xùn)練步驟完成醇疼，分類將會很快。
6. 使用算法：
   1. 首先法焰，我們需要輸入一些數(shù)據(jù)秧荆，并將其轉(zhuǎn)換成對應(yīng)的結(jié)構(gòu)化數(shù)值；
   2. 接著埃仪，基于訓(xùn)練好的回歸系數(shù)就可以對這些數(shù)值進(jìn)行簡單的額回歸計算乙濒，判定它們屬于哪個類別；
   3. 在這之后卵蛉，我們就可以在輸出的類別上做一些其他分析工作颁股。

1. 基于 Logistic 回歸和 SIgmoid 函數(shù)的分類

Logistic 回歸：

• 優(yōu)點：計算代價不高，易于理解和實現(xiàn)
• 缺點：容易欠擬合傻丝，分類精度可能不高
• 適用數(shù)據(jù)類型：數(shù)值型和標(biāo)稱型數(shù)據(jù)甘有。

我們想要的函數(shù)應(yīng)該是，能接受所有的輸入然后預(yù)測出類別葡缰。例如亏掀，在兩個雷的情況下，上述函數(shù)輸出 0 或 1泛释，該函數(shù)稱為 海維塞德階躍函數(shù)（Heaviside step function）滤愕，或者直接稱為 ** 單位階躍函數(shù)**。

海維塞德階躍函數(shù)的問題在于：該函數(shù)在跳躍點上從 0 瞬間跳躍到 1怜校，這個瞬間跳躍過程有時很難處理间影。考慮另一個函數(shù)： Sigmoid 函數(shù)：

\sigma(z) = \frac{1}{1+e^{-z}}

圖5-1 給出 Sigmoid 函數(shù)在不同坐標(biāo)尺度下的兩條曲線圖茄茁。當(dāng) x 為 0 時魂贬， Sigmoid 函數(shù)值為 0.5.隨著 x 的增大，對應(yīng)的 Sigmoid 值將逼近于 1裙顽；而隨著 x 的減小随橘，Sigmoid 值將逼近于 0。

問題：最佳回歸系數(shù)是多少锦庸？如何確定它們的大谢帷？

圖5-1

2. 基于最優(yōu)化方法的最佳回歸系數(shù)確定

Sigmoid 函數(shù)的輸入記為 z：

z=w_0x_0+w_1x_1+...+w_nx_n

寫成向量形式：

z=W^TX

表示將這兩數(shù)值向量對應(yīng)元素相乘然后全部加起來即得到 z 值甘萧。其中的向量 X 是分類器的輸入數(shù)據(jù)萝嘁，向量 W 也就是我們要找到的最佳參數(shù)（系數(shù)）

2.1 梯度上升法

梯度上升法基于的思想是：要找到某函數(shù)的最大值，最好的方法是沿著該函數(shù)的梯度方向探尋扬卷。如果 梯度記為 $$牙言，則函數(shù) f(x,y) 的梯度：

\triangledown f(x,y)=\binom{\frac{\partial f(x,y)}{\partial x}}{\frac{\partial f(x,y)}{\partial y}}

這個梯度意味著要沿 x 的方向移動

\frac{\partial f(x,y)}{\partial x}

，沿 y 方向移動

\frac{\partial f(x,y)}{\partial y}

怪得。其中咱枉，函數(shù) f(x,y) 必須在待計算的點互換有定義并且可微卑硫。

一個具體的函數(shù)例子見圖5-2

圖5-2

圖5-2 中的梯度上升算法沿梯度方向移動了一步〔隙希可以看到欢伏，梯度算法總是指向函數(shù)值增長最快的方向這里所說的是移動方向，而未提到移動量的大小亿乳。該量值稱為步長硝拧，記作 $\alpha$。用向量來表示的話葛假，梯度算法的迭代公式：


該公式將一直被迭代執(zhí)行障陶，直至達(dá)到某個停止條件為止，比如迭代次數(shù)達(dá)到某個指定值或算法達(dá)到某個可以允許的誤差范圍聊训。

梯度下降算法：與這里的梯度上升算法是一樣的抱究，只是公式中的加號變成減號：

w:=w+\alpha\triangledown{f(w)}

梯度上升算法是用來求函數(shù)的最大值，梯度下降算法用來求函數(shù)的最小值

基于上面的內(nèi)容带斑，我們來看一個 Logistic 回歸分類器的應(yīng)用例子鼓寺，從圖5-3 可以看出我們采用的數(shù)據(jù)集：

圖5-3

2.2 訓(xùn)練算法：使用梯度上升找到最佳參數(shù)

圖5-3 中有100個樣本點，每個點包含兩個數(shù)值型特征：X1 和 X2遏暴。在此數(shù)據(jù)集上，我們將通過使用梯度上升法找到最佳回歸系數(shù)指黎，也就是擬合出 Logistic 回歸模型的最佳參數(shù)朋凉。

梯度上升法的偽代碼：

每個回歸系數(shù)初始化為1
重復(fù) R 次：
    計算整個數(shù)據(jù)集的梯度
    使用 alpha * gradient 更新回歸系數(shù)的向量
    返回回歸系數(shù)

添加 loadDataSet() 函數(shù)：

def loadDataSet():
    """
    讀取數(shù)據(jù)
    :return: 
    """
    dataMat = []
    labelMat = []
    fr = open('testSet.txt')
    for line in fr.readlines():
        lineArr = line.strip().split()
        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])]) # x0=1.0,x1,x2
        labelMat.append(int(lineArr[2]))
    return dataMat, labelMat

添加 sigmoid() 函數(shù)：

def sigmoid(inX):
    """
    sigmoid 函數(shù)
    :param inX:
    :return:
    """
    return 1.0/(1+np.exp(-inX))

添加 gradAscent() 函數(shù)：

def gradAscent(dataMatIn, classLabels):
    """
    梯度上升算法
    :param dataMatIn: 2維數(shù)組，每列分別代表每個不同的特征醋安，每行代表每個訓(xùn)練樣本
    :param classLabels: 類別標(biāo)簽
    :return:
    """
    dataMatrix = np.mat(dataMatIn)
    labelMat = np.mat(classLabels).transpose() # 行向量轉(zhuǎn)置成列向量
    m,n = np.shape(dataMatrix)
    alpha = 0.001 # 步長
    maxCycles = 500 # 迭代次數(shù)
    weights = np.ones((n,1))
    for k in range(maxCycles):
        h = sigmoid(dataMatrix * weights)
        error = (labelMat - h)
        weights = weights + alpha * dataMatrix.transpose() * error
    return weights

2.3 分析數(shù)據(jù)：畫出決策邊界

上面已經(jīng)解出一組回歸系數(shù)杂彭，它確定了不同類別數(shù)據(jù)之間的分隔線。那么怎樣畫出該分隔線吓揪，從而使得優(yōu)化的過程便于理解呢亲怠？

添加 plotBestFit() 函數(shù)：

def plotBestFit(wei):
    """
    畫出數(shù)據(jù)集和 Logistic 回歸最佳擬合直線的函數(shù)
    :param wei:
    :return:
    """
    import  matplotlib.pyplot as plt
    weights = wei.getA()
    # weights = wei
    dataMat, labelMat = loadDataSet()
    dataArr = np.array(dataMat)
    n = np.shape(dataArr)[0]
    xcord1 = []
    ycord1 = []
    xcord2 = []
    ycord2 = []
    for i in range(n):
        if int(labelMat[i]) == 1:
            xcord1.append(dataArr[i,1])
            ycord1.append(dataArr[i,2])
        else:
            xcord2.append(dataArr[i,1])
            ycord2.append(dataArr[i,2])
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(xcord1, ycord1, s=30, c='red', marker='s')
    ax.scatter(xcord2, ycord2, s=30, c='green')
    x = np.arange(-3.0, 3.0, 0.1)
    y = (-weights[0] - weights[1] * x)/weights[2]
    ax.plot(x,y)
    plt.xlabel('X1')
    plt.ylabel('X2')
    plt.show()

測試：

dataArr,labelMat = loadDataSet()
weights=gradAscent(dataArr,labelMat)
print(weights) #[[ 4.12414349] [ 0.48007329] [-0.6168482 ]]
plotBestFit(weights)

圖5-5

這個分類結(jié)果相當(dāng)不錯，從圖上看只錯分了兩到四個點柠辞。但是盡管例子簡單且數(shù)據(jù)量很小团秽，卻需要大量的計算（300次乘法）。下一節(jié)對算法進(jìn)行改進(jìn)叭首。

2.4 訓(xùn)練算法：隨機(jī)梯度上升

隨機(jī)梯度上升：一次僅用一個樣本點來更新回歸系數(shù)习勤。由于可以在新樣本到來時對分類器進(jìn)行增量式更新，因此也是一個在線學(xué)習(xí)算法焙格。

隨機(jī)梯度上升的偽代碼：

所有回歸系數(shù)初始化為1
對數(shù)據(jù)集中每個樣本
    計算該樣本的梯度
    使用 alpha*gradient 更新回歸系數(shù)值
返回回歸系數(shù)值

以為是隨機(jī)梯度上升算法的實現(xiàn)代碼：

def stocGradAscent0(dataMatrix, classLabels):
    """
    隨機(jī)梯度上升算法
    :param dataMatrix: 
    :param classLabels: 
    :return: 
    """
    m,n = np.shape(dataMatrix)
    alpha = 0.01
    weights = np.ones(n)
    for i in range(m):
        h = sigmoid(sum(dataMatrix[i]*weights))
        error = classLabels[i] - h
        weights = weights + alpha * error * dataMatrix[i]
    return weights

可以看到图毕，隨機(jī)梯度上升算法與梯度上升算法在代碼上很相似，有的區(qū)別是：第一眷唉，后者的變量 h 和 誤差 error 都是向量予颤，而前者則是數(shù)值囤官；第二，前者沒有矩陣的轉(zhuǎn)換過程蛤虐，所有變量的數(shù)據(jù)類型都是 numpy 數(shù)組党饮。

測試代碼：

dataArr,labelMat = loadDataSet()
weights2 = stocGradAscent0(np.array(dataArr),labelMat)
print('weights2:',weights2) #  [ 1.01702007  0.85914348 -0.36579921]
plotBestFit(weights2)

圖5-5

一個判斷優(yōu)化算法優(yōu)劣的可靠方法是看它是否收斂，也就是說參數(shù)是否達(dá)到了穩(wěn)定值笆焰，是否還會不斷地變化劫谅？。對此嚷掠，我們在 stocGradAscent0() 函數(shù)上做了些修改捏检，使其在整個數(shù)據(jù)集上運(yùn)行 200 次。最終繪制的三個回歸系數(shù)的變化情況如圖5-6.

圖5-6

圖5-6 展示了隨機(jī)梯度上升算法在 200 次迭代過程中回歸系數(shù)的變化情況不皆。其中的系數(shù)2贯城，也就是圖5-5 中的 X2 只經(jīng)過了 50 次迭代就達(dá)到穩(wěn)定值，但系數(shù) 1 和 0 則需要更多次的迭代霹娄。值得注意的是能犯，在大的波動停止后，還以一些小的周期性波動犬耻。產(chǎn)生這種現(xiàn)象的原因是存在一些不能正確分類的樣本點（數(shù)據(jù)集并非線性可分）踩晶，在每次迭代時會引發(fā)系數(shù)的劇烈改變

改進(jìn)的隨機(jī)算法 stocGradAscent1：

def stocGradAscent1(dataMatrix, classLabels, numIter=150):
    """
    改進(jìn)的隨機(jī)梯度上升算法
    :param dataMatrix:
    :param classLabels:
    :return:
    """
    m,n = np.shape(dataMatrix)
    weights = np.ones(n)
    for j in range(numIter):
        dataIndex = range(m)
        for i in range(m):
            alpha = 4/(1.0+j+i)+0.01 # alpha 每次迭代時需要調(diào)整
            randIndex = int(np.random.uniform(0,len(dataIndex))) # 隨機(jī)選取更新
            h = sigmoid(sum(dataMatrix[i]*weights))
            error = classLabels[i] - h
            weights = weights + alpha * error * dataMatrix[randIndex]
            del(dataIndex[randIndex])
    return weights

alpha 在每次迭代的時候都會調(diào)整，這回緩解圖5-6 上的數(shù)據(jù)波動或者高頻波動枕磁。另外雖然 alpha 會隨著迭代次數(shù)不斷減小渡蜻，但永遠(yuǎn)不會減小到 0 ，這是因為調(diào)整式子中還有一個常數(shù)項计济。必須這樣做的原因是為了保證在多次迭代之后新數(shù)據(jù)仍然具有一定的影響茸苇。如果要處理的問題是動態(tài)變化的，那么可以適當(dāng)加大常數(shù)項沦寂，來確保新的值獲得更大的回歸系數(shù)学密。還有，在降低 alpha 的函數(shù)中传藏，alpha 每次減少 1/(j+1)腻暮，其中 j 是迭代次數(shù)，i 是樣本點的下標(biāo)毯侦。這樣當(dāng) j<<max(i) 時西壮，alpha 就不是嚴(yán)格下降的，避免參數(shù)的嚴(yán)格下降也常見于模擬退火算法等其他優(yōu)化算法中叫惊。

圖5-7

比較圖5-7 和圖5-6 可以看出兩點不同： 1.圖5-7中的系數(shù)沒有像圖5-6那樣出現(xiàn)周期性的波動款青，這歸功于 stocGradAscent1 里的樣本隨機(jī)選擇機(jī)制；2.圖5-7的水平軸比圖5-6 短了很多霍狰，這是由于 stocGradAscent1 可以收斂得更快抡草。這次我們僅僅對數(shù)據(jù)集做了20次遍歷饰及，而之前是500次。

程序運(yùn)行后得到圖5-8的結(jié)果圖康震，該分隔線達(dá)到了與 gradAscent() 差不多的效果燎含，但是所使用的計算量更少。

圖5-8

3. 示例：從疝氣病癥預(yù)測病馬的死亡率

使用 Logistic 回歸來預(yù)測患有疝氣的馬的存活問題腿短。

示例：使用 Logistic 回歸估計馬疝病的死亡率

1. 收集數(shù)據(jù)：給定數(shù)據(jù)文件屏箍。
2. 準(zhǔn)備數(shù)據(jù)：用 python 解析文本文件并填充缺失值。
3. 分析數(shù)據(jù)：可視化并觀察數(shù)據(jù)橘忱。
4. 訓(xùn)練算法：使用優(yōu)化算法赴魁，找到最佳的系數(shù)。
5. 測試算法：為了量化回歸的效果钝诚，需要觀察錯誤率颖御。根據(jù)錯誤率決定是否回退到訓(xùn)練階段，通過改變迭代的次數(shù)和步長等參數(shù)得到更好的回歸系數(shù)凝颇。
6. 使用算法：實現(xiàn)一個簡單的命令行程序來收集馬的癥狀并輸出預(yù)測結(jié)果潘拱。

注意：除了部分指標(biāo)主觀和難以測量外，該數(shù)據(jù)還存在一個問題拧略，數(shù)據(jù)集中有 30% 的值是缺失的芦岂。

3.1 準(zhǔn)備數(shù)據(jù)：處理數(shù)據(jù)中的缺失值

數(shù)據(jù)缺失究竟帶來了什么問題？假設(shè)有100個樣本和20 個特征垫蛆，這些數(shù)據(jù)都是機(jī)器收集回來的禽最。若機(jī)器的某個傳感器損壞導(dǎo)致一個特征無效時該怎么辦？此時是否要扔掉整個數(shù)據(jù)月褥？這種情況下弛随，另外19個特征怎么辦瓢喉？它們是否還柯永松宁赤？答案是肯定的、因為有時候數(shù)據(jù)相當(dāng)昂貴栓票，扔掉和重新獲取都是不可取的决左，所以必須采取一些方法來解決這個問題：

• 使用可用特征的均值來填補(bǔ)缺失值。
• 使用特殊值來填補(bǔ)缺失值走贪，如 -1佛猛。
• 忽略缺失值的樣本。
• 使用相似樣本的均值填補(bǔ)缺失值坠狡。
• 使用另外的機(jī)器學(xué)習(xí)算法預(yù)測缺失值继找。

在預(yù)處理階段需要做兩件事：

1. 所有的缺失值必須用一個實數(shù)值來替換，因為我們使用的 Numpy 數(shù)據(jù)類型不允許包含缺失值逃沿。這里選擇實數(shù) 0 來替換所有缺失值婴渡，恰好能適用于 Logistic 回歸幻锁。
2. 如果在測試數(shù)據(jù)集中發(fā)現(xiàn)了一條數(shù)據(jù)的類別標(biāo)簽已經(jīng)缺失，那么我們的簡單做法是將該條數(shù)據(jù)丟棄边臼。這是因為類別標(biāo)簽與特征不同哄尔，很難確定采用某個合適的值來替換。

3.2 測試算法：用 Logistic 回歸進(jìn)行分類

使用 Logistic 回歸方法進(jìn)行分類并不需要做很多工作柠并，所需做的只是把測試集上每個特征向量乘以最優(yōu)化方法得來的回歸系數(shù)岭接，再將該乘積結(jié)果求和，最后輸入到 Sigmoid 函數(shù)中即可臼予。如果對應(yīng)的 Sigmoid 值大于 0.5 就預(yù)測類別標(biāo)簽為 1鸣戴，否則為 0.

添加 classifyVector 函數(shù)：

def classifyVector(inX, weights):
    """
    根據(jù) Sigmoid 值判斷分類
    :param inX:
    :param weights:
    :return:
    """
    prob = sigmoid(sum(inX*weights))
    if prob > 0.5:return 1.0
    else: return 0.0

添加 colicTest() 函數(shù)：

def colicTest():
    """
    打開訓(xùn)練集和測試集，并對數(shù)據(jù)進(jìn)行格式化處理瘟栖，以及訓(xùn)練葵擎、測試
    :return:
    """
    frTrain = open('horseColicTraining.txt')
    frTest = open('horseColicTest.txt')
    trainingSet = []
    trainingLabels = []
    for line in frTrain.readlines():
        currLine = line.strip().split('\t')
        lineArr = []
        for i in range(21):
            lineArr.append(float(currLine[i]))
        trainingSet.append(lineArr)
        trainingLabels.append(float(currLine[21]))
    trainWeights = stocGradAscent1(np.array(trainingSet), trainingLabels, 500)
    errorCount = 0
    numTestVec = 0.0
    for line in frTest.readlines():
        numTestVec += 1.0
        currLine = line.strip().split('\t')
        lineArr = []
        for i in range(21):
            lineArr.append(float(currLine[i]))
        if int(classifyVector(np.array(lineArr), trainWeights)) != int(currLine[21]):
            errorCount += 1
    errorRate = (float(errorCount)/numTestVec)
    print('the error rate of this test is : %f' % errorRate)
    return errorRate

添加 multiTest() 函數(shù)：

def multiTest():
    """
    重復(fù)測試，并求取平均值
    :return:
    """
    numTests = 10
    errorSum = 0.0
    for k in range(numTests):
        errorSum += colicTest()
    print('after %d iterations the average error rate is : %f' % (numTests, errorSum/float(numTests)))

測試：

multiTest()

如果調(diào)整 colicTest() 中的迭代次數(shù)和 stocGradAscent1() 中的步長半哟，平均錯誤率可以降到 20% 左右酬滤。

4. 本章小結(jié)

Logistic 回歸的目的是尋找一個非線性函數(shù) Sigmoid 的最佳擬合參數(shù)，求解過程可以由最優(yōu)化算法來完成寓涨。在最優(yōu)化算法中盯串，最常用的就是梯度上升算法，而梯度上升算法又可以簡化為隨機(jī)梯度上升算法戒良。

• 機(jī)器學(xué)習(xí)實戰(zhàn)（筆記）：第 1 章 機(jī)器學(xué)習(xí)基礎(chǔ)
• 機(jī)器學(xué)習(xí)實戰(zhàn)（筆記）：第 2 章 k-近鄰算法
• 機(jī)器學(xué)習(xí)實戰(zhàn)（筆記）：第 3 章 決策樹
• 機(jī)器學(xué)習(xí)實戰(zhàn)（筆記）：第 4 章 基于概率論的分類方法：樸素貝葉斯
• 機(jī)器學(xué)習(xí)實戰(zhàn)（筆記）：第 5 章 Logistic 回歸
• 機(jī)器學(xué)習(xí)實戰(zhàn)（筆記）：第 6 章 支持向量機(jī)（未完）
• 機(jī)器學(xué)習(xí)實戰(zhàn)（筆記）：第 7 章 利用 AdaBoost 元算法提高分類性能（未完）
• 機(jī)器學(xué)習(xí)實戰(zhàn)（筆記）：第 8 章 預(yù)測數(shù)值型數(shù)據(jù)：回歸