機器學習之Logistic 迴歸算法

聲明：本篇博文是學習《機器學習實戰》一書的方式路程，系原創，若轉載請標明來源。

1 Logistic 迴歸算法的原理

1.1 須要的數學基礎

我在看機器學習實戰時對其中的代碼很是費解，說好的利用偏導數求最值怎麼代碼中沒有體現啊，就一個簡單的式子：θ= θ - α Σ [( h_θ(x⁽ⁱ⁾)-y⁽ⁱ⁾ ) ] * x_{i 。}通過查找資料才知道，書中省去了大量的理論推導過程，其中用到了線性函數、sigmoid 函數、偏導數、最大似然函數、梯度降低法。下面讓咱們一窺究竟，是站在大神的肩膀描述我本身的看法。

1.2 Logistic 迴歸的引入

Logistic 迴歸是機率非線性模型，用於處理二元分類結果，名爲迴歸實爲分類。下面給出一個二元分類的例子，以下圖所示：

圖中數據分紅兩類，打叉的一類用 y = 1 表示，另外一類爲圓圈用 y= 0 表示。如今咱們要對數據進行分類，要找到一個分類函數(邊界線)，咱們很容易得出一個線性函數對其分類：0 = θ₀ + θ₁x₁ + θ₂x₂ 。但咱們想要的函數應該是，能接受全部的輸入而後預測類別。例如：在兩個分類的狀況下，函數輸出 0 或 1。所以，咱們就須要引入Sigmoid 函數，這個函數的性質能夠知足要求。Sigmoid 函數：

Sigmoid 函數的值域爲（0,1），且具備良好的從0 到 1 的跳躍性，如在兩個不一樣的座標尺度下的函數圖形：

因此，咱們把線性方程和Sigmoid 函數結合起來就能解決問題。即 ：分類預測函數 h_θ (x) = g( θ₀ + θ₁x₁ + θ₂x₂) .咱們就能夠對樣本數據進行分類，以下圖所示：

 

對於線性的分類邊界，以下形式：

分類預測函數，以下形式：

其中，θ^T 是向量 θ (θ₀, θ₁,... ,θ_n) 的轉置，向量 x ( x₀ , x₁ ,... , x_n)，n -1爲數據的維度，x₀ =1，這是便於計算。

 1.3 分類預測函數問題的轉化成求θ

 經過上面的分析，咱們得出了分類預測函數 h_θ(x) , 但其中向量 x 是已知的（x 是未知類別號的對象數據），向量 θ 未知，即咱們把求分類函數問題轉化成求向量 θ 。由於Sigmoid 函數的取值區間（0,1），那咱們能夠看作機率 P（y = 1 | x_i ; θ）= h_θ(x) , 表示在 x_i 肯定的狀況下，類別 y = 1 的機率。由此，咱們也能夠得出在 x_i 肯定的狀況下，類別 y = 0 的機率  P（y = 0 | x_i ; θ）= 1 -  P（y = 1 | x_i ; θ）= 1 - h_θ(x) . 即 ：

咱們能夠將這兩個式子合併得：

其中的 y = 0 或 1 .

這時候咱們能夠利用最大似然函數對向量 θ 求值，能夠理解爲選取的樣本數據的機率是最大的，那麼樣本數爲 m 的似然函數爲：

經過對數的形式對似然函數進行變化，對數似然函數：

 

 這裏的最大似然函數的值就是咱們所要求的向量 θ , 而求解的方法利用梯度降低法。

1.4 梯度降低法求解θ

 在用梯度降低法時，咱們將會利用Sigmoid 函數的一個性質： g^, (z) = g(z）[ 1- g(z) ] 。

構造一個Cost函數（損失函數），該函數表示預測的輸出（h）與訓練數據類別（y）之間的誤差，能夠是兩者之間的差（h-y）或者是其餘的形式。綜合考慮全部訓練數據的「損失」，將Cost求和或者求平均，記爲J(θ)函數，表示全部訓練數據預測值與實際類別的誤差。

損失函數：

 

J(θ)代價函數：

 

其中，x⁽ⁱ⁾ 每一個樣本數據點在某一個特徵上的值，即特徵向量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