Python3《機器學習實戰》學習筆記（七）：Logistic回歸實戰篇之預測病馬死亡率

02-02

轉載請註明作者和出處： https://zhuanlan.zhihu.com/ml-jack

機器學習知乎專欄：https://zhuanlan.zhihu.com/ml-jack

CSDN博客專欄：http://blog.csdn.net/column/details/16415.html

Github代碼獲取：https://github.com/Jack-Cherish/Machine-Learning/

Python版本： Python3.x

運行平台： Windows

IDE： Sublime text3

一、前言

本文對梯度上升演算法和改進的隨機梯度上升演算法進行了對比，總結了各自的優缺點，並對sklearn.linear_model.LogisticRegression進行了詳細介紹。

二、改進的隨機梯度上升演算法

梯度上升演算法在每次更新回歸係數(最優參數)時，都需要遍歷整個數據集。可以看一下我們之前寫的梯度上升演算法：

def gradAscent(dataMatIn, classLabels):n dataMatrix = np.mat(dataMatIn) #轉換成numpy的matn labelMat = np.mat(classLabels).transpose() #轉換成numpy的mat,並進行轉置n m, n = np.shape(dataMatrix) #返回dataMatrix的大小。m為行數,n為列數。n alpha = 0.01 #移動步長,也就是學習速率,控制更新的幅度。n maxCycles = 500 #最大迭代次數n weights = np.ones((n,1))n for k in range(maxCycles):n h = sigmoid(dataMatrix * weights) #梯度上升矢量化公式n error = labelMat - hn weights = weights + alpha * dataMatrix.transpose() * errorn return weights.getA(),weights_array #將矩陣轉換為數組，返回權重數組n

假設，我們使用的數據集一共有100個樣本。那麼，dataMatrix就是一個100*3的矩陣。每次計算h的時候，都要計算dataMatrix*weights這個矩陣乘法運算，要進行100*3次乘法運算和100*2次加法運算。同理，更新回歸係數(最優參數)weights時，也需要用到整個數據集，要進行矩陣乘法運算。總而言之，該方法處理100個左右的數據集時尚可，但如果有數十億樣本和成千上萬的特徵，那麼該方法的計算複雜度就太高了。因此，需要對演算法進行改進，我們每次更新回歸係數(最優參數)的時候，能不能不用所有樣本呢？一次只用一個樣本點去更新回歸係數(最優參數)？這樣就可以有效減少計算量了，這種方法就叫做隨機梯度上升演算法。

1、隨機梯度上升演算法

讓我們直接看代碼：

def stocGradAscent1(dataMatrix, classLabels, numIter=150):n m,n = np.shape(dataMatrix) #返回dataMatrix的大小。m為行數,n為列數。n weights = np.ones(n) #參數初始化n for j in range(numIter): n dataIndex = list(range(m))n for i in range(m): n alpha = 4/(1.0+j+i)+0.01 #降低alpha的大小，每次減小1/(j+i)。n randIndex = int(random.uniform(0,len(dataIndex))) #隨機選取樣本n h = sigmoid(sum(dataMatrix[randIndex]*weights)) #選擇隨機選取的一個樣本，計算hn error = classLabels[randIndex] - h #計算誤差n weights = weights + alpha * error * dataMatrix[randIndex] #更新回歸係數n del(dataIndex[randIndex]) #刪除已經使用的樣本n return weights #返回n

該演算法第一個改進之處在於，alpha在每次迭代的時候都會調整，並且，雖然alpha會隨著迭代次數不斷減小，但永遠不會減小到0，因為這裡還存在一個常數項。必須這樣做的原因是為了保證在多次迭代之後新數據仍然具有一定的影響。如果需要處理的問題是動態變化的，那麼可以適當加大上述常數項，來確保新的值獲得更大的回歸係數。另一點值得注意的是，在降低alpha的函數中，alpha每次減少1/(j+i)，其中j是迭代次數，i是樣本點的下標。第二個改進的地方在於更新回歸係數(最優參數)時，只使用一個樣本點，並且選擇的樣本點是隨機的，每次迭代不使用已經用過的樣本點。這樣的方法，就有效地減少了計算量，並保證了回歸效果。

編寫代碼如下，看下改進的隨機梯度上升演算法分類效果如何：

