Logistic回歸

主要內(nèi)容

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

假設(shè)現(xiàn)在有一些數(shù)據(jù)點(diǎn)吁断，我們用一條直線對(duì)這些點(diǎn)進(jìn)行擬合（該線稱為最佳擬合直線），這個(gè)擬合過程就稱作回歸憾股。利用Logistic回歸進(jìn)行分類的主要思想是：根據(jù)現(xiàn)有數(shù)據(jù)對(duì)分類邊界線建立回歸公式岔留。訓(xùn)練分類器時(shí)的做法就是尋找最佳擬合參數(shù)夏哭，使用的是最優(yōu)化算法。

Logistic回歸的一般過程

1. 收集數(shù)據(jù)：采用任意方法收集數(shù)據(jù)献联。
2. 準(zhǔn)備數(shù)據(jù)：由于需要進(jìn)行距離計(jì)算方庭，因此要求數(shù)據(jù)類型為數(shù)值型。另外酱固，結(jié)構(gòu)化數(shù)據(jù)格式則最佳械念。
3. 分析數(shù)據(jù)：采用任意方法對(duì)數(shù)據(jù)進(jìn)行分析。
4. 訓(xùn)練算法： 大部分時(shí)間將用于訓(xùn)練运悲，訓(xùn)練的目的是為了找到最佳的分類回歸系數(shù)龄减。
5. 測(cè)試算法：一旦訓(xùn)練步驟完成，分類將會(huì)很快班眯。
6. 使用算法：首先希停，我們下需要一些輸入數(shù)據(jù)烁巫，并將其轉(zhuǎn)換成對(duì)應(yīng)的數(shù)據(jù)結(jié)構(gòu)化數(shù)值；接著宠能，給予訓(xùn)練好的回歸系數(shù)就可以對(duì)這些數(shù)值進(jìn)行簡(jiǎn)單的回歸計(jì)算亚隙，判定它們屬于哪個(gè)類別。

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

Logistic回歸
優(yōu)點(diǎn)：計(jì)算代價(jià)不高违崇，易于理解和實(shí)現(xiàn)阿弃。

缺點(diǎn)：容易欠擬合，分類精度可能不高羞延。

使用數(shù)據(jù)類型：數(shù)值型和標(biāo)稱型數(shù)據(jù)渣淳。

1.1 Sigmoid函數(shù)

在這里我們需要這樣一個(gè)函數(shù)，能接收所有輸入數(shù)據(jù)（特征向量）然后預(yù)測(cè)出類別伴箩。例如入愧，有兩個(gè)類別，那么這個(gè)函數(shù)的輸出可以使0或1嗤谚。
可能大家會(huì)想起單位階躍函數(shù)棺蛛，又稱海維賽德階躍函數(shù)（Heaviside step function）。然而巩步，階躍函數(shù)的問題在于：該函數(shù)在跳躍點(diǎn)上從0瞬間跳躍到1鞠值，這個(gè)瞬間跳躍過程有時(shí)很難處理。幸好渗钉，另一個(gè)函數(shù)也有類似的性質(zhì)彤恶，且數(shù)學(xué)上更好處理，這就是Sigmoid函數(shù)鳄橘。
sigmoid函數(shù)具體公式如下：
$ \sigma(z)= \frac{1}{1 + e^{-z}} $

Sigmoid函數(shù)是一種階躍函數(shù)(step function)声离。在數(shù)學(xué)中，如果實(shí)數(shù)域上的某個(gè)函數(shù)可以用半開區(qū)間上的指示函數(shù)的有限次線性組合來表示瘫怜，那么這個(gè)函數(shù)就是階躍函數(shù)术徊。而數(shù)學(xué)中指示函數(shù)（indicator function）是定義在某集合X上的函數(shù)，表示其中有哪些元素屬于某一子集A鲸湃。

圖一 給出了Sigmoid函數(shù)在不同坐標(biāo)尺度下的兩條曲線圖赠涮。當(dāng)x為0時(shí)，Sigmoid函數(shù)值為0.5暗挑。隨著x的增大笋除，對(duì)應(yīng)的Sigmoid值將逼近于1；而隨著x的減小炸裆，Sigmoid值將逼近于0垃它。如果橫坐標(biāo)刻度足夠大（圖一），Sigmoid函數(shù)將看起來很像一個(gè)階躍函數(shù)。

圖一

確定了分類器的函數(shù)形式之后损合，現(xiàn)在的問題變成了：最佳回歸系數(shù)是多少？如何確定它們的大心锓住嫁审？

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

