簡單 生成 選擇 效率 split max 坐標 opened 似然函數

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

機器學習之Logistic 回歸算法