# -*- coding:UTF-8 -*-nfrom matplotlib.font_manager import FontPropertiesnimport matplotlib.pyplot as pltnimport numpy as npnimport randomnn"""n函數說明:載入數據nnParameters:n 無nReturns:n dataMat - 數據列表n labelMat - 標籤列表nAuthor:n Jack CuinBlog:n http://blog.csdn.net/c406495762nZhihu:n https://www.zhihu.com/people/Jack--Cui/nModify:n 2017-08-28n"""ndef loadDataSet():n dataMat = [] #創建數據列表n labelMat = [] #創建標籤列表n fr = open(testSet.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 fr.close() #關閉文件n return dataMat, labelMat #返回nn"""n函數說明:sigmoid函數nnParameters:n inX - 數據nReturns:n sigmoid函數nAuthor:n Jack CuinBlog:n http://blog.csdn.net/c406495762nZhihu:n https://www.zhihu.com/people/Jack--Cui/nModify:n 2017-08-28n"""ndef sigmoid(inX):n return 1.0 / (1 + np.exp(-inX))nn"""n函數說明:繪製數據集nnParameters:n weights - 權重參數數組nReturns:n 無nAuthor:n Jack CuinBlog:n http://blog.csdn.net/c406495762nZhihu:n https://www.zhihu.com/people/Jack--Cui/nModify:n 2017-08-30n"""ndef plotBestFit(weights):n dataMat, labelMat = loadDataSet() #載入數據集n dataArr = np.array(dataMat) #轉換成numpy的array數組n n = np.shape(dataMat)[0] #數據個數n xcord1 = []; ycord1 = [] #正樣本n xcord2 = []; ycord2 = [] #負樣本n for i in range(n): #根據數據集標籤進行分類n if int(labelMat[i]) == 1:n xcord1.append(dataArr[i,1]); ycord1.append(dataArr[i,2]) #1為正樣本n else:n xcord2.append(dataArr[i,1]); ycord2.append(dataArr[i,2]) #0為負樣本n fig = plt.figure()n ax = fig.add_subplot(111) #添加subplotn ax.scatter(xcord1, ycord1, s = 20, c = red, marker = s,alpha=.5)#繪製正樣本n ax.scatter(xcord2, ycord2, s = 20, c = green,alpha=.5) #繪製負樣本n x = np.arange(-3.0, 3.0, 0.1)n y = (-weights[0] - weights[1] * x) / weights[2]n ax.plot(x, y)n plt.title(BestFit) #繪製titlen plt.xlabel(X1); plt.ylabel(X2) #繪製labeln plt.show()nn"""n函數說明:改進的隨機梯度上升演算法nnParameters:n dataMatrix - 數據數組n classLabels - 數據標籤n numIter - 迭代次數nReturns:n weights - 求得的回歸係數數組(最優參數)nAuthor:n Jack CuinBlog:n http://blog.csdn.net/c406495762nZhihu:n https://www.zhihu.com/people/Jack--Cui/nModify:n 2017-08-31n"""ndef stocGradAscent1(dataMatrix, classLabels, numIter=150):n m,n = np.shape(dataMatrix) #返回dataMatrix的大小。m為行數,n為列數。n weights = np.ones(n) #參數初始化n for j in range(numIter):n dataIndex = list(range(m))n for i in range(m):n alpha = 4/(1.0+j+i)+0.01 #降低alpha的大小，每次減小1/(j+i)。n randIndex = int(random.uniform(0,len(dataIndex))) #隨機選取樣本n h = sigmoid(sum(dataMatrix[randIndex]*weights)) #選擇隨機選取的一個樣本，計算hn error = classLabels[randIndex] - h #計算誤差n weights = weights + alpha * error * dataMatrix[randIndex] #更新回歸係數n del(dataIndex[randIndex]) #刪除已經使用的樣本n return weights #返回nnif __name__ == __main__:n dataMat, labelMat = loadDataSet()n weights = stocGradAscent1(np.array(dataMat), labelMat)n plotBestFit(weights)n

代碼運行結果：

2、回歸係數與迭代次數的關係

可以看到分類效果也是不錯的。不過，從這個分類結果中，我們不好看出迭代次數和回歸係數的關係，也就不能直觀的看到每個回歸方法的收斂情況。因此，我們編寫程序，繪製出回歸係數和迭代次數的關係曲線：