Sigmoid函數(shù)的輸入記為z,由下面公式得出：
$$ z = w_0+w_1+w_2+···+w_n$$
如果使用向量的寫法，上述公式可以寫成 $ z = w^Tx $ 赖晶。它表示將兩個(gè)數(shù)值向量的對(duì)應(yīng)元素相乘然后全部加起來即得到z值律适。 其中的向量x是分類器的輸入數(shù)據(jù)，向量w也就是我們要找到的最佳參數(shù)(參數(shù))遏插，從而使得分類盡可能地精確捂贿。

2.1 梯度上升法

梯度上升法基于的思想是：要找到某函數(shù)的最大值，最好的方法沿著該函數(shù)的梯度方向探尋胳嘲。如果梯度記為$\nabla$ 厂僧，則函數(shù)f(x,y)的梯度由下式表示：$$\nabla f(x,y) = \dbinom{\frac{ \partial f(x,y)}{ \partial x} }{ \frac{\partial f(x,y)}{\partial y}} $$

這是機(jī)器學(xué)習(xí)中最容易造成混淆的一個(gè)地方，但在數(shù)學(xué)上并不難了牛，需要做的只是牢記這些符號(hào)的意義颜屠。這個(gè)梯度意味著要沿x的方向移動(dòng) $ \frac{\partial f(x,y) }{\partial x} $，沿y的方向移動(dòng) $ \frac{\partial f(x,y)}{\partial y} $鹰祸。其中甫窟，函數(shù)$f(x,y)$必須要在待計(jì)算的點(diǎn)上有定義并且可微。

圖二

梯度下降算法

“梯度下降算法，它與這里的梯度上升算法是一樣的糠悼，只是公式中的加法需要變成減法届榄。因此，對(duì)應(yīng)的公式可以寫成：$$ w:=w+\alpha\nabla_{w}f(w) $$

梯度上升算法用來求函數(shù)的最大值倔喂，而梯度下降算法用來求函數(shù)的最小值铝条。

基于上面的內(nèi)容，我們來看一個(gè)Logistic回歸分類器的應(yīng)用例子席噩，從圖三可以看到我們采用的數(shù)據(jù)集班缰。

圖三

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

中有100個(gè)樣本點(diǎn)班挖，每個(gè)點(diǎn)包含兩個(gè)數(shù)值型特征：X1和X2鲁捏。在此數(shù)據(jù)集上，我們將通過使用梯度上升法找到最佳回歸系數(shù)萧芙，也就是擬合出Logistic回歸模型的最佳參數(shù)给梅。

梯度上升法的偽代碼如下：

每個(gè)回歸系數(shù)初始化為1

重復(fù)R次：

計(jì)算整個(gè)數(shù)據(jù)集的梯度

使用alpha × gradient更新回歸系數(shù)的向量

返回回歸系數(shù)

下面是梯度上升算法的具體實(shí)現(xiàn)：

def loadDataSet():
    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])])
        labelMat.append(int(lineArr[2]))
    return dataMat,labelMat

dataMat,labelMat = loadDataSet()
print(dataMat)
print("----------------------------")
print(labelMat)

