Python3《機器學習實戰》學習筆記（六）：Logistic迴歸基礎篇之梯度上升演算法
一 前言
本文從Logistic迴歸的原理開始講起，補充了書上省略的數學推導。本文可能會略顯枯燥，理論居多，Sklearn實戰內容會放在下一篇文章。自己慢慢推導完公式，還是蠻開心的一件事。
二 Logistic迴歸與梯度上升演算法
Logistic迴歸是眾多分類演算法中的一員。通常，Logistic迴歸用於二分類問題，例如預測明天是否會下雨。當然它也可以用於多分類問題，不過為了簡單起見，本文暫先討論二分類問題。首先，讓我們來了解一下，什麼是Logistic迴歸。
1 Logistic迴歸
假設現在有一些資料點，我們利用一條直線對這些點進行擬合(該線稱為最佳擬合直線)，這個擬合過程就稱作為迴歸，如下圖所示：
Logistic迴歸一種二分類演算法，它利用的是Sigmoid函式閾值在[0,1]這個特性。Logistic迴歸進行分類的主要思想是：根據現有資料對分類邊界線建立迴歸公式，以此進行分類。其實，Logistic本質上是一個基於條件概率的判別模型(Discriminative Model)。
所以要想了解Logistic迴歸，我們必須先看一看Sigmoid函式 ，我們也可以稱它為Logistic函式。它的公式如下：
整合成一個公式，就變成了如下公式：
下面這張圖片，為我們展示了Sigmoid函式的樣子。
z是一個矩陣，θ是引數列向量(要求解的)，x是樣本列向量(給定的資料集)。θ^T表示θ的轉置。g(z)函式實現了任意實數到[0,1]的對映，這樣我們的資料集([x0,x1,…,xn])，不管是大於1或者小於0，都可以對映到[0,1]區間進行分類。hθ(x)給出了輸出為1的概率。比如當hθ(x)=0.7，那麼說明有70%的概率輸出為1。輸出為0的概率是輸出為1的補集，也就是30%。
如果我們有合適的引數列向量θ([θ0,θ1,…θn]^T)，以及樣本列向量x([x0,x1,…,xn])，那麼我們對樣本x分類就可以通過上述公式計算出一個概率，如果這個概率大於0.5，我們就可以說樣本是正樣本，否則樣本是負樣本。
舉個例子，對於"垃圾郵件判別問題"，對於給定的郵件(樣本)，我們定義非垃圾郵件為正類，垃圾郵件為負類。我們通過計算出的概率值即可判定郵件是否是垃圾郵件。
那麼問題來了！如何得到合適的引數向量θ?
根據sigmoid函式的特性，我們可以做出如下的假設：
上式即為在已知樣本x和引數θ的情況下，樣本x屬性正樣本(y=1)和負樣本(y=0)的條件概率。理想狀態下，根據上述公式，求出各個點的概率均為1，也就是完全分類都正確。但是考慮到實際情況，樣本點的概率越接近於1，其分類效果越好。比如一個樣本屬於正樣本的概率為0.51，那麼我們就可以說明這個樣本屬於正樣本。另一個樣本屬於正樣本的概率為0.99，那麼我們也可以說明這個樣本屬於正樣本。但是顯然，第二個樣本概率更高，更具說服力。我們可以把上述兩個概率公式合二為一：
合併出來的Cost，我們稱之為代價函式(Cost Function)。當y等於1時，(1-y)項(第二項)為0；當y等於0時，y項(第一項)為0。為了簡化問題，我們對整個表示式求對數，(將指數問題對數化是處理數學問題常見的方法)：
這個代價函式，是對於一個樣本而言的。給定一個樣本，我們就可以通過這個代價函式求出，樣本所屬類別的概率，而這個概率越大越好，所以也就是求解這個代價函式的最大值。既然概率出來了，那麼最大似然估計也該出場了。假定樣本與樣本之間相互獨立，那麼整個樣本集生成的概率即為所有樣本生成概率的乘積，再將公式對數化，便可得到如下公式：
其中，m為樣本的總數，y(i)表示第i個樣本的類別，x(i)表示第i個樣本，需要注意的是θ是多維向量，x(i)也是多維向量。
綜上所述，滿足J(θ)的最大的θ值即是我們需要求解的模型。
怎麼求解使J(θ)最大的θ值呢？因為是求最大值，所以我們需要使用梯度上升演算法。如果面對的問題是求解使J(θ)最小的θ值，那麼我們就需要使用梯度下降演算法。面對我們這個問題，如果使J(θ) := -J(θ)，那麼問題就從求極大值轉換成求極小值了，使用的演算法就從梯度上升演算法變成了梯度下降演算法，它們的思想都是相同的，學會其一，就也會了另一個。本文使用梯度上升演算法進行求解。
2 梯度上升演算法
說了半天，梯度上升演算法又是啥？J(θ)太複雜，我們先看個簡單的求極大值的例子。一個看了就會想到高中生活的函式：
來吧，做高中題。這個函式的極值怎麼求？顯然這個函式開口向下，存在極大值，它的函式影象為：
求極值，先求函式的導數：
令導數為0，可求出x=2即取得函式f(x)的極大值。極大值等於f(2)=4
但是真實環境中的函式不會像上面這麼簡單，就算求出了函式的導數，也很難精確計算出函式的極值。此時我們就可以用迭代的方法來做。就像爬坡一樣，一點一點逼近極值。這種尋找最佳擬合引數的方法，就是最優化演算法。爬坡這個動作用數學公式表達即為：
其中，α為步長，也就是學習速率，控制更新的幅度。效果如下圖所示：
比如從(0,0)開始，迭代路徑就是1->2->3->4->…->n，直到求出的x為函式極大值的近似值，停止迭代。我們可以編寫Python3程式碼，來實現這一過程：
# -*- coding:UTF-8 -*-
"""
函式說明:梯度上升演算法測試函式
求函式f(x) = -x^2 + 4x的極大值
Parameters:
無
Returns:
無
Author:
Jack Cui
Blog:
http://blog.csdn.net/c406495762
Zhihu:
https://www.zhihu.com/people/Jack--Cui/
Modify:
2017-08-28
"""
def Gradient_Ascent_test():
def f_prime(x_old): #f(x)的導數
return -2 * x_old + 4
x_old = -1 #初始值，給一個小於x_new的值
x_new = 0 #梯度上升演算法初始值，即從(0,0)開始
alpha = 0.01 #步長，也就是學習速率，控制更新的幅度
presision = 0.00000001 #精度，也就是更新閾值
while abs(x_new - x_old) > presision:
x_old = x_new
x_new = x_old + alpha * f_prime(x_old) #上面提到的公式
print(x_new) #列印最終求解的極值近似值
if __name__ == '__main__':
Gradient_Ascent_test()
程式碼執行結果如下：
結果很顯然，已經非常接近我們的真實極值2了。這一過程，就是梯度上升演算法。那麼同理，J(θ)這個函式的極值，也可以這麼求解。公式可以這麼寫：
由上小節可知J(θ)為：
sigmoid函式為：
那麼，現在我只要求出J(θ)的偏導，就可以利用梯度上升演算法，求解J(θ)的極大值了。
那麼現在開始求解J(θ)對θ的偏導，求解如下(數學推導)：
其中：
再由：
可得：
接下來，就剩下第三部分：
綜上所述：
因此，梯度上升迭代公式為：
知道了，梯度上升迭代公式，我們就可以自己編寫程式碼，計算最佳擬合引數了。
三 Python3實戰
1 資料準備
這就是一個簡單的資料集，沒什麼實際意義。讓我們先從這個簡單的資料集開始學習。先看下資料集有哪些資料：
這個資料有兩維特徵，因此可以將資料在一個二維平面上展示出來。我們可以將第一列資料(X1)看作x軸上的值，第二列資料(X2)看作y軸上的值。而最後一列資料即為分類標籤。根據標籤的不同，對這些點進行分類。
那麼，先讓我們編寫程式碼，看下資料集的分佈情況：
# -*- coding:UTF-8 -*-
import matplotlib.pyplot as plt
import numpy as np
"""
函式說明:載入資料
Parameters:
無
Returns:
dataMat - 資料列表
labelMat - 標籤列表
Author:
Jack Cui
Blog:
http://blog.csdn.net/c406495762
Zhihu:
https://www.zhihu.com/people/Jack--Cui/
Modify:
2017-08-28
"""
def loadDataSet():
dataMat = [] #建立資料列表
labelMat = [] #建立標籤列表
fr = open('testSet.txt') #開啟檔案
for line in fr.readlines(): #逐行讀取
lineArr = line.strip().split() #去回車，放入列表
dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])]) #新增資料
labelMat.append(int(lineArr[2])) #新增標籤
fr.close() #關閉檔案
return dataMat, labelMat #返回
"""
函式說明:繪製資料集
Parameters:
無
Returns:
無
Author:
Jack Cui
Blog:
http://blog.csdn.net/c406495762
Zhihu:
https://www.zhihu.com/people/Jack--Cui/
Modify:
2017-08-28
"""
def plotDataSet():
dataMat, labelMat = loadDataSet() #載入資料集
dataArr = np.array(dataMat) #轉換成numpy的array陣列
n = np.shape(dataMat)[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]) #1為正樣本
else:
xcord2.append(dataArr[i,1]); ycord2.append(dataArr[i,2]) #0為負樣本
fig = plt.figure()
ax = fig.add_subplot(111) #新增subplot
ax.scatter(xcord1, ycord1, s = 20, c = 'red', marker = 's',alpha=.5)#繪製正樣本
ax.scatter(xcord2, ycord2, s = 20, c = 'green',alpha=.5) #繪製負樣本
plt.title('DataSet') #繪製title
plt.xlabel('x'); plt.ylabel('y') #繪製label
plt.show() #顯示
if __name__ == '__main__':
plotDataSet()
執行結果如下：
從上圖可以看出資料的分佈情況。假設Sigmoid函式的輸入記為z，那麼z=w0x0 + w1x1 + w2x2，即可將資料分割開。其中，x0為全是1的向量，x1為資料集的第一列資料，x2為資料集的第二列資料。另z=0，則0=w0 + w1x1 + w2x2。橫座標為x1，縱座標為x2。這個方程未知的引數為w0，w1，w2，也就是我們需要求的迴歸係數(最優引數)。
2 訓練演算法
在編寫程式碼之前，讓我們回顧下梯度上升迭代公式：
將上述公式向量化：
根據向量化的公式，編寫程式碼如下：
# -*- coding:UTF-8 -*-
import numpy as np
"""
函式說明:載入資料
Parameters:
無
Returns:
dataMat - 資料列表
labelMat - 標籤列表
Author:
Jack Cui
Blog:
http://blog.csdn.net/c406495762
Zhihu:
https://www.zhihu.com/people/Jack--Cui/
Modify:
2017-08-28
"""
def loadDataSet():
dataMat = [] #建立資料列表
labelMat = [] #建立標籤列表
fr = open('testSet.txt') #開啟檔案
for line in fr.readlines(): #逐行讀取
lineArr = line.strip().split() #去回車，放入列表
dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])]) #新增資料
labelMat.append(int(lineArr[2])) #新增標籤
fr.close() #關閉檔案
return dataMat, labelMat #返回
"""
函式說明:sigmoid函式
Parameters:
inX - 資料
Returns:
sigmoid函式
Author:
Jack Cui
Blog:
http://blog.csdn.net/c406495762
Zhihu:
https://www.zhihu.com/people/Jack--Cui/
Modify:
2017-08-28
"""
def sigmoid(inX):
return 1.0 / (1 + np.exp(-inX))
"""
函式說明:梯度上升演算法
Parameters:
dataMatIn - 資料集
classLabels - 資料標籤
Returns:
weights.getA() - 求得的權重陣列(最優引數)
Author:
Jack Cui
Blog:
http://blog.csdn.net/c406495762
Zhihu:
https://www.zhihu.com/people/Jack--Cui/
Modify:
2017-08-28
"""
def gradAscent(dataMatIn, classLabels):
dataMatrix = np.mat(dataMatIn) #轉換成numpy的mat
labelMat = np.mat(classLabels).transpose() #轉換成numpy的mat,並進行轉置
m, n = np.shape(dataMatrix) #返回dataMatrix的大小。m為行數,n為列數。
alpha = 0.001 #移動步長,也就是學習速率,控制更新的幅度。
maxCycles = 500 #最大迭代次數
weights = np.ones((n,1))
for k in range(maxCycles):
h = sigmoid(dataMatrix * weights) #梯度上升向量化公式
error = labelMat - h
weights = weights + alpha * dataMatrix.transpose() * error
return weights.getA() #將矩陣轉換為陣列，返回權重陣列
if __name__ == '__main__':
dataMat, labelMat = loadDataSet()
print(gradAscent(dataMat, labelMat))
執行結果如圖所示：
可以看出，我們已經求解出迴歸係數[w0,w1,w2]。
通過求解出的引數，我們就可以確定不同類別資料之間的分隔線，畫出決策邊界。
3 繪製決策邊界
我們已經解出了一組迴歸係數，它確定了不同類別資料之間的分隔線。現在開始繪製這個分隔線，編寫程式碼如下：
# -*- coding:UTF-8 -*-
import matplotlib.pyplot as plt
import numpy as np
"""
函式說明:載入資料
Parameters:
無
Returns:
dataMat - 資料列表
labelMat - 標籤列表
Author:
Jack Cui
Blog:
http://blog.csdn.net/c406495762
Zhihu:
https://www.zhihu.com/people/Jack--Cui/
Modify:
2017-08-28
"""
def loadDataSet():
dataMat = [] #建立資料列表
labelMat = [] #建立標籤列表
fr = open('testSet.txt') #開啟檔案
for line in fr.readlines(): #逐行讀取
lineArr = line.strip().split() #去回車，放入列表
dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])]) #新增資料
labelMat.append(int(lineArr[2])) #新增標籤
fr.close() #關閉檔案
return dataMat, labelMat #返回
"""
函式說明:sigmoid函式
Parameters:
inX - 資料
Returns:
sigmoid函式
Author:
Jack Cui
Blog:
http://blog.csdn.net/c406495762
Zhihu:
https://www.zhihu.com/people/Jack--Cui/
Modify:
2017-08-28
"""
def sigmoid(inX):
return 1.0 / (1 + np.exp(-inX))
"""
函式說明:梯度上升演算法
Parameters:
dataMatIn - 資料集
classLabels - 資料標籤
Returns:
weights.getA() - 求得的權重陣列(最優引數)
Author:
Jack Cui
Blog:
http://blog.csdn.net/c406495762
Zhihu:
https://www.zhihu.com/people/Jack--Cui/
Modify:
2017-08-28
"""
def gradAscent(dataMatIn, classLabels):
dataMatrix = np.mat(dataMatIn) #轉換成numpy的mat
labelMat = np.mat(classLabels).transpose() #轉換成numpy的mat,並進行轉置
m, n = np.shape(dataMatrix) #返回dataMatrix的大小。m為行數,n為列數。
alpha = 0.001 #移動步長,也就是學習速率,控制更新的幅度。
maxCycles = 500 #最大迭代次數
weights = np.ones((n,1))
for k in range(maxCycles):
h = sigmoid(dataMatrix * weights) #梯度上升向量化公式
error = labelMat - h
weights = weights + alpha * dataMatrix.transpose() * error
return weights.getA() #將矩陣轉換為陣列，返回權重陣列
"""
函式說明:繪製資料集
Parameters:
weights - 權重引數陣列
Returns:
無
Author:
Jack Cui
Blog:
http://blog.csdn.net/c406495762
Zhihu:
https://www.zhihu.com/people/Jack--Cui/
Modify:
2017-08-30
"""
def plotBestFit(weights):
dataMat, labelMat = loadDataSet() #載入資料集
dataArr = np.array(dataMat) #轉換成numpy的array陣列
n = np.shape(dataMat)[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]) #1為正樣本
else:
xcord2.append(dataArr[i,1]); ycord2.append(dataArr[i,2]) #0為負樣本
fig = plt.figure()
ax = fig.add_subplot(111) #新增subplot
ax.scatter(xcord1, ycord1, s = 20, c = 'red', marker = 's',alpha=.5)#繪製正樣本
ax.scatter(xcord2, ycord2, s = 20, c = 'green',alpha=.5) #繪製負樣本
x = np.arange(-3.0, 3.0, 0.1)
y = (-weights[0] - weights[1] * x) / weights[2]
ax.plot(x, y)
plt.title('BestFit') #繪製title
plt.xlabel('X1'); plt.ylabel('X2') #繪製label
plt.show()
if __name__ == '__main__':
dataMat, labelMat = loadDataSet()
weights = gradAscent(dataMat, labelMat)
plotBestFit(weights)
執行結果如下：
這個分類結果相當不錯，從上圖可以看出，只分錯了幾個點而已。但是，儘管例子簡單切資料集很小，但是這個方法卻需要大量的計算(300次乘法)。因此下篇文章將對改演算法稍作改進，從而減少計算量，使其可以應用於大資料集上。
四 總結
Logistic迴歸的一般過程：
- 收集資料：採用任意方法收集資料。
- 準備資料：由於需要進行距離計算，因此要求資料型別為數值型。另外，結構化資料格式則最佳。
- 分析資料：採用任意方法對資料進行分析。
- 訓練演算法：大部分時間將用於訓練，訓練的目的是為了找到最佳的分類迴歸係數。
- 測試演算法：一旦訓練步驟完成，分類將會很快。
- 使用演算法：首先，我們需要輸入一些資料，並將其轉換成對應的結構化數值；接著，基於訓練好的迴歸係數，就可以對這些數值進行簡單的迴歸計算，判定它們屬於哪個類別；在這之後，我們就可以在輸出的類別上做一些其他分析工作。
其他：
- Logistic迴歸的目的是尋找一個非線性函式Sigmoid的最佳擬合引數，求解過程可以由最優化演算法完成。
- 本文講述了Logistic迴歸原理以及數學推導過程。
- 下篇文章將講解Logistic迴歸的改進以及Sklearn實戰內容。
- 如有問題，請留言。如有錯誤，還望指正，謝謝！
PS： 如果覺得本篇本章對您有所幫助，歡迎關注、評論、頂！
參考文獻：