# -*- coding:UTF-8 -*-nfrom matplotlib.font_manager import FontPropertiesnimport matplotlib.pyplot as pltnimport numpy as npnimport randomnnn"""n函數說明:載入數據nnParameters:n 無nReturns:n dataMat - 數據列表n labelMat - 標籤列表nAuthor:n Jack CuinBlog:n http://blog.csdn.net/c406495762nZhihu:n https://www.zhihu.com/people/Jack--Cui/nModify:n 2017-08-28n"""ndef loadDataSet():n dataMat = [] #創建數據列表n labelMat = [] #創建標籤列表n fr = open(testSet.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 fr.close() #關閉文件n return dataMat, labelMat #返回nn"""n函數說明:sigmoid函數nnParameters:n inX - 數據nReturns:n sigmoid函數nAuthor:n Jack CuinBlog:n http://blog.csdn.net/c406495762nZhihu:n https://www.zhihu.com/people/Jack--Cui/nModify:n 2017-08-28n"""ndef sigmoid(inX):n return 1.0 / (1 + np.exp(-inX))nn"""n函數說明:梯度上升演算法nnParameters:n dataMatIn - 數據集n classLabels - 數據標籤nReturns:n weights.getA() - 求得的權重數組(最優參數)n weights_array - 每次更新的回歸係數nAuthor:n Jack CuinBlog:n http://blog.csdn.net/c406495762nZhihu:n https://www.zhihu.com/people/Jack--Cui/nModify:n 2017-08-28n"""ndef gradAscent(dataMatIn, classLabels):n dataMatrix = np.mat(dataMatIn) #轉換成numpy的matn labelMat = np.mat(classLabels).transpose() #轉換成numpy的mat,並進行轉置n m, n = np.shape(dataMatrix) #返回dataMatrix的大小。m為行數,n為列數。n alpha = 0.01 #移動步長,也就是學習速率,控制更新的幅度。n maxCycles = 500 #最大迭代次數n weights = np.ones((n,1))n weights_array = np.array([])n for k in range(maxCycles):n h = sigmoid(dataMatrix * weights) #梯度上升矢量化公式n error = labelMat - hn weights = weights + alpha * dataMatrix.transpose() * errorn weights_array = np.append(weights_array,weights)n weights_array = weights_array.reshape(maxCycles,n)n return weights.getA(),weights_array #將矩陣轉換為數組，並返回nnnn"""n函數說明:改進的隨機梯度上升演算法nnParameters:n dataMatrix - 數據數組n classLabels - 數據標籤n numIter - 迭代次數nReturns:n weights - 求得的回歸係數數組(最優參數)n weights_array - 每次更新的回歸係數nAuthor:n Jack CuinBlog:n http://blog.csdn.net/c406495762nZhihu:n https://www.zhihu.com/people/Jack--Cui/nModify:n 2017-08-31n"""ndef stocGradAscent1(dataMatrix, classLabels, numIter=150):n m,n = np.shape(dataMatrix) #返回dataMatrix的大小。m為行數,n為列數。n weights = np.ones(n) #參數初始化n weights_array = np.array([]) #存儲每次更新的回歸係數n for j in range(numIter): n dataIndex = list(range(m))n for i in range(m): n alpha = 4/(1.0+j+i)+0.01 #降低alpha的大小，每次減小1/(j+i)。n randIndex = int(random.uniform(0,len(dataIndex))) #隨機選取樣本n h = sigmoid(sum(dataMatrix[randIndex]*weights)) #選擇隨機選取的一個樣本，計算hn error = classLabels[randIndex] - h #計算誤差n weights = weights + alpha * error * dataMatrix[randIndex] #更新回歸係數n weights_array = np.append(weights_array,weights,axis=0) #添加回歸係數到數組中n del(dataIndex[randIndex]) #刪除已經使用的樣本n weights_array = weights_array.reshape(numIter*m,n) #改變維度n return weights,weights_array #返回nn"""n函數說明:繪製回歸係數與迭代次數的關係nnParameters:n weights_array1 - 回歸係數數組1n weights_array2 - 回歸係數數組2nReturns:n 無nAuthor:n Jack CuinBlog:n http://blog.csdn.net/c406495762nZhihu:n https://www.zhihu.com/people/Jack--Cui/nModify:n 2017-08-30n"""ndef plotWeights(weights_array1,weights_array2):n #設置漢字格式n font = FontProperties(fname=r"c:windowsfontssimsun.ttc", size=14)n #將fig畫布分隔成1行1列,不共享x軸和y軸,fig畫布的大小為(13,8)n #當nrow=3,nclos=2時,代表fig畫布被分為六個區域,axs[0][0]表示第一行第一列n fig, axs = plt.subplots(nrows=3, ncols=2,sharex=False, sharey=False, figsize=(20,10))n x1 = np.arange(0, len(weights_array1), 1)n #繪製w0與迭代次數的關係n axs[0][0].plot(x1,weights_array1[:,0])n axs0_title_text = axs[0][0].set_title(u梯度上升演算法：回歸係數與迭代次數關係,FontProperties=font)n axs0_ylabel_text = axs[0][0].set_ylabel(uW0,FontProperties=font)n plt.setp(axs0_title_text, size=20, weight=bold, color=black) n plt.setp(axs0_ylabel_text, size=20, weight=bold, color=black)n #繪製w1與迭代次數的關係n axs[1][0].plot(x1,weights_array1[:,1])n axs1_ylabel_text = axs[1][0].set_ylabel(uW1,FontProperties=font)n plt.setp(axs1_ylabel_text, size=20, weight=bold, color=black)n #繪製w2與迭代次數的關係n axs[2][0].plot(x1,weights_array1[:,2])n axs2_xlabel_text = axs[2][0].set_xlabel(u迭代次數,FontProperties=font)n axs2_ylabel_text = axs[2][0].set_ylabel(uW1,FontProperties=font)n plt.setp(axs2_xlabel_text, size=20, weight=bold, color=black) n plt.setp(axs2_ylabel_text, size=20, weight=bold, color=black)nnn x2 = np.arange(0, len(weights_array2), 1)n #繪製w0與迭代次數的關係n axs[0][1].plot(x2,weights_array2[:,0])n axs0_title_text = axs[0][1].set_title(u改進的隨機梯度上升演算法：回歸係數與迭代次數關係,FontProperties=font)n axs0_ylabel_text = axs[0][1].set_ylabel(uW0,FontProperties=font)n plt.setp(axs0_title_text, size=20, weight=bold, color=black) n plt.setp(axs0_ylabel_text, size=20, weight=bold, color=black)n #繪製w1與迭代次數的關係n axs[1][1].plot(x2,weights_array2[:,1])n axs1_ylabel_text = axs[1][1].set_ylabel(uW1,FontProperties=font)n plt.setp(axs1_ylabel_text, size=20, weight=bold, color=black)n #繪製w2與迭代次數的關係n axs[2][1].plot(x2,weights_array2[:,2])n axs2_xlabel_text = axs[2][1].set_xlabel(u迭代次數,FontProperties=font)n axs2_ylabel_text = axs[2][1].set_ylabel(uW1,FontProperties=font)n plt.setp(axs2_xlabel_text, size=20, weight=bold, color=black) n plt.setp(axs2_ylabel_text, size=20, weight=bold, color=black)nn plt.show() nnif __name__ == __main__:n dataMat, labelMat = loadDataSet() n weights1,weights_array1 = stocGradAscent1(np.array(dataMat), labelMat)nn weights2,weights_array2 = gradAscent(dataMat, labelMat)n plotWeights(weights_array1, weights_array2)n