[[1.0, -0.017612, 14.053064], [1.0, -1.395634, 4.662541], [1.0, -0.752157, 6.53862], [1.0, -1.322371, 7.152853], [1.0, 0.423363, 11.054677], [1.0, 0.406704, 7.067335], [1.0, 0.667394, 12.741452], [1.0, -2.46015, 6.866805], [1.0, 0.569411, 9.548755], [1.0, -0.026632, 10.427743], [1.0, 0.850433, 6.920334], [1.0, 1.347183, 13.1755], [1.0, 1.176813, 3.16702], [1.0, -1.781871, 9.097953], [1.0, -0.566606, 5.749003], [1.0, 0.931635, 1.589505], [1.0, -0.024205, 6.151823], [1.0, -0.036453, 2.690988], [1.0, -0.196949, 0.444165], [1.0, 1.014459, 5.754399], [1.0, 1.985298, 3.230619], [1.0, -1.693453, -0.55754], [1.0, -0.576525, 11.778922], [1.0, -0.346811, -1.67873], [1.0, -2.124484, 2.672471], [1.0, 1.217916, 9.597015], [1.0, -0.733928, 9.098687], [1.0, -3.642001, -1.618087], [1.0, 0.315985, 3.523953], [1.0, 1.416614, 9.619232], [1.0, -0.386323, 3.989286], [1.0, 0.556921, 8.294984], [1.0, 1.224863, 11.58736], [1.0, -1.347803, -2.406051], [1.0, 1.196604, 4.951851], [1.0, 0.275221, 9.543647], [1.0, 0.470575, 9.332488], [1.0, -1.889567, 9.542662], [1.0, -1.527893, 12.150579], [1.0, -1.185247, 11.309318], [1.0, -0.445678, 3.297303], [1.0, 1.042222, 6.105155], [1.0, -0.618787, 10.320986], [1.0, 1.152083, 0.548467], [1.0, 0.828534, 2.676045], [1.0, -1.237728, 10.549033], [1.0, -0.683565, -2.166125], [1.0, 0.229456, 5.921938], [1.0, -0.959885, 11.555336], [1.0, 0.492911, 10.993324], [1.0, 0.184992, 8.721488], [1.0, -0.355715, 10.325976], [1.0, -0.397822, 8.058397], [1.0, 0.824839, 13.730343], [1.0, 1.507278, 5.027866], [1.0, 0.099671, 6.835839], [1.0, -0.344008, 10.717485], [1.0, 1.785928, 7.718645], [1.0, -0.918801, 11.560217], [1.0, -0.364009, 4.7473], [1.0, -0.841722, 4.119083], [1.0, 0.490426, 1.960539], [1.0, -0.007194, 9.075792], [1.0, 0.356107, 12.447863], [1.0, 0.342578, 12.281162], [1.0, -0.810823, -1.466018], [1.0, 2.530777, 6.476801], [1.0, 1.296683, 11.607559], [1.0, 0.475487, 12.040035], [1.0, -0.783277, 11.009725], [1.0, 0.074798, 11.02365], [1.0, -1.337472, 0.468339], [1.0, -0.102781, 13.763651], [1.0, -0.147324, 2.874846], [1.0, 0.518389, 9.887035], [1.0, 1.015399, 7.571882], [1.0, -1.658086, -0.027255], [1.0, 1.319944, 2.171228], [1.0, 2.056216, 5.019981], [1.0, -0.851633, 4.375691], [1.0, -1.510047, 6.061992], [1.0, -1.076637, -3.181888], [1.0, 1.821096, 10.28399], [1.0, 3.01015, 8.401766], [1.0, -1.099458, 1.688274], [1.0, -0.834872, -1.733869], [1.0, -0.846637, 3.849075], [1.0, 1.400102, 12.628781], [1.0, 1.752842, 5.468166], [1.0, 0.078557, 0.059736], [1.0, 0.089392, -0.7153], [1.0, 1.825662, 12.693808], [1.0, 0.197445, 9.744638], [1.0, 0.126117, 0.922311], [1.0, -0.679797, 1.22053], [1.0, 0.677983, 2.556666], [1.0, 0.761349, 10.693862], [1.0, -2.168791, 0.143632], [1.0, 1.38861, 9.341997], [1.0, 0.317029, 14.739025]]
----------------------------
[0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0]

函數(shù)loadDataSet()，它的主要功能是打開文本文件testSet.txt并逐行讀取双揪。每行前兩個(gè)值分別是X1和X2动羽，第三個(gè)值是數(shù)據(jù)對(duì)應(yīng)的類別標(biāo)簽。此外渔期，為了方便計(jì)算运吓，該函數(shù)還將X0的值設(shè)為1.0渴邦。

點(diǎn)擊 testSet.txt獲取文件。

sigmoid函數(shù)

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

gradAscent函數(shù)

def gradAscent(dataMatIn,classLabels):
    #?（以下兩行）轉(zhuǎn)換為numpy矩陣數(shù)據(jù)類型
    dataMatrix = np.mat(dataMatIn)
    labelMat = np.mat(classLabels).transpose()
    
    m,n = np.shape(dataMatrix)
    alpha = 0.001
    maxCycles = 500
    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

梯度上升算法的實(shí)際工作是在函數(shù)gradAscent()里完成的拘哨，該函數(shù)有兩個(gè)參數(shù)谋梭。第一個(gè)參數(shù)是dataMatIn，它是一個(gè)2維NumPy數(shù)組倦青，每列分別代表每個(gè)不同的特征瓮床，每行則代表每個(gè)訓(xùn)練樣本。我們現(xiàn)在采用的是100個(gè)樣本的簡(jiǎn)單數(shù)據(jù)集产镐，它包含了兩個(gè)特征X1和X2隘庄，再加上第0維特征X0，所以dataMath里存放的將是100×3的矩陣癣亚。在?處丑掺，我們獲得輸入數(shù)據(jù)并將它們轉(zhuǎn)換成NumPy矩陣。這是本書首次使用NumPy矩陣述雾，如果你對(duì)矩陣數(shù)學(xué)不太熟悉街州，那么一些運(yùn)算可能就會(huì)不易理解。比如绰咽，NumPy對(duì)2維數(shù)組和矩陣都提供一些操作支持菇肃，如果混淆了數(shù)據(jù)類型和對(duì)應(yīng)的操作地粪，執(zhí)行結(jié)果將與預(yù)期截然不同取募。第二個(gè)參數(shù)是類別標(biāo)簽，它是一個(gè)1×100的行向量蟆技。為了便于矩陣運(yùn)算玩敏，需要將該行向量轉(zhuǎn)換為列向量，做法是將原向量轉(zhuǎn)置质礼，再將它賦值給labelMat旺聚。接下來的代碼是得到矩陣大小，再設(shè)置一些梯度上升算法所需的參數(shù)眶蕉。

