機器學習之Logistic回歸(五)

01-27

主要內容

Sigmoid函數和Logistoc回歸分類器
最優化理論初步
梯度下降最優化演算法
數據中的缺失項處理

我們將介紹最優化演算法，並利用他們訓練出一個非線性函數用於分類。

利用Logistic回歸的主要思想是：根據現有的數據對分類邊界線建立回歸公式，以此進行分類。

「回歸」一詞源於最佳擬合，表示要找到最佳擬合參數，使用的是最優化演算法。

Logistic回歸一般過程：

收集數據：採用任意方法收集數據
準備數據：由於需要進行距離計算，因此要求數據類型為數值型，另外，結構化數據為最佳。
分析數據：採用任意方法對於數據進行分析
訓練演算法：大部分時間將用於訓練，並將其轉換為對應的結構化數據
測試演算法：一旦訓練步驟完成，分類將會很快；
使用演算法：第一，我們需要輸入一些數據，並將其轉換為對應的結構化數值；第二，基於訓練好的回歸係數就可以對這些數值進行簡單的回歸計算，判斷他們屬於哪個類別，在這之後，我們就可以輸出的類別上做一些其他分析工作。

5.1 基於logistic回歸和sigmoid函數的分類

logistic回歸：

優點：計算代價不高，易於理解和實現；

缺點：容易欠擬合，分類精度可能不高。

適用數據類型：數值型和標稱型數據。

下面進入數學理論部分：

我們想要的是能接受所有的輸入然後預測出類別。該函數被稱之為海維塞德階躍遷函數或者單位躍遷函數。

有一個有類似於可以輸出0或者1這種性質的函數，就是Sigmoid函數。

Sigmoid計算公式如下：

下圖給出了Sigmoid函數在不同坐標尺下的曲線圖，當x為0時，Sigmoid函數值為0.5 ，隨著x增大，Sigmoid函數的值無限接近於1，隨著x的減小，Sigmoid函數的值無限接近於0 。實際上，Sigmoid函數看起來像一個階躍函數。

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

確定了分類器的函數形式之後，問題就變成了：最佳回歸係數為多少？如何確定大小？

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

Sigmoid函數的輸入記為z，公式如下：

其中x是分類器的輸入數據，向量w是我們要找的最佳參數。

為了尋找最佳參數，我們介紹一下梯度上升的最優化方法。我們將使用梯度上升的最優化方法來尋找最佳參數。

5.2.1 梯度上升法

梯度上升法基本思想：要找到某函數的最大值，最好的方法就是沿著該函數的梯度方向探尋。如果梯度為 ▽ ，則函數f(x,y)的梯度如下表示：

梯度要沿著x的方向移動 (?f(x+y))/?x

梯度要沿著y的方向移動 (?f(x+y))/?y

梯度上升演算法的迭代公式：

w: = w + α▽f(w)

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

圖5-3 有100個樣本點：每個點包含兩個數值型特徵：X1和X2.。我們將通過使用梯度上升法找到最佳回歸係數，也就是擬合出 Logistic回歸模型的最佳參數。

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

每個回歸係數初始化為1

重複R次：

計算整個數據集的梯度

使用alpha * gradient 更新回歸係數的向量

返回回歸係數

下面的代碼是梯度上升演算法的具體實現，創建名為 logRegres.py的文件，輸入代碼：

#Logistic回歸梯度上升最優化演算法nimport mathnfrom numpy import *nndef loadDataSet():n dataMat = []; labelMat = []n fr = open(rE:MLML_source_codemliaCh05testSet.txt)n for line in fr.readlines():n lineArr = line.strip().split()n dataMat.append([1.0,float(lineArr[0]),float(lineArr[1])])n labelMat.append(int(lineArr[2]))n return dataMat,labelMatnndef sigmoid(inX):n return 1.0/(1+exp(-inX))nndef gradAscent(dataMatIn,classLabels):n #調用numpy函數可以將數組轉行為矩陣數據類型n dataMatrix = mat(dataMatIn)n #將矩陣轉置n labelMat = mat(classLabels).transpose()n #看一下dataMatrix的大小，3*100的矩陣n m,n = shape(dataMatrix)n #目標移動步長n alpha = 0.01n #迭代次數n maxCycles =500n #n已經賦值為3n weights = ones((n,1))n for k in range(maxCycles):n h = sigmoid(dataMatrix*weights)n error = (labelMat - h)n weights = weights + alpha * dataMatrix.transpose() * errorn return weightsn

接下來測試代碼，載入這個函數