運行結果如下：

由於改進的隨機梯度上升演算法，隨機選取樣本點，所以每次的運行結果是不同的。但是大體趨勢是一樣的。我們改進的隨機梯度上升演算法收斂效果更好。為什麼這麼說呢？讓我們分析一下。我們一共有100個樣本點，改進的隨機梯度上升演算法迭代次數為150。而上圖顯示15000次迭代次數的原因是，使用一次樣本就更新一下回歸係數。因此，迭代150次，相當於更新回歸係數150*100=15000次。簡而言之，迭代150次，更新1.5萬次回歸參數。從上圖左側的改進隨機梯度上升演算法回歸效果中可以看出，其實在更新2000次回歸係數的時候，已經收斂了。相當於遍歷整個數據集20次的時候，回歸係數已收斂。訓練已完成。

再讓我們看看上圖右側的梯度上升演算法回歸效果，梯度上升演算法每次更新回歸係數都要遍歷整個數據集。從圖中可以看出，當迭代次數為300多次的時候，回歸係數才收斂。湊個整，就當它在遍歷整個數據集300次的時候已經收斂好了。

沒有對比就沒有傷害，改進的隨機梯度上升演算法，在遍曆數據集的第20次開始收斂。而梯度上升演算法，在遍曆數據集的第300次才開始收斂。想像一下，大量數據的情況下，誰更牛逼？

三、從疝氣病癥狀預測病馬的死亡率

1、實戰背景

本次實戰內容，將使用Logistic回歸來預測患疝氣病的馬的存活問題。原始數據集下載地址：http://archive.ics.uci.edu/ml/datasets/Horse+Colic

這裡的數據包含了368個樣本和28個特徵。這種病不一定源自馬的腸胃問題，其他問題也可能引發馬疝病。該數據集中包含了醫院檢測馬疝病的一些指標，有的指標比較主觀，有的指標難以測量，例如馬的疼痛級別。另外需要說明的是，除了部分指標主觀和難以測量外，該數據還存在一個問題，數據集中有30%的值是缺失的。下面將首先介紹如何處理數據集中的數據缺失問題，然後再利用Logistic回歸和隨機梯度上升演算法來預測病馬的生死。

2、準備數據

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

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

預處理數據做兩件事：

如果測試集中一條數據的特徵值已經確實，那麼我們選擇實數0來替換所有缺失值，因為本文使用Logistic回歸。因此這樣做不會影響回歸係數的值。sigmoid(0)=0.5，即它對結果的預測不具有任何傾向性。
如果測試集中一條數據的類別標籤已經缺失，那麼我們將該類別數據丟棄，因為類別標籤與特徵不同，很難確定採用某個合適的值來替換。

原始的數據集經過處理，保存為兩個文件：horseColicTest.txt和horseColicTraining.txt。已經處理好的「乾淨」可用的數據集下載地址：

https://github.com/Jack-Cherish/Machine-Learning/blob/master/Logistic/horseColicTraining.txt
https://github.com/Jack-Cherish/Machine-Learning/blob/master/Logistic/horseColicTest.txt

有了這些數據集，我們只需要一個Logistic分類器，就可以利用該分類器來預測病馬的生死問題了。

3、使用Python構建Logistic回歸分類器

在使用Sklearn構建Logistic回歸分類器之前，我們先用自己寫的改進的隨機梯度上升演算法進行預測，先熱熱身。使用Logistic回歸方法進行分類並不需要做很多工作，所需做的只是把測試集上每個特徵向量乘以最優化方法得來的回歸係數，再將乘積結果求和，最後輸入到Sigmoid函數中即可。如果對應的Sigmoid值大於0.5就預測類別標籤為1，否則為0。