變量alpha是向目標(biāo)移動(dòng)的步長(zhǎng)砰粹，maxCycles是迭代次數(shù)。在for循環(huán)迭代完成后造挽，將返回訓(xùn)練好的回歸系數(shù)碱璃。需要強(qiáng)調(diào)的是，在?處的運(yùn)算是矩陣運(yùn)算饭入。變量h不是一個(gè)數(shù)而是一個(gè)列向量嵌器，列向量的元素個(gè)數(shù)等于樣本個(gè)數(shù)，這里是100谐丢。對(duì)應(yīng)地爽航，運(yùn)算dataMatrix * weights代表的不是一次乘積計(jì)算蚓让，事實(shí)上該運(yùn)算包含了300次的乘積。

?中公式的前兩行讥珍，是在計(jì)算真實(shí)類別與預(yù)測(cè)類別的差值历极，接下來就是按照該差值的方向調(diào)整回歸系數(shù)。

import numpy as np
weights = gradAscent(dataMat,labelMat)
print weights

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

調(diào)用gradAscent函數(shù)衷佃，可以得到調(diào)整后的回歸系數(shù)执解。

畫出數(shù)據(jù)集和logistic回歸最佳擬合直線的函數(shù)

def plotBestFit(dataMat1,labelMat1,weights):
    import matplotlib.pyplot as plt  
    dataArr = np.array(dataMat1)
    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)
    # ?  最佳擬合直線
    y1 = (-weights[0]-weights[1]*x)/weights[2] 
    y = y1.reshape(60,1)
    ax.plot(x,y)
    
    plt.xlabel('X1')
    plt.ylabel('X2')
    plt.show()

?處設(shè)置了sigmoid函數(shù)為0「傩铮回憶5.2節(jié)衰腌，0是兩個(gè)類別（類別1和類別0）的分界處。因此觅赊，我們?cè)O(shè)定 $ 0 = w0x0+ w1x1 + w2x2 $右蕊，然后解出X2和X1的關(guān)系式（即分隔線的方程，注意X0=1）吮螺。

%matplotlib inline

plotBestFit(np.array(dataMat),np.array(labelMat),weights)

500次迭代擬合直線

圖四 梯度上升算法在500次迭代后得到的Logistic回歸最佳擬合直線
從這個(gè)分類結(jié)果相當(dāng)不錯(cuò)饶囚，從圖上看只錯(cuò)分了兩到四個(gè)點(diǎn)。但是鸠补，盡管雷子簡(jiǎn)單且數(shù)據(jù)集很小萝风，這個(gè)方法卻需要大量的計(jì)算（300次乘法）。因此下一節(jié)將對(duì)該算法稍作改進(jìn)紫岩，從而使它可以用在真實(shí)數(shù)據(jù)集上规惰。

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

梯度上升算法在每次更新回歸系數(shù)時(shí)都需要遍歷整個(gè)數(shù)據(jù)集，該方法在處理100個(gè)左右的數(shù)據(jù)集時(shí)尚可泉蝌，但如果有數(shù)十億樣本和成千上萬的特征歇万，那么該方案的計(jì)算復(fù)雜度就太高了。一種改進(jìn)方法是一次僅用一個(gè)樣本點(diǎn)來更新回歸系數(shù)勋陪，該方法稱為隨機(jī)梯度上升算法贪磺。由于可以在新樣本到來時(shí)對(duì)分類器進(jìn)行增量式更新，因而隨機(jī)梯度上升算法是一個(gè)在線學(xué)習(xí)算法诅愚。與“在線學(xué)習(xí)”相對(duì)應(yīng)寒锚，一次處理所有數(shù)據(jù)被稱作是“批處理”。

隨機(jī)梯度上升算法可以寫成如下的偽代碼：

所有回歸系數(shù)初始化1

對(duì)數(shù)據(jù)集中每個(gè)樣本

   計(jì)算該樣本的梯度

   使用alpha x gradient 更新回歸系數(shù)值

返回回歸系數(shù)值

以下是隨機(jī)梯度上升算法實(shí)現(xiàn)代碼违孝。