In [1]: import osnnIn [2]: os.chdir(E:MLML_IALogistic)nnIn [3]: import logRegresnnIn [4]: import logRegresnnIn [5]: reload(logRegres)nOut[5]: <module logRegres from logRegres.py>nnIn [6]: a,b= logRegres.loadDataSet()nnIn [40]: logRegres.gradAscent(a,b)nOut[40]:nmatrix([[ 4.12414349],n [ 0.48007329],n [-0.6168482 ]])n

5.3.2 分析數據：畫出決策邊界

上面解出的回歸係數，確定了不同類別數據之間的分割線。那麼如何畫出分割線呢？

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

我們最終需要解出X1和X2之間的關係式，運行代碼，

In [24]: import logRegresnnIn [25]: reload(logRegres)nOut[25]: <module logRegres from logRegres.py>nnIn [26]: a,b = logRegres.loadDataSet()nnIn [27]: logRegres.gradAscent(a,b)nOut[27]:nmatrix([[ 4.12414349],n [ 0.48007329],n [-0.6168482 ]])nnIn [28]: weights = logRegres.gradAscent(a,b)nnIn [29]: logRegres.plotBestFit(weights.getA())n

梯度上升演算法在500次迭代後得到的Logistic回歸最佳擬合曲線

最終圖上看上去是錯分了2個點。

儘管例子比較簡單，但是卻需要大量的計算。後面我們將對演算法作以改進。

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

梯度上升演算法在每次更新回歸係數時都需要去遍歷整個數據集，該方法在處理100個左右的數據集時速度可以，如果去處理數億或者幾千萬的特徵，這個方法的複雜度就太高了。

我們改進的方法就是一次僅用一個樣本點來更新回歸係數，這個方法叫做隨機梯度上升演算法。

由於可以在新樣本到來時對分類器進行增量式更新，因而隨機梯度上升演算法是一個在線學習演算法。與「在線學習」相對應，一次處理的所有數據被稱之為「批處理」。

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

所有的回歸係數初始都為1

對數據集中每個樣本

計算該樣本的梯度

使用alpha*gradient更新回歸係數值

返回歸係數數值

# 隨機梯度上升演算法ndef stocGradAscent0(dataMatrix,classLabels):n m,n = shape(dataMatrix)n alpha = 0.01n weights = ones(n)n for i in range(m):n h = sigmoid(sum(dataMatrix[i]*weights))n error = classLabels[i] - hn weights = weights + alpha * error * dataMatrix[i]n return weightsn

隨機梯度演算法與梯度上升演算法有一定的區別：

1.後者的h和誤差error都是向量，而前者則是數值；

2.前者沒有矩陣的轉換過程，所有變數的數據類型都是Numpy數組。

In [38]: import logRegresnnIn [39]: reload(logRegres)nOut[39]: <module logRegres from logRegres.pyc>nnIn [40]: a,b = logRegres.loadDataSet()nnIn [41]: weights = logRegres.stocGradAscent1(array(a),b)nnIn [42]: weightsnOut[42]: array([ 1.01702007, 0.85914348, -0.36579921])nnIn [43]: logRegres.plotBestFit(weights)n

隨機梯度上升演算法在上述數據集上執行結果，最佳擬合直線並非最佳分類線

本次運行後的曲線如果與上面的圖進行比較，其實不如上面的完美。那麼主要是因為上面的圖迭代了500次。

因為數據波動或者高頻波動，且收斂速度較慢，所以我們需要對演算法改進。

值得注意的是：主要做了三個方面的改進：

alpha在每次迭代的時候都會調整，這會緩解數據波動或者高頻波動。
通過隨機選取樣本來更新回歸係數，這樣可以減少周期性波動
增加了一個迭代參數作為第三個參數，演算法將按照給的新的參數值進行迭代

# 改進的隨機梯度演算法ndef stocGradAscent1(dataMatrix,classLabels,numIter=150):n m,n = shape(dataMatrix)n weights = ones(n)n"""n i 是樣本點的下標,j 是迭代次數n """n for j in range(numIter): n dataIndex = range(m)n for i in range(m):n alpha = 4/(1.0 + j + i) + 0.01 #alpha每次迭代時需要調整n randIndex = int(random.uniform(0,len(dataIndex)))n h = sigmoid(sum(dataMatrix[randIndex]*weights))n error = classLabels[randIndex] - hn weights = weights + alpha *error * dataMatrix[randIndex]n del(dataIndex[randIndex])n return weightsn

下面所用的計算量更少：

In [43]: reload(logRegres)nOut[43]: <module logRegres from logRegres.pyc>nnIn [44]: a,b = logRegres.loadDataSet()nnIn [45]: weights = logRegres.stocGradAscent1(array(a),b)nnIn [46]: logRegres.plotBestFit(weights)n

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