# -*- coding:UTF-8 -*-nimport numpy as npnimport randomnn"""n函數說明:sigmoid函數nnParameters:n inX - 數據nReturns:n sigmoid函數nAuthor:n Jack CuinBlog:n http://blog.csdn.net/c406495762nZhihu:n https://www.zhihu.com/people/Jack--Cui/nModify:n 2017-09-05n"""ndef sigmoid(inX):n return 1.0 / (1 + np.exp(-inX))nn"""n函數說明:改進的隨機梯度上升演算法nnParameters:n dataMatrix - 數據數組n classLabels - 數據標籤n numIter - 迭代次數nReturns:n weights - 求得的回歸係數數組(最優參數)nAuthor:n Jack CuinBlog:n http://blog.csdn.net/c406495762nZhihu:n https://www.zhihu.com/people/Jack--Cui/nModify:n 2017-09-05n"""ndef stocGradAscent1(dataMatrix, classLabels, numIter=150):n m,n = np.shape(dataMatrix) #返回dataMatrix的大小。m為行數,n為列數。n weights = np.ones(n) #參數初始化 #存儲每次更新的回歸係數n for j in range(numIter): n dataIndex = list(range(m))n for i in range(m): n alpha = 4/(1.0+j+i)+0.01 #降低alpha的大小，每次減小1/(j+i)。n randIndex = int(random.uniform(0,len(dataIndex))) #隨機選取樣本n h = sigmoid(sum(dataMatrix[randIndex]*weights)) #選擇隨機選取的一個樣本，計算hn error = classLabels[randIndex] - h #計算誤差n weights = weights + alpha * error * dataMatrix[randIndex] #更新回歸係數n del(dataIndex[randIndex]) #刪除已經使用的樣本n return weights #返回nn"""n函數說明:使用Python寫的Logistic分類器做預測nnParameters:n 無nReturns:n 無nAuthor:n Jack CuinBlog:n http://blog.csdn.net/c406495762nZhihu:n https://www.zhihu.com/people/Jack--Cui/nModify:n 2017-09-05n"""ndef colicTest():n frTrain = open(horseColicTraining.txt) #打開訓練集n frTest = open(horseColicTest.txt) #打開測試集n trainingSet = []; trainingLabels = []n for line in frTrain.readlines():n currLine = line.strip().split(t)n lineArr = []n for i in range(len(currLine)-1):n lineArr.append(float(currLine[i]))n trainingSet.append(lineArr)n trainingLabels.append(float(currLine[-1]))n trainWeights = stocGradAscent1(np.array(trainingSet), trainingLabels, 500) #使用改進的隨即上升梯度訓練n errorCount = 0; numTestVec = 0.0n for line in frTest.readlines():n numTestVec += 1.0n currLine = line.strip().split(t)n lineArr =[]n for i in range(len(currLine)-1):n lineArr.append(float(currLine[i]))n if int(classifyVector(np.array(lineArr), trainWeights))!= int(currLine[-1]):n errorCount += 1n errorRate = (float(errorCount)/numTestVec) * 100 #錯誤率計算n print("測試集錯誤率為: %.2f%%" % errorRate)nn"""n函數說明:分類函數nnParameters:n inX - 特徵向量n weights - 回歸係數nReturns:n 分類結果nAuthor:n Jack CuinBlog:n http://blog.csdn.net/c406495762nZhihu:n https://www.zhihu.com/people/Jack--Cui/nModify:n 2017-09-05n"""ndef classifyVector(inX, weights):n prob = sigmoid(sum(inX*weights))n if prob > 0.5: return 1.0n else: return 0.0nnif __name__ == __main__:n colicTest()n

運行結果如下：

錯誤率還是蠻高的，而且耗時1.9s，並且每次運行的錯誤率也是不同的，錯誤率高的時候可能達到40%多。為啥這樣？首先，因為數據集本身有30%的數據缺失，這個是不能避免的。另一個主要原因是，我們使用的是改進的隨機梯度上升演算法，因為數據集本身就很小，就幾百的數據量。用改進的隨機梯度上升演算法顯然不合適。讓我們再試試梯度上升演算法，看看它的效果如何？