def stocGradAscent0(dataMatrix, classLabels):
    m,n = shape(dataMatrix)
    alpha = 0.01
    weights = 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ù)組衅鹿。

下面驗(yàn)證該方法的結(jié)果:

from numpy import *
dataArr,labelMat=loadDataSet()
weights=stocGradAscent0(array(dataArr),labelMat)
plotBestFit(dataArr,labelMat,weights)

image

執(zhí)行完畢后將得到 圖五 所示的最佳擬合直線圖，該圖與圖四有一些相似之處过咬〈蟛常可以看到，擬合出來的直線效果還不錯(cuò)掸绞，但并不像圖四那樣完美泵三。這里的分類器錯(cuò)分了三分之一的樣本。

直接比較程序清單3和程序清單1的代碼結(jié)果是不公平的衔掸，后者的結(jié)果是在整個(gè)數(shù)據(jù)集上迭代了500次才得到的烫幕。一個(gè)判斷優(yōu)化算法優(yōu)劣的可靠方法是看它是否收斂，也就是說參數(shù)是否達(dá)到了穩(wěn)定值敞映，是否還會(huì)不斷地變化较曼？對(duì)此，我們?cè)诔绦蚯鍐?-3中隨機(jī)梯度上升算法上做了些修改振愿，使其在整個(gè)數(shù)據(jù)集上運(yùn)行200次捷犹。最終繪制的三個(gè)回歸系數(shù)的變化情況如圖六所示。

image

圖六展示了隨機(jī)梯度上升算法在200次迭代過程中回歸系數(shù)的變化情況。其中的系數(shù)2胳蛮，也就是圖5-5中的X2只經(jīng)過了50次迭代就達(dá)到了穩(wěn)定值销凑，但系數(shù)1和0則需要更多次的迭代丛晌。另外值得注意的是仅炊，在大的波動(dòng)停止后，還有一些小的周期性波動(dòng)澎蛛。不難理解抚垄，產(chǎn)生這種現(xiàn)象的原因是存在一些不能正確分類的樣本點(diǎn)（數(shù)據(jù)集并非線性可分），在每次迭代時(shí)會(huì)引發(fā)系數(shù)的劇烈改變谋逻。我們期望算法能避免來回波動(dòng)呆馁，從而收斂到某個(gè)值。另外毁兆，收斂速度也需要加快浙滤。

對(duì)于圖六存在的問題，可以通過修改程序stocGradAscent0 的隨機(jī)梯度上升算法來解決气堕，代碼如下：

def stocGradAscent1(dataMatrix, classLabels, numIter=150):
    import numpy as np
    m,n = np.shape(dataMatrix)
    weights = np.ones(n)
    for j in range(numIter):         
        dataIndex = range(m)
        for i in range(m):
            #?  alpha每次迭代時(shí)需要調(diào)整
            alpha = 4/(1.0+j+i)+0.01
            #?  隨機(jī)選取更新 
            randIndex = int(np.random.uniform(0,len(dataIndex)))
            h = sigmoid(sum(dataMatrix[randIndex]*weights))
            error = classLabels[randIndex] - h
            weights = weights + alpha * error * dataMatrix[randIndex]
            del(dataIndex[randIndex])
    return weights

stocGradAscent1 與 stocGradAscent0 類似纺腊，但增加了兩處代碼來進(jìn)行改進(jìn)畔咧。

第一處改進(jìn)在?處。一方面揖膜，alpha在每次迭代的時(shí)候都會(huì)調(diào)整誓沸，這會(huì)緩解圖六上的數(shù)據(jù)波動(dòng)或者高頻波動(dòng)。另外壹粟。雖然alpha會(huì)隨著迭代次數(shù)不斷減小拜隧，但永遠(yuǎn)不會(huì)減小到0，這是因?yàn)?中還存在一個(gè)常數(shù)項(xiàng)趁仙。必須這樣做的原因是為了保證多次迭代之后新數(shù)據(jù)仍然具有一定得影響洪添。如果要處理的問題是動(dòng)態(tài)變化的，那么可以適當(dāng)加大上述常數(shù)項(xiàng)雀费，來確保新的值獲得更大的回歸系數(shù)薇组。另一點(diǎn)值得注意的是，在降低alpha的函數(shù)中坐儿，alpha每次減少** 1/(j+i)** 律胀，其中j時(shí)迭代次數(shù)，i是樣本點(diǎn)的下標(biāo)貌矿。這樣當(dāng)j<<max(i)時(shí)炭菌，alpha就不是嚴(yán)格下降的。避免參數(shù)的嚴(yán)格下降也常見于模擬退火算法等其他優(yōu)化算法逛漫。