我們將使用logistics回歸來預測患有疝病的馬的存活問題。該數據有部分指標主觀和難以測量，而且存在30%的缺失值。

步驟：

1.收集數據：給定數據文件；

2.準備數據：用python解析文本文件並填充缺失值；

3.分析數據：可視化並觀察數據；

4.訓練演算法：使用最優化演算法，得到最佳的係數；

5.測試演算法：為了量化回歸的效果，我們需要觀察錯誤率。再去根據錯誤率判斷，我們是否需要回退到訓練階段，通過改變迭代的次數和步長等參數來得到更好的回歸係數；

6.使用演算法：實現一個簡單的命令行程序來手機馬的癥狀並輸出預測結果。

5.3.1 準備數據：處理數據中的缺失數據

數據的缺失值處理是一個非常棘手的問題，有時候複雜業務環境下數據相當昂貴，我們不能簡單的去丟棄掉，或者重新獲取。

方法：

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

我們現在要對要用的數據集做預處理。

1.所有的缺失值必須使用一個實數值來替換，因為我們使用numpy數據類型不允許包含缺失值，這裡選擇實數0來替換所有的缺失值，恰好使用會回歸。這樣做是因為需要一個在更新的時候不會影響係數的值。

回歸係數更新公式如下：

weights = weights + alpha *error * dataMatrix[randIndex]

2.如果數據集中類別標籤已經缺失，那麼就丟棄。

5.3.2 測試演算法：用logistic回歸進行分類

程序清單5-5 logistic 回歸分類函數

# 疝氣病預測兵馬死亡率：Logistic回歸分類器nndef classifyVector(inX,weights):n """ n 回歸係數：inX n 特徵向量：weights n """n prob = sigmoid(sum(inX*weights))n if prob > 0.5:n return 1.0n else:n return 0.0nndef colicTest():n frTrain = open(tem_dir + "" + "horseColicTraining.txt")n frTest = open(tem_dir + "" + "horseColicTest.txt")n trainingSet = [] ; trainingLabels = []nn for line in frTrain.readlines():n currLine = line.strip().split(t)n lineArr = []n for i in range(21):n lineArr.append(float(currLine[i]))n trainingSet.append(lineArr)n trainingLabels.append(float(currLine[21]))n #stocGradAscent1計算回歸係數n trainWeights = stocGradAscent1(array(trainingSet),trainingLabels,500)n errorCount = 0; numTestVec = 0.0nn for line in frTest.readlines():n numTestVec += 1.0n currLine = line.strip().split(t)n lineArr = []n for i in range(21):n lineArr.append(float(currLine[i]))n if int(classifyVector(array(lineArr),trainWeights)) != int(currLine[21]):n errorCount += 1n errorRate = (float(errorCount)/numTestVec)n print "the error rate of this test is: %f" %errorRaten return errorRatenndef multiTest():n numTests = 10;errorSum = 0.0n for k in range(numTests):n errorSum += colicTest()n print "after %d iterations the average error rate is:%f" n% (numTests,errorSum/float(numTests))n

測試代碼，結果如下：

In [46]: reload(logRegres)nOut[46]: <module logRegres from logRegres.py>nnIn [47]: logRegres.multiTest()nthe error rate of this test is: 0.313433nthe error rate of this test is: 0.313433nthe error rate of this test is: 0.313433nthe error rate of this test is: 0.283582nthe error rate of this test is: 0.417910nthe error rate of this test is: 0.328358nthe error rate of this test is: 0.462687nthe error rate of this test is: 0.313433nthe error rate of this test is: 0.358209nthe error rate of this test is: 0.388060nafter 10 iterations the average error rate is:0.349254n

10次迭代之後的結果是34.9% 。因為數據本身就有30%的缺失，得到的錯誤率結果並不差。

如果手動去調整 colicTest() 函數中的迭代次數和 stocGradAscent1() 中的步長，那麼最終的平均錯誤率就在20%左右。

5.4 小結

Logistic 回歸的目的是尋找一個非線性函數sigmoid的最佳擬合參數。求解過程可以由最優化演算法來完成。在最優化的演算法中，最常用的就是梯度上升演算法，梯度上升演算法又可以簡化為隨機梯度上升演算法。

隨機梯度上升演算法和梯度上升演算法的效果相當，但佔用更少的計算資源。此外，隨機梯度是一種在線演算法，可以在數據到來時就完成參數的更新，而不需要重新讀取整個數據集來進行批處理運算。

機器學習的一個重要的問題就是如何處理缺失數據，這個問題沒有標準答案，取決於實際的需求。