# -*- coding:UTF-8 -*-nimport numpy as npnimport randomnn"""n函數說明:sigmoid函數nnParameters:n inX - 數據nReturns:n sigmoid函數nAuthor:n Jack CuinBlog:n http://blog.csdn.net/c406495762nZhihu:n https://www.zhihu.com/people/Jack--Cui/nModify:n 2017-09-05n"""ndef sigmoid(inX):n return 1.0 / (1 + np.exp(-inX))nn"""n函數說明:梯度上升演算法nnParameters:n dataMatIn - 數據集n classLabels - 數據標籤nReturns:n weights.getA() - 求得的權重數組(最優參數)nAuthor:n Jack CuinBlog:n http://blog.csdn.net/c406495762nZhihu:n https://www.zhihu.com/people/Jack--Cui/nModify:n 2017-08-28n"""ndef gradAscent(dataMatIn, classLabels):n dataMatrix = np.mat(dataMatIn) #轉換成numpy的matn labelMat = np.mat(classLabels).transpose() #轉換成numpy的mat,並進行轉置n m, n = np.shape(dataMatrix) #返回dataMatrix的大小。m為行數,n為列數。n alpha = 0.01 #移動步長,也就是學習速率,控制更新的幅度。n maxCycles = 500 #最大迭代次數n weights = np.ones((n,1))n for k in range(maxCycles):n h = sigmoid(dataMatrix * weights) #梯度上升矢量化公式n error = labelMat - hn weights = weights + alpha * dataMatrix.transpose() * errorn return weights.getA() #將矩陣轉換為數組，並返回nnnn"""n函數說明:使用Python寫的Logistic分類器做預測nnParameters:n 無nReturns:n 無nAuthor:n Jack CuinBlog:n http://blog.csdn.net/c406495762nZhihu:n https://www.zhihu.com/people/Jack--Cui/nModify:n 2017-09-05n"""ndef colicTest():n frTrain = open(horseColicTraining.txt) #打開訓練集n frTest = open(horseColicTest.txt) #打開測試集n trainingSet = []; trainingLabels = []n for line in frTrain.readlines():n currLine = line.strip().split(t)n lineArr = []n for i in range(len(currLine)-1):n lineArr.append(float(currLine[i]))n trainingSet.append(lineArr)n trainingLabels.append(float(currLine[-1]))n trainWeights = gradAscent(np.array(trainingSet), trainingLabels) #使用改進的隨即上升梯度訓練n errorCount = 0; numTestVec = 0.0n for line in frTest.readlines():n numTestVec += 1.0n currLine = line.strip().split(t)n lineArr =[]n for i in range(len(currLine)-1):n lineArr.append(float(currLine[i]))n if int(classifyVector(np.array(lineArr), trainWeights[:,0]))!= int(currLine[-1]):n errorCount += 1n errorRate = (float(errorCount)/numTestVec) * 100 #錯誤率計算n print("測試集錯誤率為: %.2f%%" % errorRate)nn"""n函數說明:分類函數nnParameters:n inX - 特徵向量n weights - 回歸係數nReturns:n 分類結果nAuthor:n Jack CuinBlog:n http://blog.csdn.net/c406495762nZhihu:n https://www.zhihu.com/people/Jack--Cui/nModify:n 2017-09-05n"""ndef classifyVector(inX, weights):n prob = sigmoid(sum(inX*weights))n if prob > 0.5: return 1.0n else: return 0.0nnif __name__ == __main__:n colicTest()n

運行結果如下：

可以看到演算法耗時減少了，錯誤率穩定且較低。很顯然，使用隨機梯度上升演算法，反而得不償失了。所以可以得到如下結論：

當數據集較小時，我們使用梯度上升演算法
當數據集較大時，我們使用改進的隨機梯度上升演算法

對應的，在Sklearn中，我們就可以根據數據情況選擇優化演算法，比如數據較小的時候，我們使用liblinear，數據較大時，我們使用sag和saga。

四、使用Sklearn構建Logistic回歸分類器

開始新一輪的征程，讓我們看下Sklearn的Logistic回歸分類器！

官方英文文檔地址：http://scikit-learn.org/dev/modules/generated/sklearn.linear_model.LogisticRegression.html#sklearn.linear_model.LogisticRegression

sklearn.linear_model模塊提供了很多模型供我們使用，比如Logistic回歸、Lasso回歸、貝葉斯脊回歸等，可見需要學習的東西還有很多很多。本篇文章，我們使用LogisticRegressioin。

1、LogisticRegression

讓我們先看下LogisticRegression這個函數，一共有14個參數：

參數說明如下：

penalty：懲罰項，str類型，可選參數為l1和l2，默認為l2。用於指定懲罰項中使用的規範。newton-cg、sag和lbfgs求解演算法只支持L2規範。L1G規範假設的是模型的參數滿足拉普拉斯分布，L2假設的模型參數滿足高斯分布，所謂的範式就是加上對參數的約束，使得模型更不會過擬合(overfit)，但是如果要說是不是加了約束就會好，這個沒有人能回答，只能說，加約束的情況下，理論上應該可以獲得泛化能力更強的結果。
dual：對偶或原始方法，bool類型，默認為False。對偶方法只用在求解線性多核(liblinear)的L2懲罰項上。當樣本數量>樣本特徵的時候，dual通常設置為False。
tol：停止求解的標準，float類型，默認為1e-4。就是求解到多少的時候，停止，認為已經求出最優解。
c：正則化係數λ的倒數，float類型，默認為1.0。必須是正浮點型數。像SVM一樣，越小的數值表示越強的正則化。
fit_intercept：是否存在截距或偏差，bool類型，默認為True。
intercept_scaling：僅在正則化項為"liblinear"，且fit_intercept設置為True時有用。float類型，默認為1。
class_weight：用於標示分類模型中各種類型的權重，可以是一個字典或者balanced字元串，默認為不輸入，也就是不考慮權重，即為None。如果選擇輸入的話，可以選擇balanced讓類庫自己計算類型權重，或者自己輸入各個類型的權重。舉個例子，比如對於0,1的二元模型，我們可以定義class_weight={0:0.9,1:0.1}，這樣類型0的權重為90%，而類型1的權重為10%。如果class_weight選擇balanced，那麼類庫會根據訓練樣本量來計算權重。某種類型樣本量越多，則權重越低，樣本量越少，則權重越高。當class_weight為balanced時，類權重計算方法如下：n_samples / (n_classes * np.bincount(y))。n_samples為樣本數，n_classes為類別數量，np.bincount(y)會輸出每個類的樣本數，例如y=[1,0,0,1,1],則np.bincount(y)=[2,3]。

