目錄
5.1 基于logistic回歸和sigmoid函數的分類
5.2基于最優化方法的最佳回歸系數确定
5.2.1 梯度上升法
5.2.2 訓練算法：使用梯度上升法找到最佳參數
5.2.3 分析資料：畫出決策邊界
5.2.4 訓練算法：随機梯度上升
5.3 示例：從疝氣病症預測病馬的死亡率
5.3.1 準備資料：處理資料中的缺失值
5.3.2 測試算法：用logistic回歸進行分類
5.4 本章小結
假設現在有一些資料點，我們用一條直線對這些點進行拟合（改線稱為最佳拟合直線），這個拟合過程就叫作回歸。利用logistic回歸進行分類的主要思想是：根據現有資料對分類邊界線建立回歸公式，以此進行分類。這裡的“回歸”一詞源于最佳拟合，表示要找到最佳拟合參數集。訓練分類器時的做法就是尋找最佳拟合參數，使用的是最優化算法。
logistic回歸的一般過程
（1）收集資料：采用任意方法收集資料。
（2）準備資料：由于需要進行距離計算，是以要求資料類型為數值型。另外，結構化資料格式則最佳。
（3）分析資料：采用任意方式對資料進行分析。
（4）訓練算法：大部分時間将用于訓練，訓練的目的是為了找到最佳的分類回歸系數。
（5）測試算法：一旦訓練步驟完成，分類将會很快。
（6）使用算法：首先，需要輸入一些資料，并将其轉換成對應的結構化數值；接着，基于訓練好的回歸系數就可以對這些數值進行簡單的回歸計算，判定它們屬于哪個類别；在這之後，就可以在輸出的類别上做一些其他分析工作。
5.1 基于logistic回歸和sigmoid函數的分類
logistic回歸
優點：計算代價不高，易于了解和實作。
缺點：容易欠拟合，分類精度可能不高。
适用資料類型：數值型和标稱型資料。
我們想要的函數應該是，能接受所有的輸入然後預測出類别。例如，在兩個類的情況下，上述函數輸出0或1，該函數稱為‘海維賽德階躍函數（heaviside step function）’或直接稱為機關階躍函數。然而，海維賽德階躍函數的問題在于：該函數在階躍點上從0瞬間跳躍到1，這個瞬間跳躍過程有時很難處理。幸好，另一個函數也有類似的性質，且數學上更易處理，這就是sigmoid函數。sigmoid函數具體的計算公式如下：
下圖給出了sigmoid函數在不同坐标尺度下的兩條曲線圖。當x為0時，sigmoid函數值為0.5.随着x的增大，對應的sigmoid值将逼近于1；而随着x的減小，sigmoid值将逼近于0.如果橫坐标刻度足夠大，sigmoid函數看起來很像一個階躍函數：
兩種坐标尺度下的sigmoid函數圖。上圖的橫坐标是-6到6，這時的曲線變化較為平滑；下圖橫坐标尺度足夠大，可以看到，在x=0點處sigmoid函數看起來很像階躍函數。
是以，為了實作logistic回歸分類器，可以在每個特征上都乘于一個回歸系數，然後把所有的結果值相加，将這個總和帶入sigmoid函數中，進而得到一個範圍在0~1之間的數值。任何大于0.5的資料被分為1類，小于0.5即被歸入0類。是以，logistic回歸也可以被看成是一種機率估計。
确定了分類器的函數形式之後，現在的問題變成了：最佳回歸系數是多少？如何确定它們的大小？
5.2基于最優化方法的最佳回歸系數确定
sigmoid函數的輸入記為z，由下面公式得出：
如果采用向量的寫法，上述公式可以寫成
，它表示将這兩個數值向量對應元素相乘然後全部加起來即得到z值。其中的向量x是分類器的輸入資料，向量w就是要找的最佳參數（系數），進而使得分類器盡可能地精确。為了尋找該最佳參數，需要用到最優化理論的一些知識。
5.2.1 梯度上升法
梯度上升法基于的思想是：要找到某函數的最大值，最好的方法就是沿着該函數的梯度方向探尋。如果梯度記為
，則函數f（x,y）的梯度由下式表示：
這是機器學習中最易造成混淆的一個地方，但在數學上并不難，需要做的隻是牢記這些符号的意義。這個梯度意味着要沿x的方向移動
,沿y的方向移動
。其中，函數f(x,y)必須要在待計算的點上有定義并且可微。
梯度上升算法到達每個點後都會重新估計移動的方向。從p0開始，計算完該點的梯度，函數就根據梯度移動到下一點p1.在p1點，梯度再次被重新計算，并沿新的梯度方向移動到p2.如此循環疊代，知道滿足停止條件。疊代的過程中，梯度算子總是保證我們能選取到最佳的移動方向。
從上圖中可以看到，梯度算子總是指向函數值增長最快的方向。這裡所說的是移動方向，而未提到移動量的大小。該量值稱為步長，記作
。用向量來表示的話，梯度算法的疊代公式如下：
該公式将一直被疊代執行，直至達到某個停止條件為止，比如疊代次數達到某個指定值或算法達到某個可以允許的誤差範圍。
梯度下降法：梯度下降法和梯度上升法是一樣的。梯度上升法用來求函數的最大值，而梯度下降法用來求函數的最小值。
5.2.2 訓練算法：使用梯度上升法找到最佳參數
上圖中有100個樣本點，每個點包含兩個數值型特征：x1和x2.在此資料集上，我們将通過使用梯度上升法找到最佳回歸系數，也就是拟合出logistic回歸模型的最佳參數。
梯度上升法的僞代碼：
每個回歸系數初始化為1
重複R次：
。。計算整個資料集的梯度
。。使用alphaxgradient更新回歸系數的向量
。。傳回回歸系數
下面用代碼實作梯度上升法：
首先做一個包含100個樣本點的資料集放到testSet.txt檔案中（其中第一個數是x1，第二個數是x2，第三個數代表類别）：
建立名為logRegres.py的檔案，輸入下列代碼：
from numpy import *
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
def sigmoid(inX):
return 1.0/(1+exp(-inX))
def gradAscent(dataMatIn,classLabels):
dataMatrix=mat(dataMatIn)
labelMat=mat(classLabels).transpose()
m,n=shape(dataMatrix)
alpha=0.001
maxCycles=500
weights=ones((n,1))
for k in range(maxCycles):
h=sigmoid(dataMatrix*weights)
error=labelMat-h
weights=weights+alpha*dataMatrix.transpose()*error
return weights
程式開頭提供了一個便利函數loadDataSet()，它的主要功能是打開文本檔案testSet.txt并逐行讀取。每行前兩個值分别是x1和x2，第三個值是資料對應的類别标簽。此外，為了友善計算，該函數還将x0值設為1.0。
梯度上升法的實際工作是在函數gradAscent()裡完成的，該函數有兩個參數。第一個參數是dataMatIn，它是一個2維的numpy數組，每列分别代表每個不同的特征，每行則代表每個訓練樣本。現采用100個樣本的簡單資料集，它包含了兩個特征x1和x2，再加上第0維特征x0，是以dataMatIn裡存放的将是100x3的矩陣。mat（）将獲得的輸入資料轉換成numpy矩陣。第二個參數是類别标簽，它是一個1x100的行向量。為了便于矩陣運算，需要将該行向量轉換成列向量，做法是将原向量轉置，再将它指派給labelMat。
變量alpha是向目标移動的步長，maxCycle是疊代次數。在for循環疊代完成後，将傳回訓練好的回歸系數。需要強調的是，h=sigmoid(dataMatrix*weights)中的運算是矩陣運算。變量h不是一個數而是一個列向量，列向量的元素個數等于樣本個數，這裡是100.對應地，運算dataMatrix*weights代表的不止一次乘積運算，事實上該運算包含了300次的乘積。
運作代碼：
if __name__=="__main__":
dataArr,labelMat=loadDataSet()
print(gradAscent(dataArr,labelMat))
運作結果：
[[ 3.8218897 ]
[ 0.06112853]
[-0.5704236 ]]
5.2.3 分析資料：畫出決策邊界
上面解出的一組回歸系數，确定了不同類别資料之間的分隔線，現畫出分隔線，打開logRegres.py添加以下代碼：
#畫出資料集和logistic回歸最佳拟合直線的函數
def plotBestFit(wei):
import matplotlib.pyplot as plt
weights=wei.getA() #将矩陣轉換為數組
dataMat,labelMat=loadDataSet()
dataArr=array(dataMat)
n=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,color='red',marker='s')
ax.scatter(xcord2, ycord2, s=30, color='green')
x=arange(-3,3,0.1)
y=(-weights[0]-weights[1]*x)/weights[2]
ax.plot(x,y)
plt.xlabel('x1')
plt.ylabel('x2')
plt.show()
運作程式：
if __name__=="__main__":
dataArr,labelMat=loadDataSet()
weights=gradAscent(dataArr,labelMat)
plotBestFit(weights)
運作結果：
可以看出雖然分類結果很好，但是盡管資料集很小，這個方法卻需要大量的計算（300次乘法）。
5.2.4 訓練算法：随機梯度上升
梯度上升法是每次更新回歸系數時都需要周遊整個資料集，該方法在處理100個左右的資料及時尚可，但如果有數十億樣本和成千上萬的特征，那麼該方法的計算複雜度就太高了。一種改進方法是一次僅用一個樣本點來更新回歸系數，該方法稱為随機梯度上升算法。由于可以在新樣本到來時對分類器進行增量式更新，因而随機梯度上升法是一個線上學習算法。與“線上學習”相對應，一次處理所有資料被稱作是“批處理”。
随機梯度上升法可以寫成以下的僞代碼：
所有回歸系數初始化為1
對資料集中每個樣本
。。計算該樣本的梯度
。。使用alphaxgradient更新回歸系數值
傳回回歸系數值
以下是随機梯度上升法的實作代碼：
#随機梯度上升算法
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
可以看到，随機梯度上升法與梯度上升法在代碼上很相似，但也有一些差別：第一，後者的變量h和誤差error都是向量，而前者則全是數值；第二，前者沒有矩陣的轉換過程，所有變量的資料類型都是numpy數組。
為了驗證該方法的結果，将上述代碼添加到logRegres.py中，運作程式：
if __name__=="__main__":
dataArr,labelMat=loadDataSet()
weights=stocGradAscent0(array(dataArr),labelMat)
plotBestFit(weights)
運作結果：
可以看到随機梯度上升法拟合出來的結果不夠完美，但是直接比較兩個結果是不公平的，梯度上升法的結果是在整個資料集上疊代500次才得到的。一個判斷優化算法優劣的可靠方法是看它是否收斂，也就是參數是否達到了穩定值，是否還會不斷地變化？對此，在程式中對随機梯度上升法做一些修改，使其在整個資料集上運作150次。
#改進的随機梯度上升算法
def stocGradAscent1(dataMatrix,classLabels,numIter=150):
import random
m,n=shape(dataMatrix)
weights=ones(n)
for j in range(numIter):
dataIndex=range(m)
for i in range(m):
alpha=4/(1+j+i)+0.01
randIndex=int(random.uniform(0,len(dataIndex)))
h=sigmoid(sum(dataMatrix[randIndex]*weights))
error=classLabels[randIndex]-h
weights=weights+alpha*error*dataMatrix[randIndex]
return weights
在改進的算法中，‘alpha=4/(1+j+i)+0.01’使得alpha在每次疊代的時候都會調整，這回緩解回歸系數波動或者高頻波動。另外，雖然alpha會随着疊代次數不斷減小，但永遠不會減小到0，因為式子中還存在一個常數項。必須這樣做的原因是為了保證在多次疊代之後新資料仍然具有一定的影響。如果要處理的問題是動态變化的，那麼可以适當加大上述常數項，來確定新的值獲更大的回歸系數。另一點值得注意的是，在降低alpha的函數中，alpha每次減少1/(j+i),其中j是疊代次數，i是樣本點的下标。這樣當j<<max(i)時，alpha就不是嚴格下降的。避免參數的嚴格下降也常見于模拟退火算法等其它優化算法中。
在改進的程式中，通過随機選取樣本來更新回歸系數。這種方法将減少周期性的波動。這種方法每次随機從清單中選出一個值，然後從清單中删掉該值（再進行下一次疊代）。
運作結果：
可以看到随機梯度上升法疊代150次達到了和梯度上升算法相似的效果。
5.3 示例：從疝氣病症預測病馬的死亡率
資料集來源：https://archive.ics.uci.edu/ml/datasets/Horse+Colic：包含368個樣本和28個特征，該資料集中包含了醫院檢測馬疝病的一些名額，有的名額比較主觀，有的名額難以測量，例如馬的疼痛級别。
示例：使用logistic回歸估計馬疝病的死亡率
（1）收集資料：給定資料檔案。
（2）準備資料：用python解析文本檔案并填充缺失值。
（3）分析資料：可視化并觀察資料。
（4）訓練算法：使用優化算法，找到最佳的系數。
（5）測試算法：為了量化回歸的效果，需要觀察錯誤率。根據錯誤率決定是否回退到訓練階段，通過改變疊代的次數和步長等參數來得到更好的回歸系數。
（6）使用算法：實作一個簡單的指令行程式來收集馬的症狀并輸出預測結果。
除了部分名額主觀和難以測量外，該資料還存在一個問題，資料集中有30%的值是缺失的。下面将首先介紹如何處理資料集中的資料缺失問題，然後再利用logistic回歸和随機梯度上升算法來預測病馬的生死。
5.3.1 準備資料：處理資料中的缺失值
資料中的缺失值是個非常棘手的問題，有很多文獻都緻力于解決這個問題。那麼，資料缺失究竟帶來了什麼問題？假設有100個樣本和20個特征，這些資料都是機器學習收集回來的。若機器上的某個傳感器損壞導緻一個特征無效時該怎麼辦？此時是否要扔掉整個資料？這種情況下，另外19個特征怎麼辦？它們是否還可用？答案是肯定的。因為有時候資料相當昂貴，扔掉和重新擷取都是不可取的，是以必須采用一些方法來解決這個問題。
下面是一些可選的做法：
*使用可用特征的均值來填補缺失值；
*使用特殊值來填補缺失值，如-1；
*忽略有缺失值的樣本；
*使用相似樣本的均值填補缺失值；
*使用另外的機器學習算法預測缺失值。
首先，需要對資料集進行預處理，使其可以順利地使用分類算法。在預處理階段需要做兩件事：第一，所有缺失的值必須用一個實數值來替換，因為使用的numpy資料類型不允許包含缺失值。這裡選擇實數0來替換所有缺失值，恰好能适用于logistic回歸。這樣做不會在更新時不會影響系數的值。回歸系數的更新公式如下：
如果dataMatrix的某特征對應值為0，那麼該特征的系數将不做更新，即：
另外，由于sigmoid（0）=0.5，即它對結果的預測不具有任何傾向性，是以上述做法也不會對誤差項造成任何影響。基于上述原因，缺失值用0代替可以保留現有資料，也不需要對優化算法進行修改。此外，資料集中的特征取值一般不為0，是以在某種意義上說它也滿足“特殊值”這個要求。
預進行中做的第二件事是，如果在測試資料集中發現了一條資料的類别标簽已經缺失，那麼簡單做法是将該條資料丢棄。這是因為類别标簽與特征不同，很難确定采用某個合适的值來替換。采用logistic回歸進行分類時這種做法是合理的，而如果采用類似kNN的方法就可能不太可行。
原始的資料集經過預處理之後儲存成兩個檔案：horseColicTest.txt和horseColicTraining.txt。
5.3.2 測試算法：用logistic回歸進行分類
使用logistic回歸方法進行分類要做的是把測試集上每個特征向量乘以最優化方法得來的回歸系數，再将該乘積結果求和，最後輸入到sigmoid函數中即可。如果對應的sigmoid值大于0.5就預測類别标簽為1，否則為0.
将下列代碼添加到logRegres.py中：
def classifyVector(inX,weights):
prob=sigmoid(sum(inX*weights))
if prob>0.5:
return 1.0
else:
return 0.0
def colicTest():
frTrain=open('./data/horseColicTraining.txt')
frTest=open('./data/horseColicTest.txt')
trainingSet=[]
trainingLabels=[]
for line in frTrain.readlines():
currLine=line.strip().split(' ')
lineArr=[]
for i in range(27):
lineArr.append(float(currLine[i]))
trainingSet.append(lineArr)
trainingLabels.append(float(currLine[27]))
trainWeights=stocGradAscent1(array(trainingSet),trainingLabels,500)
errorCount=0
numTestVec=0.0
for line in frTest.readlines():
numTestVec+=1.0
currLine=line.strip().split(' ')
lineArr=[]
for i in range(27):
lineArr.append(float(currLine[i]))
if int(classifyVector(array(lineArr),trainWeights))!=int(currLine[27]):
errorCount+=1
errorRate=(float(errorCount)/numTestVec)
print('the error rate of this test is: %f'%errorRate)
return errorRate
def multiTest():
numTests=10
errorSum=0.0
for k in range(numTests):
errorSum+=colicTest()
print('after %d iterations the average error rate is: %f'%(numTests,float(errorSum)/float(numTests)))
第一個函數classifyVector（），它以回歸系數和特征向量作為輸入來計算對應的sigmoid值。如果sigmoid值大于0.5函數傳回1，否則傳回0.
接下來的函數colicTest（）是用于打開測試集和訓練集，并對資料進行格式處理的函數。該函數首先導入訓練集，最後一列是類别标簽。資料最初有三個類别标簽，分别代表馬的三種情況：‘仍存活’、‘已經死亡’和‘已經安樂死’，這裡為了友善，将‘已經死亡’和‘已經安樂死’合并成‘未能存活’這個标簽。資料導入之後，便可以使用stocGradAscent1（）來計算回歸系數向量。這裡可以自由設定疊代的次數，例如在訓練集上使用500次疊代，實驗結果表明這比預設疊代150次的效果更好。在系數計算完成之後，導入測試集并計算分類錯誤率。整體看來，colicTest（）具有完全獨立的功能，多次運作得到的結果可能稍有不同，這是因為其中有随機的成分在裡面。如果在stocGradAscent1（）函數中回歸系數已經完全收斂，那麼結果才将是确定的。
最後一個函數是multiTest（），其功能是調用函數colicTest（）10次并求結果的平均值。
運作代碼：
if __name__=="__main__":
multiTest()
運作結果（可以通過調整colicTest（）中的疊代次數和stocGradAscent（）中的步長改變平均錯誤率）：
the error rate of this test is: 0.632353
the error rate of this test is: 0.632353
the error rate of this test is: 0.632353
the error rate of this test is: 0.632353
the error rate of this test is: 0.632353
the error rate of this test is: 0.632353
the error rate of this test is: 0.632353
the error rate of this test is: 0.632353
the error rate of this test is: 0.632353
the error rate of this test is: 0.632353
after 10 iterations the average error rate is: 0.632353
※由于每次的結果都一樣，是以可能有些錯誤，但是還未找到具體錯誤。
5.4 本章小結
logistic回歸的目的是尋找一個非線性函數sigmoid（）的最佳拟合參數，求解過程可以由最優化算法來完成。在最優化算法中，最常用的是梯度上升算法，而梯度上升算法又可以簡化為随機梯度上升算法。
随機梯度上升算法與梯度上升算法的效果相當，但占用更少的計算資源。此外，随機梯度上升是一個線上算法，它可以在新資料到來時就完成參數更新，而不需要重新讀取整個資料集來進行批處理運算。
機器學習的一個重要問題就是如何處理缺失資料。這個問題沒有标準答案，取決于實際應用中的需求。現有一些解決方案，每種方案都各有優缺點。