第5章 Logistic回歸
Logistic 回歸 概述
Logistic 回歸雖然名字叫回歸，但是它是用來做分類的。其主要思想是: 根據現有數據對分類邊界線建立回歸公式，以此進行分類。
須知概念
Sigmoid 函數
回歸 概念
假設現在有一些數據點，我們用一條直線對這些點進行擬合（這條直線稱為最佳擬合直線），這個擬合的過程就叫做回歸。進而可以得到對這些點的擬合直線方程，那么我們根據這個回歸方程，怎么進行分類呢？請看下面。
二值型輸出分類函數
我們想要的函數應該是: 能接受所有的輸入然后預測出類別。例如，在兩個類的情況下，上述函數輸出 0 或 1.或許你之前接觸過具有這種性質的函數，該函數稱為 海維塞得階躍函數(Heaviside step function)，或者直接稱為 單位階躍函數。然而，海維塞得階躍函數的問題在於: 該函數在跳躍點上從 0 瞬間跳躍到 1，這個瞬間跳躍過程有時很難處理。幸好，另一個函數也有類似的性質（可以輸出 0 或者 1 的性質），且數學上更易處理，這就是 Sigmoid 函數。 Sigmoid 函數具體的計算公式如下:
下圖給出了 Sigmoid 函數在不同坐標尺度下的兩條曲線圖。當 x 為 0 時，Sigmoid 函數值為 0.5 。隨着 x 的增大，對應的 Sigmoid 值將逼近於 1 ; 而隨着 x 的減小， Sigmoid 值將逼近於 0 。如果橫坐標刻度足夠大， Sigmoid 函數看起來很像一個階躍函數。
因此，為了實現 Logistic 回歸分類器，我們可以在每個特征上都乘以一個回歸系數（如下公式所示），然后把所有結果值相加，將這個總和代入 Sigmoid 函數中，進而得到一個范圍在 0~1 之間的數值。任何大於 0.5 的數據被分入 1 類，小於 0.5 即被歸入 0 類。所以， Logistic 回歸也可以被看成是一種概率估計。
基於最優化方法的回歸系數確定
Sigmoid 函數的輸入記為 z ，由下面公式得到:
如果采用向量的寫法，上述公式可以寫成 
，它表示將這兩個數值向量對應元素相乘然后全部加起來即得到 z 值。其中的向量 x 是分類器的輸入數據，向量 w 也就是我們要找到的最佳參數（系數），從而使得分類器盡可能地精確。為了尋找該最佳參數，需要用到最優化理論的一些知識。我們這里使用的是——梯度上升法。
梯度上升法
梯度的介紹
需要一點點向量方面的數學基礎
向量 = 值 + 方向
梯度 = 向量
梯度 = 梯度值 + 梯度方向
梯度上升法的思想
要找到某函數的最大值，最好的方法是沿着該函數的梯度方向探尋。如果梯度記為 ▽ ，則函數 f(x, y) 的梯度由下式表示:
這個梯度意味着要沿 x 的方向移動 
，沿 y 的方向移動 
。其中，函數f(x, y) 必須要在待計算的點上有定義並且可微。下圖是一個具體的例子。
上圖展示的，梯度上升算法到達每個點后都會重新估計移動的方向。從 P0 開始，計算完該點的梯度，函數就根據梯度移動到下一點 P1。在 P1 點，梯度再次被重新計算，並沿着新的梯度方向移動到 P2 。如此循環迭代，直到滿足停止條件。迭代過程中，梯度算子總是保證我們能選取到最佳的移動方向。
上圖中的梯度上升算法沿梯度方向移動了一步。可以看到，梯度算子總是指向函數值增長最快的方向。這里所說的是移動方向，而未提到移動量的大小。該量值稱為步長，記作 α 。用向量來表示的話，梯度上升算法的迭代公式如下:
該公式將一直被迭代執行，直至達到某個停止條件為止，比如迭代次數達到某個指定值或者算法達到某個可以允許的誤差范圍。
介紹一下幾個相關的概念：
例如：y = w1x1 + w2x2 + ... + wnxn
梯度：參考上圖的例子，二維圖像，x方向代表第一個系數，也就是 w1，y方向代表第二個系數也就是 w2，這樣的向量就是梯度。
α：上面的梯度算法的迭代公式中的阿爾法，這個代表的是移動步長。移動步長會影響最終結果的擬合程度，最好的方法就是隨着迭代次數更改移動步長。
步長通俗的理解，100米，如果我一步走10米，我需要走10步；如果一步走20米，我只需要走5步。這里的一步走多少米就是步長的意思。
▽f(w)：代表沿着梯度變化的方向。
Note: 我們常聽到的是梯度下降算法，它與這里的梯度上升算法是一樣的，只是公式中的加法需要變成減法。因此，對應的公式可以寫成
梯度上升算法用來求函數的最大值，而梯度下降算法用來求函數的最小值。
局部最優現象
上圖表示參數 θ 與誤差函數 J(θ) 的關系圖，紅色的部分是表示 J(θ) 有着比較高的取值，我們需要的是，能夠讓 J(θ) 的值盡量的低。也就是深藍色的部分。θ0，θ1 表示 θ 向量的兩個維度。
可能梯度下降的最終點並非是全局最小點，可能是一個局部最小點，如我們上圖中的右邊的梯度下降曲線，描述的是最終到達一個局部最小點，這是我們重新選擇了一個初始點得到的。
看來我們這個算法將會在很大的程度上被初始點的選擇影響而陷入局部最小點。
Logistic 回歸 原理
Logistic 回歸 工作原理
每個回歸系數初始化為 1
重復 R 次:
計算整個數據集的梯度
使用 步長 x 梯度 更新回歸系數的向量
返回回歸系數
Logistic 回歸 開發流程
收集數據: 采用任意方法收集數據
准備數據: 由於需要進行距離計算，因此要求數據類型為數值型。另外，結構化數據格式則最佳。
分析數據: 采用任意方法對數據進行分析。
訓練算法: 大部分時間將用於訓練，訓練的目的是為了找到最佳的分類回歸系數。
測試算法: 一旦訓練步驟完成，分類將會很快。
使用算法: 首先，我們需要輸入一些數據，並將其轉換成對應的結構化數值；接着，基於訓練好的回歸系數就可以對這些數值進行簡單的回歸計算，判定它們屬於哪個類別；在這之后，我們就可以在輸出的類別上做一些其他分析工作。
Logistic 回歸 算法特點
優點: 計算代價不高，易於理解和實現。
缺點: 容易欠擬合，分類精度可能不高。
適用數據類型: 數值型和標稱型數據。
附加 方向導數與梯度
Logistic 回歸 項目案例
項目案例1: 使用 Logistic 回歸在簡單數據集上的分類
項目概述
在一個簡單的數據集上，采用梯度上升法找到 Logistic 回歸分類器在此數據集上的最佳回歸系數
開發流程
收集數據: 可以使用任何方法
准備數據: 由於需要進行距離計算，因此要求數據類型為數值型。另外，結構化數據格式則最佳
分析數據: 畫出決策邊界
訓練算法: 使用梯度上升找到最佳參數
測試算法: 使用 Logistic 回歸進行分類
使用算法: 對簡單數據集中數據進行分類
收集數據: 可以使用任何方法
我們采用存儲在 TestSet.txt 文本文件中的數據，存儲格式如下:
-0.017612 14.053064 0
-1.395634 4.662541 1
-0.752157 6.538620 0
-1.322371 7.152853 0
0.423363 11.054677 0
繪制在圖中，如下圖所示:
准備數據: 由於需要進行距離計算，因此要求數據類型為數值型。另外，結構化數據格式則最佳
分析數據: 畫出決策邊界
畫出數據集和 Logistic 回歸最佳擬合直線的函數
def plotBestFit(dataArr, labelMat, weights): ''' Desc: 將我們得到的數據可視化展示出來 Args: dataArr:樣本數據的特征 labelMat:樣本數據的類別標簽，即目標變量 weights:回歸系數 Returns: None ''' 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, c='red', marker='s') ax.scatter(xcord2, ycord2, s=30, c='green') x = arange(-3.0, 3.0, 0.1) """ y的由來，卧槽，是不是沒看懂？ 首先理論上是這個樣子的。 dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])]) w0*x0+w1*x1+w2*x2=f(x) x0最開始就設置為1叻， x2就是我們畫圖的y值，而f(x)被我們磨合誤差給算到w0,w1,w2身上去了 所以： w0+w1*x+w2*y=0 => y = (-w0-w1*x)/w2 """ y = (-weights[0]-weights[1]*x)/weights[2] ax.plot(x, y) plt.xlabel('X'); plt.ylabel('Y') plt.show()
訓練算法: 使用梯度上升找到最佳參數
Logistic 回歸梯度上升優化算法
# 正常的處理方案 # 兩個參數：第一個參數==> dataMatIn 是一個2維NumPy數組，每列分別代表每個不同的特征，每行則代表每個訓練樣本。 # 第二個參數==> classLabels 是類別標簽，它是一個 1*100 的行向量。為了便於矩陣計算，需要將該行向量轉換為列向量，做法是將原向量轉置，再將它賦值給labelMat。 def gradAscent(dataMatIn, classLabels): # 轉化為矩陣[[1,1,2],[1,1,2]....] dataMatrix = mat(dataMatIn) # 轉換為 NumPy 矩陣 # 轉化為矩陣[[0,1,0,1,0,1.....]]，並轉制[[0],[1],[0].....] # transpose() 行列轉置函數 # 將行向量轉化為列向量 => 矩陣的轉置 labelMat = mat(classLabels).transpose() # 首先將數組轉換為 NumPy 矩陣，然后再將行向量轉置為列向量 # m->數據量，樣本數 n->特征數 m,n = shape(dataMatrix) # print m, n, '__'*10, shape(dataMatrix.transpose()), '__'*100 # alpha代表向目標移動的步長 alpha = 0.001 # 迭代次數 maxCycles = 500 # 生成一個長度和特征數相同的矩陣，此處n為3 -> [[1],[1],[1]] # weights 代表回歸系數， 此處的 ones((n,1)) 創建一個長度和特征數相同的矩陣，其中的數全部都是 1 weights = ones((n,1)) for k in range(maxCycles): #heavy on matrix operations # m*3 的矩陣 * 3*1 的單位矩陣 ＝ m*1的矩陣 # 那么乘上單位矩陣的意義，就代表：通過公式得到的理論值 # 參考地址： 矩陣乘法的本質是什么？ https://www.zhihu.com/question/21351965/answer/31050145 # print 'dataMatrix====', dataMatrix # print 'weights====', weights # n*3 * 3*1 = n*1 h = sigmoid(dataMatrix*weights) # 矩陣乘法 # print 'hhhhhhh====', h # labelMat是實際值 error = (labelMat - h) # 向量相減 # 0.001* (3*m)*(m*1) 表示在每一個列上的一個誤差情況，最后得出 x1,x2,xn的系數的偏移量 weights = weights + alpha * dataMatrix.transpose() * error # 矩陣乘法，最后得到回歸系數 return array(weights)
大家看到這兒可能會有一些疑惑，就是，我們在迭代中更新我們的回歸系數，后邊的部分是怎么計算出來的？為什么會是 alpha * dataMatrix.transpose() * error ?因為這就是我們所求的梯度，也就是對 f(w) 對 w 求一階導數。具體推導如下:
測試算法: 使用 Logistic 回歸進行分類
def testLR(): # 1.收集並准備數據 dataMat, labelMat = loadDataSet("input/5.Logistic/TestSet.txt") # print dataMat, '---\n', labelMat # 2.訓練模型， f(x)=a1*x1+b2*x2+..+nn*xn中 (a1,b2, .., nn).T的矩陣值 # 因為數組沒有是復制n份， array的乘法就是乘法 dataArr = array(dataMat) # print dataArr weights = gradAscent(dataArr, labelMat) # weights = stocGradAscent0(dataArr, labelMat) # weights = stocGradAscent1(dataArr, labelMat) # print '*'*30, weights # 數據可視化 plotBestFit(dataArr, labelMat, weights)
使用算法: 對簡單數據集中數據進行分類
注意
梯度上升算法在每次更新回歸系數時都需要遍歷整個數據集，該方法在處理 100 個左右的數據集時尚可，但如果有數十億樣本和成千上萬的特征，那么該方法的計算復雜度就太高了。一種改進方法是一次僅用一個樣本點來更新回歸系數，該方法稱為 隨機梯度上升算法。由於可以在新樣本到來時對分類器進行增量式更新，因而隨機梯度上升算法是一個在線學習算法。與 “在線學習” 相對應，一次處理所有數據被稱作是 “批處理”。
隨機梯度上升算法可以寫成如下的偽代碼:
所有回歸系數初始化為 1
對數據集中每個樣本
計算該樣本的梯度
使用 alpha x gradient 更新回歸系數值
返回回歸系數值
以下是隨機梯度上升算法的實現代碼:
# 隨機梯度上升 # 梯度上升優化算法在每次更新數據集時都需要遍歷整個數據集，計算復雜都較高 # 隨機梯度上升一次只用一個樣本點來更新回歸系數 def stocGradAscent0(dataMatrix, classLabels): m,n = shape(dataMatrix) alpha = 0.01 # n*1的矩陣 # 函數ones創建一個全1的數組 weights = ones(n) # 初始化長度為n的數組，元素全部為 1 for i in range(m): # sum(dataMatrix[i]*weights)為了求 f(x)的值， f(x)=a1*x1+b2*x2+..+nn*xn,此處求出的 h 是一個具體的數值，而不是一個矩陣 h = sigmoid(sum(dataMatrix[i]*weights)) # print 'dataMatrix[i]===', dataMatrix[i] # 計算真實類別與預測類別之間的差值，然后按照該差值調整回歸系數 error = classLabels[i] - h # 0.01*(1*1)*(1*n) print weights, "*"*10 , dataMatrix[i], "*"*10 , error weights = weights + alpha * error * dataMatrix[i] return weights
可以看到，隨機梯度上升算法與梯度上升算法在代碼上很相似，但也有一些區別: 第一，后者的變量 h 和誤差 error 都是向量，而前者則全是數值；第二，前者沒有矩陣的轉換過程，所有變量的數據類型都是 NumPy 數組。
判斷優化算法優劣的可靠方法是看它是否收斂，也就是說參數是否達到了穩定值，是否還會不斷地變化？下圖展示了隨機梯度上升算法在 200 次迭代過程中回歸系數的變化情況。其中的系數2，也就是 X2 只經過了 50 次迭代就達到了穩定值，但系數 1 和 0 則需要更多次的迭代。如下圖所示:
針對這個問題，我們改進了之前的隨機梯度上升算法，如下:
# 隨機梯度上升算法（隨機化） def stocGradAscent1(dataMatrix, classLabels, numIter=150): m,n = shape(dataMatrix) weights = ones(n) # 創建與列數相同的矩陣的系數矩陣，所有的元素都是1 # 隨機梯度, 循環150,觀察是否收斂 for j in range(numIter): # [0, 1, 2 .. m-1] dataIndex = range(m) for i in range(m): # i和j的不斷增大，導致alpha的值不斷減少，但是不為0 alpha = 4/(1.0+j+i)+0.0001 # alpha 會隨着迭代不斷減小，但永遠不會減小到0，因為后邊還有一個常數項0.0001 # 隨機產生一個 0～len()之間的一個值 # random.uniform(x, y) 方法將隨機生成下一個實數，它在[x,y]范圍內,x是這個范圍內的最小值，y是這個范圍內的最大值。 randIndex = int(random.uniform(0,len(dataIndex))) # sum(dataMatrix[i]*weights)為了求 f(x)的值， f(x)=a1*x1+b2*x2+..+nn*xn h = sigmoid(sum(dataMatrix[randIndex]*weights)) error = classLabels[randIndex] - h # print weights, '__h=%s' % h, '__'*20, alpha, '__'*20, error, '__'*20, dataMatrix[randIndex] weights = weights + alpha * error * dataMatrix[randIndex] del(dataIndex[randIndex]) return weights
上面的改進版隨機梯度上升算法，我們修改了兩處代碼。
第一處改進為 alpha 的值。alpha 在每次迭代的時候都會調整，這回緩解上面波動圖的數據波動或者高頻波動。另外，雖然 alpha 會隨着迭代次數不斷減少，但永遠不會減小到 0，因為我們在計算公式中添加了一個常數項。
第二處修改為 randIndex 更新，這里通過隨機選取樣本拉來更新回歸系數。這種方法將減少周期性的波動。這種方法每次隨機從列表中選出一個值，然后從列表中刪掉該值（再進行下一次迭代）。
程序運行之后能看到類似於下圖的結果圖。
完整代碼地址: https://github.com/apachecn/MachineLearning/blob/master/src/python/5.Logistic/logistic.py
項目案例2: 從疝氣病症預測病馬的死亡率
項目概述
使用 Logistic 回歸來預測患有疝病的馬的存活問題。疝病是描述馬胃腸痛的術語。然而，這種病不一定源自馬的胃腸問題，其他問題也可能引發馬疝病。這個數據集中包含了醫院檢測馬疝病的一些指標，有的指標比較主觀，有的指標難以測量，例如馬的疼痛級別。
開發流程
收集數據: 給定數據文件
准備數據: 用 Python 解析文本文件並填充缺失值
分析數據: 可視化並觀察數據
訓練算法: 使用優化算法，找到最佳的系數
測試算法: 為了量化回歸的效果，需要觀察錯誤率。根據錯誤率決定是否回退到訓練階段，
通過改變迭代的次數和步長的參數來得到更好的回歸系數
使用算法: 實現一個簡單的命令行程序來手機馬的症狀並輸出預測結果並非難事，
這可以作為留給大家的一道習題
收集數據: 給定數據文件
病馬的訓練數據已經給出來了，如下形式存儲在文本文件中:
1.000000 1.000000 39.200000 88.000000 20.000000 0.000000 0.000000 4.000000 1.000000 3.000000 4.000000 2.000000 0.000000 0.000000 0.000000 4.000000 2.000000 50.000000 85.000000 2.000000 2.000000 0.000000
2.000000 1.000000 38.300000 40.000000 24.000000 1.000000 1.000000 3.000000 1.000000 3.000000 3.000000 1.000000 0.000000 0.000000 0.000000 1.000000 1.000000 33.000000 6.700000 0.000000 0.000000 1.000000
准備數據: 用 Python 解析文本文件並填充缺失值
處理數據中的缺失值
假設有100個樣本和20個特征，這些數據都是機器收集回來的。若機器上的某個傳感器損壞導致一個特征無效時該怎么辦？此時是否要扔掉整個數據？這種情況下，另外19個特征怎么辦？ 它們是否還可以用？答案是肯定的。因為有時候數據相當昂貴，扔掉和重新獲取都是不可取的，所以必須采用一些方法來解決這個問題。
下面給出了一些可選的做法：
- 使用可用特征的均值來填補缺失值；
- 使用特殊值來填補缺失值，如 -1；
- 忽略有缺失值的樣本；
- 使用有相似樣本的均值添補缺失值；
- 使用另外的機器學習算法預測缺失值。
現在，我們對下一節要用的數據集進行預處理，使其可以順利地使用分類算法。在預處理需要做兩件事:
-
所有的缺失值必須用一個實數值來替換，因為我們使用的 NumPy 數據類型不允許包含缺失值。我們這里選擇實數 0 來替換所有缺失值，恰好能適用於 Logistic 回歸。這樣做的直覺在於，我們需要的是一個在更新時不會影響系數的值。回歸系數的更新公式如下:
weights = weights + alpha * error * dataMatrix[randIndex]
如果 dataMatrix 的某個特征對應值為 0，那么該特征的系數將不做更新，即:
weights = weights
另外，由於 Sigmoid(0) = 0.5 ，即它對結果的預測不具有任何傾向性，因此我們上述做法也不會對誤差造成任何影響。基於上述原因，將缺失值用 0 代替既可以保留現有數據，也不需要對優化算法進行修改。此外，該數據集中的特征取值一般不為 0，因此在某種意義上說它也滿足 “特殊值” 這個要求。
-
如果在測試數據集中發現了一條數據的類別標簽已經缺失，那么我們的簡單做法是將該條數據丟棄。這是因為類別標簽與特征不同，很難確定采用某個合適的值來替換。采用 Logistic 回歸進行分類時這種做法是合理的，而如果采用類似 kNN 的方法就可能不太可行。
原始的數據集經過預處理后，保存成兩個文件: horseColicTest.txt 和 horseColicTraining.txt 。
分析數據: 可視化並觀察數據
將數據使用 MatPlotlib 打印出來，觀察數據是否是我們想要的格式
訓練算法: 使用優化算法，找到最佳的系數
下面給出 原始的梯度上升算法，隨機梯度上升算法，改進版隨機梯度上升算法 的代碼:
# 正常的處理方案 # 兩個參數：第一個參數==> dataMatIn 是一個2維NumPy數組，每列分別代表每個不同的特征，每行則代表每個訓練樣本。 # 第二個參數==> classLabels 是類別標簽，它是一個 1*100 的行向量。為了便於矩陣計算，需要將該行向量轉換為列向量，做法是將原向量轉置，再將它賦值給labelMat。 def gradAscent(dataMatIn, classLabels): # 轉化為矩陣[[1,1,2],[1,1,2]....] dataMatrix = mat(dataMatIn) # 轉換為 NumPy 矩陣 # 轉化為矩陣[[0,1,0,1,0,1.....]]，並轉制[[0],[1],[0].....] # transpose() 行列轉置函數 # 將行向量轉化為列向量 => 矩陣的轉置 labelMat = mat(classLabels).transpose() # 首先將數組轉換為 NumPy 矩陣，然后再將行向量轉置為列向量 # m->數據量，樣本數 n->特征數 m,n = shape(dataMatrix) # print m, n, '__'*10, shape(dataMatrix.transpose()), '__'*100 # alpha代表向目標移動的步長 alpha = 0.001 # 迭代次數 maxCycles = 500 # 生成一個長度和特征數相同的矩陣，此處n為3 -> [[1],[1],[1]] # weights 代表回歸系數， 此處的 ones((n,1)) 創建一個長度和特征數相同的矩陣，其中的數全部都是 1 weights = ones((n,1)) for k in range(maxCycles): #heavy on matrix operations # m*3 的矩陣 * 3*1 的單位矩陣 ＝ m*1的矩陣 # 那么乘上單位矩陣的意義，就代表：通過公式得到的理論值 # 參考地址： 矩陣乘法的本質是什么？ https://www.zhihu.com/question/21351965/answer/31050145 # print 'dataMatrix====', dataMatrix # print 'weights====', weights # n*3 * 3*1 = n*1 h = sigmoid(dataMatrix*weights) # 矩陣乘法 # print 'hhhhhhh====', h # labelMat是實際值 error = (labelMat - h) # 向量相減 # 0.001* (3*m)*(m*1) 表示在每一個列上的一個誤差情況，最后得出 x1,x2,xn的系數的偏移量 weights = weights + alpha * dataMatrix.transpose() * error # 矩陣乘法，最后得到回歸系數 return array(weights) # 隨機梯度上升 # 梯度上升優化算法在每次更新數據集時都需要遍歷整個數據集，計算復雜都較高 # 隨機梯度上升一次只用一個樣本點來更新回歸系數 def stocGradAscent0(dataMatrix, classLabels): m,n = shape(dataMatrix) alpha = 0.01 # n*1的矩陣 # 函數ones創建一個全1的數組 weights = ones(n) # 初始化長度為n的數組，元素全部為 1 for i in range(m): # sum(dataMatrix[i]*weights)為了求 f(x)的值， f(x)=a1*x1+b2*x2+..+nn*xn,此處求出的 h 是一個具體的數值，而不是一個矩陣 h = sigmoid(sum(dataMatrix[i]*weights)) # print 'dataMatrix[i]===', dataMatrix[i] # 計算真實類別與預測類別之間的差值，然后按照該差值調整回歸系數 error = classLabels[i] - h # 0.01*(1*1)*(1*n) print weights, "*"*10 , dataMatrix[i], "*"*10 , error weights = weights + alpha * error * dataMatrix[i] return weights # 隨機梯度上升算法（隨機化） def stocGradAscent1(dataMatrix, classLabels, numIter=150): m,n = shape(dataMatrix) weights = ones(n) # 創建與列數相同的矩陣的系數矩陣，所有的元素都是1 # 隨機梯度, 循環150,觀察是否收斂 for j in range(numIter): # [0, 1, 2 .. m-1] dataIndex = range(m) for i in range(m): # i和j的不斷增大，導致alpha的值不斷減少，但是不為0 alpha = 4/(1.0+j+i)+0.0001 # alpha 會隨着迭代不斷減小，但永遠不會減小到0，因為后邊還有一個常數項0.0001 # 隨機產生一個 0～len()之間的一個值 # random.uniform(x, y) 方法將隨機生成下一個實數，它在[x,y]范圍內,x是這個范圍內的最小值，y是這個范圍內的最大值。 randIndex = int(random.uniform(0,len(dataIndex))) # sum(dataMatrix[i]*weights)為了求 f(x)的值， f(x)=a1*x1+b2*x2+..+nn*xn h = sigmoid(sum(dataMatrix[randIndex]*weights)) error = classLabels[randIndex] - h # print weights, '__h=%s' % h, '__'*20, alpha, '__'*20, error, '__'*20, dataMatrix[randIndex] weights = weights + alpha * error * dataMatrix[randIndex] del(dataIndex[randIndex]) return weights
測試算法: 為了量化回歸的效果，需要觀察錯誤率。根據錯誤率決定是否回退到訓練階段，通過改變迭代的次數和步長的參數來得到更好的回歸系數
Logistic 回歸分類函數
# 分類函數，根據回歸系數和特征向量來計算 Sigmoid的值 def classifyVector(inX, weights): ''' Desc: 最終的分類函數，根據回歸系數和特征向量來計算 Sigmoid 的值，大於0.5函數返回1，否則返回0 Args: inX -- 特征向量，features weights -- 根據梯度下降/隨機梯度下降 計算得到的回歸系數 Returns: 如果 prob 計算大於 0.5 函數返回 1 否則返回 0 ''' prob = sigmoid(sum(inX*weights)) if prob > 0.5: return 1.0 else: return 0.0 # 打開測試集和訓練集,並對數據進行格式化處理 def colicTest(): ''' Desc: 打開測試集和訓練集，並對數據進行格式化處理 Args: None Returns: errorRate -- 分類錯誤率 ''' frTrain = open('input/5.Logistic/horseColicTraining.txt') frTest = open('input/5.Logistic/horseColicTest.txt') trainingSet = [] trainingLabels = [] # 解析訓練數據集中的數據特征和Labels # trainingSet 中存儲訓練數據集的特征，trainingLabels 存儲訓練數據集的樣本對應的分類標簽 for line in frTrain.readlines(): currLine = line.strip().split('\t') lineArr = [] for i in range(21): lineArr.append(float(currLine[i])) trainingSet.append(lineArr) trainingLabels.append(float(currLine[21])) # 使用 改進后的 隨機梯度下降算法 求得在此數據集上的最佳回歸系數 trainWeights trainWeights = stocGradAscent1(array(trainingSet), trainingLabels, 500) errorCount = 0 numTestVec = 0.0 # 讀取 測試數據集 進行測試，計算分類錯誤的樣本條數和最終的錯誤率 for line in frTest.readlines(): numTestVec += 1.0 currLine = line.strip().split('\t') lineArr = [] for i in range(21): lineArr.append(float(currLine[i])) if int(classifyVector(array(lineArr), trainWeights)) != int(currLine[21]): errorCount += 1 errorRate = (float(errorCount) / numTestVec) print "the error rate of this test is: %f" % errorRate return errorRate # 調用 colicTest() 10次並求結果的平均值 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, errorSum/float(numTests))
使用算法: 實現一個簡單的命令行程序來手機馬的症狀並輸出預測結果並非難事，這可以作為留給大家的一道習題
完整代碼地址: https://github.com/apachecn/MachineLearning/blob/master/src/python/5.Logistic/logistic.py
- 作者：羊三 小瑤
- GitHub地址: https://github.com/apachecn/MachineLearning
- 版權聲明：歡迎轉載學習 => 請標注信息來源於 ApacheCN