那麼class_weight有什麼作用呢？

在分類模型中，我們經常會遇到兩類問題：
1.第一種是誤分類的代價很高。比如對合法用戶和非法用戶進行分類，將非法用戶分類為合法用戶的代價很高，我們寧願將合法用戶分類為非法用戶，這時可以人工再甄別，但是卻不願將非法用戶分類為合法用戶。這時，我們可以適當提高非法用戶的權重。
2. 第二種是樣本是高度失衡的，比如我們有合法用戶和非法用戶的二元樣本數據10000條，裡面合法用戶有9995條，非法用戶只有5條，如果我們不考慮權重，則我們可以將所有的測試集都預測為合法用戶，這樣預測準確率理論上有99.95%，但是卻沒有任何意義。這時，我們可以選擇balanced，讓類庫自動提高非法用戶樣本的權重。提高了某種分類的權重，相比不考慮權重，會有更多的樣本分類劃分到高權重的類別，從而可以解決上面兩類問題。

random_state：隨機數種子，int類型，可選參數，默認為無，僅在正則化優化演算法為sag,liblinear時有用。
solver：優化演算法選擇參數，只有五個可選參數，即newton-cg,lbfgs,liblinear,sag,saga。默認為liblinear。solver參數決定了我們對邏輯回歸損失函數的優化方法，有四種演算法可以選擇，分別是：

liblinear：使用了開源的liblinear庫實現，內部使用了坐標軸下降法來迭代優化損失函數。
lbfgs：擬牛頓法的一種，利用損失函數二階導數矩陣即海森矩陣來迭代優化損失函數。
newton-cg：也是牛頓法家族的一種，利用損失函數二階導數矩陣即海森矩陣來迭代優化損失函數。
sag：即隨機平均梯度下降，是梯度下降法的變種，和普通梯度下降法的區別是每次迭代僅僅用一部分的樣本來計算梯度，適合於樣本數據多的時候。
saga：線性收斂的隨機優化演算法的的變重。
總結：

liblinear適用於小數據集，而sag和saga適用於大數據集因為速度更快。
對於多分類問題，只有newton-cg,sag,saga和lbfgs能夠處理多項損失，而liblinear受限於一對剩餘(OvR)。啥意思，就是用liblinear的時候，如果是多分類問題，得先把一種類別作為一個類別，剩餘的所有類別作為另外一個類別。一次類推，遍歷所有類別，進行分類。
newton-cg,sag和lbfgs這三種優化演算法時都需要損失函數的一階或者二階連續導數，因此不能用於沒有連續導數的L1正則化，只能用於L2正則化。而liblinear和saga通吃L1正則化和L2正則化。
同時，sag每次僅僅使用了部分樣本進行梯度迭代，所以當樣本量少的時候不要選擇它，而如果樣本量非常大，比如大於10萬，sag是第一選擇。但是sag不能用於L1正則化，所以當你有大量的樣本，又需要L1正則化的話就要自己做取捨了。要麼通過對樣本採樣來降低樣本量，要麼回到L2正則化。
從上面的描述，大家可能覺得，既然newton-cg, lbfgs和sag這麼多限制，如果不是大樣本，我們選擇liblinear不就行了嘛！錯，因為liblinear也有自己的弱點！我們知道，邏輯回歸有二元邏輯回歸和多元邏輯回歸。對於多元邏輯回歸常見的有one-vs-rest(OvR)和many-vs-many(MvM)兩種。而MvM一般比OvR分類相對準確一些。鬱悶的是liblinear只支持OvR，不支持MvM，這樣如果我們需要相對精確的多元邏輯回歸時，就不能選擇liblinear了。也意味著如果我們需要相對精確的多元邏輯回歸不能使用L1正則化了。

max_iter：演算法收斂最大迭代次數，int類型，默認為10。僅在正則化優化演算法為newton-cg, sag和lbfgs才有用，演算法收斂的最大迭代次數。
multi_class：分類方式選擇參數，str類型，可選參數為ovr和multinomial，默認為ovr。ovr即前面提到的one-vs-rest(OvR)，而multinomial即前面提到的many-vs-many(MvM)。如果是二元邏輯回歸，ovr和multinomial並沒有任何區別，區別主要在多元邏輯回歸上。

OvR和MvM有什麼不同？

