機器學習之Logistic 回歸算法
1 Logistic 回歸算法的原理
1.1 需要的數學基礎
我在看機器學習實戰時對其中的代碼非常費解,說好的利用偏導數求最值怎麽代碼中沒有體現啊,就一個簡單的式子:θ= θ - α Σ [( hθ(x(i))-y(i) ) ] * xi 。經過查找資料才知道,書中省去了大量的理論推導過程,其中用到了線性函數、sigmoid 函數、偏導數、最大似然函數、梯度下降法。下面讓我們一窺究竟,是站在大神的肩膀描述我自己的見解。
1.2 Logistic 回歸的引入
Logistic 回歸是概率非線性模型,用於處理二元分類結果,名為回歸實為分類。下面給出一個二元分類的例子,如下圖所示:
圖中數據分成兩類,打叉的一類用 y = 1 表示,另一類為圓圈用 y= 0 表示。現在我們要對數據進行分類,要找到一個分類函數(邊界線),我們很容易得出一個線性函數對其分類:0 = θ0
Sigmoid 函數的值域為(0,1),且具有良好的從0 到 1 的跳躍性,如在兩個不同的坐標尺度下的函數圖形:
所以,我們把線性方程和Sigmoid 函數結合起來就能解決問題。即 :分類預測函數 hθ (x) = g( θ0 + θ1x1 + θ2x2) .我們就可以對樣本數據進行分類,如下圖所示:
對於線性的分類邊界,如下形式:
分類預測函數,如下形式:
其中,θT
1.3 分類預測函數問題的轉化成求θ
通過上面的分析,我們得出了分類預測函數 hθ(x) , 但其中向量 x 是已知的(x 是未知類別號的對象數據),向量 θ 未知,即我們把求分類函數問題轉化成求向量 θ 。因為Sigmoid 函數的取值區間(0,1),那我們可以看做概率 P(y = 1 | xi ; θ)= hθ(x) , 表示在 xi 確定的情況下,類別 y = 1 的概率。由此,我們也可以得出在 xi 確定的情況下,類別 y = 0 的概率 P(y = 0 | xi
我們可以將這兩個式子合並得:
其中的 y = 0 或 1 .
這時候我們可以利用最大似然函數對向量 θ 求值,可以理解為選取的樣本數據的概率是最大的,那麽樣本數為 m 的似然函數為:
通過對數的形式對似然函數進行變化,對數似然函數:
這裏的最大似然函數的值就是我們所要求的向量 θ , 而求解的方法利用梯度下降法。
1.4 梯度下降法求解θ
在用梯度下降法時,我們將會利用Sigmoid 函數的一個性質: g, (z) = g(z)[ 1- g(z) ] 。
構造一個Cost函數(損失函數),該函數表示預測的輸出(h)與訓練數據類別(y)之間的偏差,可以是二者之間的差(h-y)或者是其他的形式。綜合考慮所有訓練數據的“損失”,將Cost求和或者求平均,記為J(θ)函數,表示所有訓練數據預測值與實際類別的偏差。
損失函數:
J(θ)代價函數:
其中,x(i) 每個樣本數據點在某一個特征上的值,即特征向量x的某個值,y(i) 是類別號,m 是樣本對象個數。
梯度下降法含義:
梯度下降法,就是利用負梯度方向來決定每次叠代的新的搜索方向,使得每次叠代能使待優化的目標函數逐步減小。梯度其實就是函數的偏導數。
這裏對用梯度下降法對 J (θ) 求最小值,與求似然函數的最大值是一樣的。則 J(θ) 最小值的求解過程:
其中 α 是步長。
則可以得出:
因為 α 是個常量,所以一般情況可以把 1/m 省去,省去不是沒有用1/m ,只是看成 α 和1/m 合並成 α 。
最終式子為:
這就是一開始,我對代碼中公式困惑的地方。在這裏我在補充一點,以上的梯度下降法可以認為是批量梯度下降法(Batch Gradient Descent),由於我們有m個樣本,這裏求梯度的時候就用了所有m個樣本的梯度數據。下面介紹隨機梯度下降法(Stochastic Gradient Descent):
隨機梯度下降法,其實和批量梯度下降法原理類似,區別在與求梯度時沒有用所有的m個樣本的數據,而是僅僅選取一個樣本j來求梯度。隨機梯度下降法,和4.1的批量梯度下降法是兩個極端,一個采用所有數據來梯度下降,一個用一個樣本來梯度下降。自然各自的優缺點都非常突出。對於訓練速度來說,隨機梯度下降法由於每次僅僅采用一個樣本來叠代,訓練速度很快,而批量梯度下降法在樣本量很大的時候,訓練速度不能讓人滿意。對於準確度來說,隨機梯度下降法用於僅僅用一個樣本決定梯度方向,導致解很有可能不是最優。對於收斂速度來說,由於隨機梯度下降法一次叠代一個樣本,導致叠代方向變化很大,不能很快的收斂到局部最優解。公式為:
到這裏已經把Logistics 回歸的原理推導完成。這裏對以上推導做一次總結:
邊界線 ——> Sigmoid 函數 ——>求向量 θ ——> P(y=1 | x ; θ) = hθ(x) —— > 似然函數 l (θ) ——> 代價函數 J (θ) ——> 梯度下降法求解向量 θ ——> 最終公式 θj
2 使用python 實現Logistic 回歸
2.1 Logistic 回歸梯度上升優化算法
2.1.1 收集數據和處理數據
1 # 從文件中提取數據 2 def loadDataSet(): 3 dataMat = [] ; labelMat = [] 4 fr = open(‘testSet.txt‘) 5 for line in fr.readlines(): # 對文件的數據進行按行遍歷 6 lineArr = line.strip().split() 7 dataMat.append([1.0, float(lineArr [0]), float(lineArr[1])]) 8 labelMat.append(int(lineArr[2])) # 數據的類別號列表 9 return dataMat , labelMat
2.1.2 Sigmoid 函數
1 # 定義sigmoid 函數 2 def sigmoid(inX): 3 return longfloat( 1.0 / (1 + exp(-inX)))
使用longfloat() 是為防止溢出。
2.1.3 批量梯度上升算法的實現
1 # Logistic 回歸梯度上升優化算法 2 def gradAscent(dataMatIn, classLabels): 3 dataMatrix = mat(dataMatIn) # 將數列轉化成矩陣 4 labelMat = mat(classLabels).transpose() # 將類標號轉化成矩陣並轉置成列向量 5 m,n = shape(dataMatrix) # 獲得矩陣dataMatrix 的行、列數 6 alpha = 0.001 # 向目標移動的步長 7 maxCycles = 500 # 叠代次數 8 weights = ones((n ,1)) # 生成 n 行 1 列的 矩陣且值為1 9 for k in range(maxCycles): 10 h = sigmoid(dataMatrix*weights) # dataMatrix*weights 是 m * 1 的矩陣,其每一個元素都會調用sigmoid()函數,h 也是一個 m * 1 的矩陣 11 error = (labelMat - h) # 損失函數 12 weights = weights + alpha + dataMatrix.transpose()* error # 每步 weights 該變量 13 return weights.getA()
解釋第 13 行代碼:矩陣通過這個getA()這個方法可以將自身返回成一個n維數組對象 ,因為plotBestFit()函數中有計算散點x,y坐標的部分,其中計算y的時候用到了weights,如果weights是矩陣的話,weights[1]就是[0.48007329](註意這裏有中括號!),就不是一個數了,最終你會發現y的計算結果的len()只有1,而x的len()則是60。
2.2 根據批量梯度上升法分析數據:畫出決策邊界
1 # 對數據分類的邊界圖形顯示 2 def plotBestFit(weights): 3 import matplotlib.pyplot as plt 4 dataMat,labelMat=loadDataSet() 5 dataArr = array(dataMat) 6 n = shape(dataArr)[0] # 數據的行數,即對象的個數 7 xcord1 = []; ycord1 = [] # 對類別號為 1 的對象,分 X 軸和 Y 軸的數據 8 xcord2 = []; ycord2 = [] # 對類別號為 0 的對象,分 X 軸和 Y 軸的數據 9 for i in range(n): # 對所有的對象進行遍歷 10 if int(labelMat[i])== 1: # 對象的類別為:1 11 xcord1.append(dataArr[i,1]); ycord1.append(dataArr[i,2]) 12 else: 13 xcord2.append(dataArr[i,1]); ycord2.append(dataArr[i,2]) 14 fig = plt.figure() 15 ax = fig.add_subplot(111) 16 ax.scatter(xcord1, ycord1, s=30, c=‘red‘, marker=‘s‘) # 對散點的格式的設置,坐標號、點的大小、顏色、點的圖形(方塊) 17 ax.scatter(xcord2, ycord2, s=30, c=‘green‘) # 點的圖形默認為圓 18 x = arange(-3.0, 3.0, 0.1) 19 y = (-weights[0]-weights[1]*x)/weights[2] # 線性方程 y = aX + b,y 是數據第三列的特征,X 是數據第二列的特征 20 ax.plot(x, y) 21 plt.xlabel(‘X1‘); plt.ylabel(‘X2‘) 22 plt.show()
運行結果:
有結果效果圖可知,邊界線基本可以對樣本進行較好的分類。
2.3 利用隨機梯度上升法訓練樣本
1 # 隨機梯度上升法 2 def stocGradAscent0(dataMatrix, classLabels, numIter=150): 3 m,n = shape(dataMatrix ) 4 weights = ones(n) 5 for j in range(numIter): # 叠代次數 6 dataIndex = range(m) # 7 for i in range(m): # 對所有對象的遍歷 8 alpha = 4/(1.0+j+i)+0.01 # 對步長的調整 9 randIndex = int (random.uniform(0,len(dataIndex))) # 隨機生成一個整數,介於 0 到 m 10 h = sigmoid(sum(dataMatrix[randIndex]*weights)) # 對隨機選擇的對象計算類別的數值(回歸系數值) 11 error = classLabels[randIndex] - h # 根據實際類型與計算類型值的誤差,損失函數 12 weights = weights + alpha * error * dataMatrix[randIndex] # 每步 weights 改變值 13 del(dataIndex [randIndex]) # 去除已經選擇過的對象,避免下次選中 14 return weights
在隨機梯度上升法求解的 θ ,對樣本進行分類,並畫出其邊界線,運行的結果圖:
利用隨機梯度上升法,通過150 次叠代得出的效果圖和批量梯度上升法的效果圖差不過,但隨機梯度的效率比批量梯度上升法快很多。
3 實例:從疝氣病癥預測病馬的死亡率
3.1 未知對象的預測
1 # 對未知對象進行預測類別號 2 def classifyVector(inX , weights): 3 prob = sigmoid(sum(inX * weights)) # 計算回歸系數 4 if prob > 0.5: 5 return 1.0 6 else: 7 return 0.0
3.2 樣本數據和測試函數
1 # 實例:從疝氣病預測病馬的死亡率 2 def colicTest(): 3 frTrain = open(‘horseColicTraining.txt‘); frTest = open(‘horseColicTest.txt‘)# 數據文件 4 trainingSet = []; trainingLabels = [] # 樣本集和類標號集的初始化 5 for line in frTrain.readlines(): 6 currLine = line.strip().split(‘\t‘) # 根據制表符進行字符串的分割 7 lineArr =[] 8 for i in range(21): 9 lineArr.append(float(currLine[i])) 10 trainingSet.append(lineArr) 11 trainingLabels.append(float(currLine[21])) 12 #trainWeights = stocGradAscent0(array(trainingSet), trainingLabels, 500) # 通過對訓練樣本計算出回歸系數 13 trainWeights = gradAscent(array(trainingSet), trainingLabels) # 通過對訓練樣本計算出回歸系數 14 errorCount = 0; numTestVec = 0.0 # 錯誤個數和錯誤率的初始化 15 # 預測樣本的遍歷 16 for line in frTest.readlines(): 17 numTestVec += 1.0 18 currLine = line.strip().split(‘\t‘) 19 lineArr =[] 20 for i in range(21): 21 lineArr.append(float(currLine[i])) 22 if int(classifyVector(array(lineArr), trainWeights))!= int(currLine[21]): # 預測的類別號與真實類別號的比較 23 errorCount += 1 24 errorRate = (float(errorCount)/numTestVec) 25 print u‘本次測試的錯誤率是; %f‘ % errorRate 26 return errorRate
3.3 預測結果函數
1 # 預測結果函數 2 def multiTest(): 3 numTests = 10; errorSum = 0.0 4 for k in range(numTests): 5 errorSum += colicTest() 6 print u‘經過 %d 次測試結果的平均錯誤率是: %f ‘ % (numTests ,errorSum /float(numTests))
利用隨機梯度上升法訓練的樣本集,測試結果:
利用批量梯度上升法訓練的樣本集,測試結果:
附 完整程序
1 # coding:utf-8 2 from numpy import * 3 4 # 從文件中提取數據 5 def loadDataSet(): 6 dataMat = [] ; labelMat = [] 7 fr = open(‘testSet.txt‘) 8 for line in fr.readlines(): # 對文件的數據進行按行遍歷 9 lineArr = line.strip().split() 10 dataMat.append([1.0, float(lineArr [0]), float(lineArr[1])]) 11 labelMat.append(int(lineArr[2])) # 數據的類別號列表 12 return dataMat , labelMat 13 14 # 定義sigmoid 函數 15 def sigmoid(inX): 16 return longfloat( 1.0 / (1 + exp(-inX))) 17 18 # Logistic 回歸批量梯度上升優化算法 19 def gradAscent(dataMatIn, classLabels): 20 dataMatrix = mat(dataMatIn) # 將數列轉化成矩陣 21 labelMat = mat(classLabels).transpose() # 將類標號轉化成矩陣並轉置成列向量 22 m,n = shape(dataMatrix) # 獲得矩陣dataMatrix 的行、列數 23 alpha = 0.001 # 向目標移動的步長 24 maxCycles = 500 # 叠代次數 25 weights = ones((n ,1)) # 生成 n 行 1 列的 矩陣且值為1 26 for k in range(maxCycles): 27 h = sigmoid(dataMatrix*weights) # dataMatrix*weights 是 m * 1 的矩陣,其每一個元素都會調用sigmoid()函數,h 也是一個 m * 1 的矩陣 28 error = (labelMat - h) # 損失函數 29 weights = weights + alpha + dataMatrix.transpose()* error # 每步 weights 該變量 30 return weights.getA() 31 32 # 對數據分類的邊界圖形顯示 33 def plotBestFit(weights): 34 import matplotlib.pyplot as plt 35 dataMat,labelMat=loadDataSet() 36 dataArr = array(dataMat) 37 n = shape(dataArr)[0] # 數據的行數,即對象的個數 38 xcord1 = []; ycord1 = [] # 對類別號為 1 的對象,分 X 軸和 Y 軸的數據 39 xcord2 = []; ycord2 = [] # 對類別號為 0 的對象,分 X 軸和 Y 軸的數據 40 for i in range(n): # 對所有的對象進行遍歷 41 if int(labelMat[i])== 1: # 對象的類別為:1 42 xcord1.append(dataArr[i,1]); ycord1.append(dataArr[i,2]) 43 else: 44 xcord2.append(dataArr[i,1]); ycord2.append(dataArr[i,2]) 45 fig = plt.figure() 46 ax = fig.add_subplot(111) 47 ax.scatter(xcord1, ycord1, s=30, c=‘red‘, marker=‘s‘) # 對散點的格式的設置,坐標號、點的大小、顏色、點的圖形(方塊) 48 ax.scatter(xcord2, ycord2, s=30, c=‘green‘) # 點的圖形默認為圓 49 x = arange(-3.0, 3.0, 0.1) 50 y = (-weights[0]-weights[1]*x)/weights[2] # 線性方程 y = aX + b,y 是數據第三列的特征,X 是數據第二列的特征 51 ax.plot(x, y) 52 plt.xlabel(‘X1‘); plt.ylabel(‘X2‘) 53 plt.show() 54 55 # 隨機梯度上升法 56 def stocGradAscent0(dataMatrix, classLabels, numIter=150): 57 m,n = shape(dataMatrix ) 58 weights = ones(n) 59 for j in range(numIter): # 叠代次數 60 dataIndex = range(m) # 61 for i in range(m): # 對所有對象的遍歷 62 alpha = 4/(1.0+j+i)+0.01 # 對步長的調整 63 randIndex = int (random.uniform(0,len(dataIndex))) # 隨機生成一個整數,介於 0 到 m 64 h = sigmoid(sum(dataMatrix[randIndex]*weights)) # 對隨機選擇的對象計算類別的數值(回歸系數值) 65 error = classLabels[randIndex] - h # 根據實際類型與計算類型值的誤差,損失函數 66 weights = weights + alpha * error * dataMatrix[randIndex] # 每步 weights 改變值 67 del(dataIndex [randIndex]) # 去除已經選擇過的對象,避免下次選中 68 return weights 69 70 # 對未知對象進行預測類別號 71 def classifyVector(inX , weights): 72 prob = sigmoid(sum(inX * weights)) # 計算回歸系數 73 if prob > 0.5: 74 return 1.0 75 else: 76 return 0.0 77 78 # 實例:從疝氣病預測病馬的死亡率 79 def colicTest(): 80 frTrain = open(‘horseColicTraining.txt‘); frTest = open(‘horseColicTest.txt‘)# 數據文件 81 trainingSet = []; trainingLabels = [] # 樣本集和類標號集的初始化 82 for line in frTrain.readlines(): 83 currLine = line.strip().split(‘\t‘) # 根據制表符進行字符串的分割 84 lineArr =[] 85 for i in range(21): 86 lineArr.append(float(currLine[i])) 87 trainingSet.append(lineArr) 88 trainingLabels.append(float(currLine[21])) 89 trainWeights = stocGradAscent0(array(trainingSet), trainingLabels, 500) # 通過隨機梯度上升法計算回歸系數 90 #trainWeights = gradAscent(array(trainingSet), trainingLabels) # 通過批量梯度上升法計算出回歸系數 91 errorCount = 0; numTestVec = 0.0 # 錯誤個數和錯誤率的初始化 92 # 預測樣本的遍歷 93 for line in frTest.readlines(): 94 numTestVec += 1.0 95 currLine = line.strip().split(‘\t‘) 96 lineArr =[] 97 for i in range(21): 98 lineArr.append(float(currLine[i])) 99 if int(classifyVector(array(lineArr), trainWeights))!= int(currLine[21]): # 預測的類別號與真實類別號的比較 100 errorCount += 1 101 errorRate = (float(errorCount)/numTestVec) 102 print u‘本次測試的錯誤率是; %f‘ % errorRate 103 return errorRate 104 105 # 預測結果函數 106 def multiTest(): 107 numTests = 10; errorSum = 0.0 108 for k in range(numTests): 109 errorSum += colicTest() 110 print u‘經過 %d 次測試結果的平均錯誤率是: %f ‘ % (numTests ,errorSum /float(numTests)) 111 112 if __name__ == ‘__main__‘: 113 # 用批量梯度上升法畫邊界線 114 ‘‘‘ 115 dataSet, labelSet = loadDataSet() 116 #weights = gradAscent(dataSet, labelSet) 117 #plotBestFit(weights) 118 ‘‘‘ 119 # 用隨機梯度上升法畫邊界線 120 ‘‘‘ 121 dataSet, labelSet = loadDataSet() 122 #weightss = stocGradAscent0(array(dataSet), labelSet ) 123 #plotBestFit(weightss) 124 ‘‘‘ 125 # 實例的預測 126 multiTest()View Code
機器學習之Logistic 回歸算法