第二個(gè)改進(jìn)的地方在?處黑低，這里通過隨機(jī)選取樣本來更新回歸系數(shù)。這種方法將減少周期性的波動(dòng)（如圖五中的波動(dòng)）酌毡。這種方法每次隨機(jī)從列表中選出一個(gè)值克握，然后從列表中刪掉該值（再進(jìn)行下一次迭代）。

此外枷踏，改進(jìn)算法還增加了一個(gè)迭代次數(shù)作為第3個(gè)參數(shù)菩暗。如果該參數(shù)沒有給定的話，算法將默認(rèn)迭代150次旭蠕。如果給定停团，那么算法將按照新的參數(shù)值進(jìn)行迭代。

image

該方法比采用固定alpha的方法收斂速度更快佑稠。

比較圖七和圖六可以看到兩點(diǎn)不同。第一點(diǎn)是旗芬，圖七中的系數(shù)沒有像圖六里那樣出現(xiàn)周期性的波動(dòng)舌胶，這歸功于stocGradAscent1()里的樣本隨機(jī)選擇機(jī)制；第二點(diǎn)是疮丛，圖七的水平軸比圖六短了很多幔嫂，這是由于stocGradAscent1()可以收斂得更快漱办。這次我們僅僅對(duì)數(shù)據(jù)集做了20次遍歷，而之前的方法是500次婉烟。

看一下在同一個(gè)數(shù)據(jù)集上的分類效果：

%matplotlib inline
dataMat,labelMat = loadDataSet()
weights = stocGradAscent1(np.array(dataMat),np.array(labelMat))
plotBestFit(np.array(dataMat),np.array(labelMat),weights)

png

從結(jié)果圖中可以看出娩井，分割線達(dá)到了與GradientAscent()差不多的效果，但是所使用的計(jì)算量更少似袁。默認(rèn)迭代次數(shù)是150洞辣，可以通過stocGradAscent()的第3個(gè)參數(shù)來對(duì)此進(jìn)行修改。

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

下面將使用Logistic回歸來預(yù)測(cè)患有疝病的馬的存活問題昙衅。這里的數(shù)據(jù)1包含368個(gè)樣本和28個(gè)特征扬霜。我并非育馬專家，從一些文獻(xiàn)中了解到而涉，疝病是描述馬胃腸痛的術(shù)語著瓶。然而，這種病不一定源自馬的胃腸問題啼县，其他問題也可能引發(fā)馬疝病材原。該數(shù)據(jù)集中包含了醫(yī)院檢測(cè)馬疝病的一些指標(biāo)，有的指標(biāo)比較主觀季眷，有的指標(biāo)難以測(cè)量余蟹，例如馬的疼痛級(jí)別。

數(shù)據(jù)集來自2010年1月11日的UCI機(jī)器學(xué)習(xí)數(shù)據(jù)庫(http://archive.ics.uci.edu/ml/datasets/Horse+Colic)子刮。該數(shù)據(jù)最早由加拿大安大略省圭爾夫大學(xué)計(jì)算機(jī)系的Mary McLeish和Matt Cecile收集威酒。

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

1. 收集數(shù)據(jù)：給定數(shù)據(jù)文件。
2. 準(zhǔn)備數(shù)據(jù)：用Python解析文本文件并填充缺失值挺峡。
3. 分析數(shù)據(jù)：可視化并觀察數(shù)據(jù)葵孤。
4. 訓(xùn)練算法：使用優(yōu)化算法，找到最佳的系數(shù)橱赠。
5. 測(cè)試算法：為了量化回歸的效果尤仍，需要觀察錯(cuò)誤率。根據(jù)錯(cuò)誤率決定是否回退到訓(xùn)練階段病线，通過改變迭代的次數(shù)和步長(zhǎng)等參數(shù)來得到更好的回歸系數(shù)吓著。
6. 使用算法：實(shí)現(xiàn)一個(gè)簡(jiǎn)單的命令行程序來收集馬的癥狀并輸出預(yù)測(cè)結(jié)果。

另外需要說明的是送挑，除了部分指標(biāo)主觀和難以測(cè)量之外，該數(shù)據(jù)集還存在一個(gè)問題暖眼，數(shù)據(jù)集中有30%的數(shù)據(jù)值是缺失的惕耕。下面將首先介紹如何處理數(shù)據(jù)集中的數(shù)據(jù)缺失問題，然后再利用Logistic回歸和隨機(jī)梯度上升算法來預(yù)測(cè)病馬的生死诫肠。

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

數(shù)據(jù)中的缺失值是個(gè)非常棘手的問題司澎，有很多文獻(xiàn)都致力于解決這個(gè)問題欺缘。那么，數(shù)據(jù)缺失究竟帶來了什么問題挤安？假設(shè)有100個(gè)樣本和20個(gè)特征谚殊，這些數(shù)據(jù)都是機(jī)器收集回來的。若機(jī)器上的某個(gè)傳感器損壞導(dǎo)致一個(gè)特征無效時(shí)該怎么辦蛤铜？此時(shí)是否要扔掉整個(gè)數(shù)據(jù)嫩絮？這種情況下，另外19個(gè)特征怎么辦围肥？它們是否還可用剿干？答案是肯定的。因?yàn)橛袝r(shí)候數(shù)據(jù)相當(dāng)昂貴穆刻，扔掉和重新獲取都是不可取的置尔，所以必須采用一些方法來解決這個(gè)問題。