OvR的思想很簡單，無論你是多少元邏輯回歸，我們都可以看做二元邏輯回歸。具體做法是，對於第K類的分類決策，我們把所有第K類的樣本作為正例，除了第K類樣本以外的所有樣本都作為負例，然後在上面做二元邏輯回歸，得到第K類的分類模型。其他類的分類模型獲得以此類推。
而MvM則相對複雜，這裡舉MvM的特例one-vs-one(OvO)作講解。如果模型有T類，我們每次在所有的T類樣本裡面選擇兩類樣本出來，不妨記為T1類和T2類，把所有的輸出為T1和T2的樣本放在一起，把T1作為正例，T2作為負例，進行二元邏輯回歸，得到模型參數。我們一共需要T(T-1)/2次分類。
可以看出OvR相對簡單，但分類效果相對略差（這裡指大多數樣本分布情況，某些樣本分布下OvR可能更好）。而MvM分類相對精確，但是分類速度沒有OvR快。如果選擇了ovr，則4種損失函數的優化方法liblinear，newton-cg,lbfgs和sag都可以選擇。但是如果選擇了multinomial,則只能選擇newton-cg, lbfgs和sag了。

verbose：日誌冗長度，int類型。默認為0。就是不輸出訓練過程，1的時候偶爾輸出結果，大於1，對於每個子模型都輸出。
warm_start：熱啟動參數，bool類型。默認為False。如果為True，則下一次訓練是以追加樹的形式進行（重新使用上一次的調用作為初始化）。
n_jobs：並行數。int類型，默認為1。1的時候，用CPU的一個內核運行程序，2的時候，用CPU的2個內核運行程序。為-1的時候，用所有CPU的內核運行程序。

累死我了....終於寫完所有參數了。

除此之外，LogisticRegression也有一些方法供我們使用：

有一些方法和MultinomialNB的方法都是類似的，因此不再累述。

對於每個函數的具體使用，可以看下官方文檔：http://scikit-learn.org/dev/modules/generated/sklearn.linear_model.LogisticRegression.html#sklearn.linear_model.LogisticRegression

同時，如果對於過擬合、正則化、L1範數、L2範數不了解的，可以看這位大牛的博客：http://blog.csdn.net/zouxy09/article/details/24971995

2、編寫代碼

了解到這些，我們就可以編寫Sklearn分類器的代碼了。代碼非常短：

# -*- coding:UTF-8 -*-nfrom sklearn.linear_model import LogisticRegressionnn"""n函數說明:使用Sklearn構建Logistic回歸分類器nnParameters:n 無nReturns:n 無nAuthor:n Jack CuinBlog:n http://blog.csdn.net/c406495762nZhihu:n https://www.zhihu.com/people/Jack--Cui/nModify:n 2017-09-05n"""ndef colicSklearn():n frTrain = open(horseColicTraining.txt) #打開訓練集n frTest = open(horseColicTest.txt) #打開測試集n trainingSet = []; trainingLabels = []n testSet = []; testLabels = []n for line in frTrain.readlines():n currLine = line.strip().split(t)n lineArr = []n for i in range(len(currLine)-1):n lineArr.append(float(currLine[i]))n trainingSet.append(lineArr)n trainingLabels.append(float(currLine[-1]))n for line in frTest.readlines():n currLine = line.strip().split(t)n lineArr =[]n for i in range(len(currLine)-1):n lineArr.append(float(currLine[i]))n testSet.append(lineArr)n testLabels.append(float(currLine[-1]))n classifier = LogisticRegression(solver=liblinear,max_iter=10).fit(trainingSet, trainingLabels)n test_accurcy = classifier.score(testSet, testLabels) * 100n print(正確率:%f%% % test_accurcy)nnif __name__ == __main__:n colicSklearn()n

運行結果如下：

可以看到，正確率又高一些了。更改solver參數，比如設置為sag，使用隨機平均梯度下降演算法，看一看效果。你會發現，有警告了。

顯而易見，警告是因為演算法還沒有收斂。更改max_iter=5000，再運行代碼：

可以看到，對於我們這樣的小數據集，sag演算法需要迭代上千次才收斂，而liblinear只需要不到10次。

還是那句話，我們需要根據數據集情況，選擇最優化演算法。

五、總結

1、Logistic回歸的優缺點

優點：

實現簡單，易於理解和實現；計算代價不高，速度很快，存儲資源低。

缺點：

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

2、其他

Logistic回歸的目的是尋找一個非線性函數Sigmoid的最佳擬合參數，求解過程可以由最優化演算法完成。
改進的一些最優化演算法，比如sag。它可以在新數據到來時就完成參數更新，而不需要重新讀取整個數據集來進行批量處理。
機器學習的一個重要問題就是如何處理缺失數據。這個問題沒有標準答案，取決於實際應用中的需求。現有一些解決方案，每種方案都各有優缺點。
我們需要根據數據的情況，這是Sklearn的參數，以期達到更好的分類效果。
下篇文章將講解支持向量機SVM。
如有問題，請留言。如有錯誤，還望指正，謝謝！

PS：如果覺得本篇本章對您有所幫助，歡迎關注、評論、贊！

本文出現的所有代碼和數據集，均可在我的github上下載，歡迎Follow、Star：https://github.com/Jack-Cherish/Machine-Learning

參考文獻：

《機器學習實戰》的第五章內容。