下面是一些可選的做法：

• 使用可用特征的均值來填補(bǔ)缺失值氢伟；
• 使用特殊值來填補(bǔ)缺失值榜轿，如-1；
• 忽略有缺失值的樣本朵锣；
• 使用相似樣本的均值填補(bǔ)缺失值差导；
• 使用另外的機(jī)器學(xué)習(xí)算法預(yù)測(cè)缺失值。

在數(shù)據(jù)預(yù)處理階段需要做兩件事：第一猪勇，所有的缺失值必須用一個(gè)實(shí)數(shù)來替換设褐，因?yàn)槲覀兪褂玫腘umPy數(shù)據(jù)類型不允許包含缺失值。這里選擇實(shí)數(shù)0來替換所有缺失值泣刹，恰好能適用于Logistic回歸助析。原因在于，我們需要的是一個(gè)在更新時(shí)不會(huì)影響系數(shù)的值椅您⊥饧剑回歸系數(shù)的更新公式如下：
$$ weights = weights + alpha * error * dataMatrix[randIndex] $$

如果dataMatrix的某特征對(duì)應(yīng)值為0，那么該特征的系數(shù)將不做更新掀泳，即：
$$ weights = weights $$

另外雪隧，由于sigmoid(0)=0.5，即它對(duì)結(jié)果的預(yù)測(cè)不具有任何傾向性员舵，因此上述做法也不會(huì)對(duì)誤差項(xiàng)造成任何影響脑沿。基于上述原因马僻，將缺失值用0代替既可以保留現(xiàn)有數(shù)據(jù)庄拇，也不需要對(duì)優(yōu)化算法進(jìn)行修改。此外，該數(shù)據(jù)集中的特征取值一般不為0措近，因此在某種意義上說它也滿足“特殊值”這個(gè)要求溶弟。

預(yù)處理中做的第二件事是，如果在測(cè)試數(shù)據(jù)集中發(fā)現(xiàn)了一條數(shù)據(jù)的類別標(biāo)簽已經(jīng)缺失瞭郑，那么我們的簡(jiǎn)單做法是將該條數(shù)據(jù)丟棄辜御。這是因?yàn)轭悇e標(biāo)簽與特征不同，很難確定采用某個(gè)合適的值來替換屈张。采用Logistic回歸進(jìn)行分類時(shí)這種做法是合理的擒权，而如果采用類似kNN的方法就可能不太可行。

原始的數(shù)據(jù)集經(jīng)過預(yù)處理之后保存成兩個(gè)文件：horseColicTest.txt和horseColicTraining.txt⊥嗉耄現(xiàn)在我們有一個(gè)“干凈”可用的數(shù)據(jù)集和一個(gè)不錯(cuò)的優(yōu)化算法菜拓，下面將把這些部分融合在一起訓(xùn)練出一個(gè)分類器，然后利用該分類器來預(yù)測(cè)病馬的生死問題笛厦。

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

前面介紹了優(yōu)化算法纳鼎，但目前為止還沒有在分類上做任何實(shí)際嘗試。使用Logistic回歸方法進(jìn)行分類并不需要做很多工作裳凸，所需做的只是把測(cè)試集上每個(gè)特征向量乘以最優(yōu)化方法得來的回歸系數(shù)贱鄙，再將該乘積結(jié)果求和，最后輸入到Sigmoid函數(shù)中即可姨谷。如果對(duì)應(yīng)的Sigmoid值大于0.5就預(yù)測(cè)類別標(biāo)簽為1逗宁，否則為0。

Logistic回歸分類函數(shù)：

def classifyVector(inX,weights):
    prob = sigmoid(sum(inX*weights))
    if prob > 0.5:
        return 1.0
    else:
        return 0.0

函數(shù)classifyVector()梦湘，它以回歸系數(shù)和特征向量作為輸入來計(jì)算對(duì)應(yīng)的Sigmoid值瞎颗。如果Sigmoid值大于0.5函數(shù)返回1，否則返回0捌议。

def colicTest():
    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),np.array(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 "錯(cuò)誤率 = %f" %errorRate
    return errorRate

函數(shù)colicTest()哼拔，是用于打開測(cè)試集和訓(xùn)練集，并對(duì)數(shù)據(jù)進(jìn)行格式化處理的函數(shù)瓣颅。該函數(shù)首先導(dǎo)入訓(xùn)練集倦逐，同前面一樣，數(shù)據(jù)的最后一列仍然是類別標(biāo)簽宫补。數(shù)據(jù)最初有三個(gè)類別標(biāo)簽檬姥，分別代表馬的三種情況：“仍存活”、“已經(jīng)死亡”和“已經(jīng)安樂死”粉怕。這里為了方便健民，將“已經(jīng)死亡”和“已經(jīng)安樂死”合并成“未能存活”這個(gè)標(biāo)簽 。數(shù)據(jù)導(dǎo)入之后斋荞，便可以使用函數(shù)stocGradAscent1()來計(jì)算回歸系數(shù)向量荞雏。這里可以自由設(shè)定迭代的次數(shù)，例如在訓(xùn)練集上使用500次迭代平酿，實(shí)驗(yàn)結(jié)果表明這比默認(rèn)迭代150次的效果更好凤优。在系數(shù)計(jì)算完成之后，導(dǎo)入測(cè)試集并計(jì)算分類錯(cuò)誤率蜈彼。整體看來筑辨，colicTest()具有完全獨(dú)立的功能，多次運(yùn)行得到的結(jié)果可能稍有不同幸逆，這是因?yàn)槠渲杏须S機(jī)的成分在里面棍辕。如果在stocGradAscent1()函數(shù)中回歸系數(shù)已經(jīng)完全收斂，那么結(jié)果才將是確定的还绘。

def multiTest():
    numTests = 10;errorSum = 0.0
    for k in range(numTests):
        errorSum += colicTest()
    print " %d次平均錯(cuò)誤率是:%f" % (numTests,errorSum/float(numTests))

最后一個(gè)函數(shù)是multiTest()楚昭，其功能是調(diào)用函數(shù)colicTest()10次并求結(jié)果的平均值。調(diào)用multiTest()函數(shù)拍顷，可以看到打印的10次錯(cuò)誤率：

multiTest()

/usr/local/lib/python2.7/site-packages/ipykernel_launcher.py:2: RuntimeWarning: overflow encountered in exp
  


錯(cuò)誤率 = 0.358209
錯(cuò)誤率 = 0.313433
錯(cuò)誤率 = 0.462687
錯(cuò)誤率 = 0.402985
錯(cuò)誤率 = 0.358209
錯(cuò)誤率 = 0.402985
錯(cuò)誤率 = 0.313433
錯(cuò)誤率 = 0.373134
錯(cuò)誤率 = 0.373134
錯(cuò)誤率 = 0.402985
 10次平均錯(cuò)誤率是:0.376119

從上面的結(jié)果可以看到抚太，10次迭代之后的平均錯(cuò)誤率為37%。事實(shí)上昔案，這個(gè)結(jié)果并不差尿贫，因?yàn)閿?shù)據(jù)集有30%的數(shù)據(jù)已經(jīng)缺失。當(dāng)然踏揣，如果調(diào)整colicTest()中的迭代次數(shù)和stochGradAscent1()中的步長(zhǎng)庆亡，平均錯(cuò)誤率可以降到20%左右。

小結(jié)

Logistic回歸的目的是尋找一個(gè)非線性函數(shù)Sigmoid的最佳擬合參數(shù)捞稿，求解過程可以由最優(yōu)化算法來完成又谋。在最優(yōu)化算法中，最常用的就是梯度上升算法娱局，而梯度上升算法又可以簡(jiǎn)化為隨機(jī)梯度上升算法彰亥。
隨機(jī)梯度上升算法與梯度上升算法的效果相當(dāng)，但占用更少的計(jì)算資源铃辖。此外剩愧，隨機(jī)梯度上升是一個(gè)在線算法，它可以在新數(shù)據(jù)到來時(shí)就完成參數(shù)更新娇斩，而不需要重新讀取整個(gè)數(shù)據(jù)集來進(jìn)行批處理運(yùn)算仁卷。
機(jī)器學(xué)習(xí)的一個(gè)重要問題就是如何處理缺失數(shù)據(jù)。這個(gè)問題沒有標(biāo)準(zhǔn)答案犬第，取決于實(shí)際應(yīng)用中的需求〗趸現(xiàn)有一些解決方案，每種方案都各有優(yōu)缺點(diǎn)歉嗓。