【機器學習】logistic迴歸原理分析及python實現

【機器學習】logistic迴歸原理分析及python實現

【機器學習】logistic迴歸原理分析及python實現

1.sigmoid函式和logistic迴歸分類器

2.梯度上升最優化演算法

3.資料中的缺失項處理

4.logistic實現馬疝氣病預測

        首先闡述logistic迴歸的定義,然後介紹一些最優化演算法,其中包括基本的梯度上升演算法和改進的隨機梯度上升演算法,這些最優化演算法用於分類器的訓練,最後給出logistic迴歸例項,預測一匹有疝氣病的馬是否被治癒(二分類)。

一.sigmoid函式和logistic迴歸分類器

1.什麼是迴歸?

    假設現在有一些資料點,我們用一條直線來對這些點進行擬合(該線稱為最佳擬合直線),這個擬合過程就稱做迴歸。

2.sigmoid函式

    我們想要的函式應該是,接受所有的輸入,然後預測出類別。例如,在二類情況下,輸出0和1。像單位階躍函式就可以實現,但該函式在跳躍點上從0->1是瞬間跳躍,這個瞬間跳躍過程有時很難處理。於是另一個有類似性質且數學上更易處理的函式出現了—–sigmoid函式。

         sigmoid函式表示式:

                            

        實現:    

"""
函式:sigmoid函式
"""
def sigmoid(z):
return 1.0/(1 np.exp(-z) )

3.sigmoid函式如何用於二分類?

        為了實現logistic迴歸分類器,可以在每個特徵上都乘以一個迴歸係數w,然後把所有的結果值相加得到z值,將這個z值帶入sigmoid函式中,會輸出一個在【0,1】內的數值。分類:z>0.5,輸出1;z<0.5,輸出0。

    1)輸入樣本:  X=(x0,x1……xn)  

    2)如何將樣本值轉化為sigmoid的輸入?x-> z

       相應的迴歸係數W=(w0,w1……wn),樣本特徵值與相應係數相乘求和:

    3)帶入sigmoid函式: 


4.logistic 迴歸的優點與缺點?

        優點:計算代價不高,易於理解與實現(簡單)。 
        缺點:容易欠擬合,分類精度可能不高。 

        適用資料型別:數值型和標稱型資料。

5.logistic迴歸的一般過程? 

    1、收集資料:任何方式 
    2、準備資料:由於要計算距離,因此要求資料都是數值型的,另外結構化資料格式最佳。 
    3、分析資料:採用任一方是對資料進行分析 
    4、訓練演算法:大部分時間將用於訓練,訓練的目的為了找到最佳的分類迴歸係數 
    5、測試演算法:一旦訓練步驟完成,分類將會很快 
    6、使用演算法:首先,我們需要輸入一些資料,並將其轉化成對應的結構化數值;接著基於訓練好的迴歸係數就可以對這些數值進行簡單的迴歸計算,判定它們屬於哪一類別;在這之後,我們就可以在輸出的類別上做一些其他的分析工作。

二.梯度上升最優化演算法

1.如何確定最佳迴歸係數?

        確定了分類器的函式形式之後,現在的問題變成了:最佳迴歸係數W=(w0,w1……wn)是多少?如何確定它們的大小?

為了尋找最佳引數,需要用到最優化理論的知識。

1.1梯度上升法

        梯度上升法的基本思想是:要找到某函式的最大值,最好的方法是沿著該函式的梯度方向去探尋。

        函式f(x,y)的梯度記為:

                

        這個梯度意味著f沿著x方向移動,沿著y方向移動。其中,函式f(x,y)必須在待計算的點上有定義並且可微。

梯度上升演算法公式  其中,α為學習步長。

        梯度上升演算法用來求函式的最大值,梯度下降演算法(將上式的 改為 – 號)用來求函式的最小值。

何時停止迭代?

        該公式一直被迭代,直到達到某個停止條件為止,比如迭代次數達到某指定值或者誤差達到允許的範圍內。

    

    梯度上升演算法的python實現:

"""
函式:梯度上升演算法
"""
def gradAscent(dataMat,labelMat):
dataSet=np.mat(dataMat)                          # m*n
labelSet=np.mat(labelMat).transpose()            # 1*m->m*1
m,n=np.shape(dataSet)                            # m*n: m個樣本,n個特徵
alpha=0.001                                      # 學習步長
maxCycles=500                                    # 最大迭代次數
weights=np.ones( (n,1) )
for i in range(maxCycles):
y=sigmoid(dataSet * weights)                 # 預測值
error=labelSet -y
weights=weights  alpha *dataSet.transpose()*error
#print(type(weights))
return weights.getA(),weights  ##getA():將Mat轉化為ndarray,因為mat不能用index

    畫出決策邊界:

"""
函式:畫出決策邊界
"""
def plotBestFit(weight):
dataMat, labelMat = loadDataSet()
dataArr=np.array(dataMat)
m,n=np.shape(dataArr)
x1=[]           #x1,y1:類別為1的特徵
x2=[]           #x2,y2:類別為2的特徵
y1=[]
y2=[]
for i in range(m):
if (labelMat[i])==1:
x1.append(dataArr[i,1])
y1.append(dataArr[i,2])
else:
x2.append(dataArr[i,1])
y2.append(dataArr[i,2])
fig=plt.figure()
ax=fig.add_subplot(111)
ax.scatter(x1,y1,s=30,c='red',marker='s')
ax.scatter(x2,y2,s=30,c='green')
#畫出擬合直線
x=np.arange(-3.0, 3.0, 0.1)
y=(-weights[0]-weights[1]*x)/weights[2]    #直線滿足關係:0=w0*1.0 w1*x1 w2*x2
ax.plot(x,y)
plt.xlabel('X1')
plt.ylabel('X2')
plt.show()

    梯度上升演算法的實際效果:

                        

    分析:

            梯度上升演算法的分類效果挺好,但這個方法需要大量的計算。因為每次更新迴歸係數時都需要遍歷整個資料集,該方法在處理小樣本、小資料集時尚可,但如果有數十億樣本和成千上萬的特徵,那麼該方法的計算複雜度就太高了。下面將介紹一種改進方法—-隨機梯度上升。

1.2一種改進方法—-隨機梯度上升

    隨機梯度上升,一次僅用一個樣本點來更新迴歸係數。

    

實現:

"""
函式:隨機梯度上升演算法0.0
改進:每次用一個樣本來更新迴歸係數
"""
def stocGradAscent0(dataMat,labelMat):
m, n = np.shape(dataMat)  # m*n: m個樣本,n個特徵
alpha = 0.001  # 學習步長
maxCycles=500
weights = np.ones(n)
for cycle in range(maxCycles):
for i in range(m):
y = sigmoid(sum(dataMat[i] * weights) )  # 預測值
error = labelMat[i] - y
weights = weights   alpha  * error* dataMat[i]
# print(type(weights))
return weights

    梯度上升演算法的實際效果:

                    

    分析:

    一個判斷優化演算法優劣的可靠方法是看它是否收斂,即引數是否達到了穩定值。在隨機上升演算法中,在500次迭代後差不多達到穩定,即收斂。

1.3改進的隨機梯度上升法

def stocGradAscent1(dataMat,labelMat):
m, n = np.shape(dataMat)  # m*n: m個樣本,n個特徵
maxCycles = 150
weights = np.ones(n)
for cycle in range(maxCycles):
dataIndex=list( range(m))
for i in range(m):
alpha = 4 / (1.0   cycle   i)   0.01                 # 學習步長
randIndex=int(np.random.uniform(0,len(dataIndex) ))  #隨機選取樣本
y = sigmoid(sum(dataMat[randIndex] * weights ))          # 預測值
error = labelMat[randIndex] - y
weights = weights   alpha  * error * dataMat[i]
del(dataIndex[randIndex])
# print(type(weights))
return weights

預設150次迭代,結果如下:

                        

2.視覺化:迭代次數與迴歸係數的關係

    此部分參考:https://blog.csdn.net/c406495762/article/details/77851973(感謝好帖!)

    分別展示梯度上升法、隨機梯度上升法、改進的隨機梯度上升法的迭代次數與迴歸次數的關係如下:

    疑問:實現的迭代過程中發現前兩種方法未實現收斂,但從分類結果來看,是挺好的,這是為什麼呢?若有人復現出與我不同的結果,請留言或私信告知,感謝!

    程式碼如下

# -*- coding:utf-8 -*-
"""
@author:Lisa
@file:visual_iteration_process.py
@note:視覺化迭代次數與迴歸係數的關係
@time:2018/7/12 0012下午 8:08
"""
import numpy as np
import matplotlib.pyplot as plt
"""
函式:載入資料集
"""
def loadDataSet():
dataMat=[]  #列表list
labelMat=[]
txt=open('testSet.txt')
for line in txt.readlines():
lineArr=line.strip().split()        #strip():返回一個帶前導和尾隨空格的字串的副本
#split():預設以空格為分隔符,空字串從結果中刪除
dataMat.append( [1.0, float(lineArr[0]), float(lineArr[1]) ])  #將二維特徵擴充套件到三維,第一維都設定為1.0
labelMat.append(int(lineArr[2]) )
return dataMat,labelMat
"""
函式:sigmoid函式
"""
def sigmoid(z):
return 1.0/(1 np.exp(-z) )
"""
函式:梯度上升演算法
"""
def gradAscent(dataMat,labelMat):
dataSet=np.mat(dataMat)                          # m*n
labelSet=np.mat(labelMat).transpose()            # 1*m->m*1
m,n=np.shape(dataSet)                            # m*n: m個樣本,n個特徵
alpha=0.001                                      # 學習步長
maxCycles=500                                    # 最大迭代次數
weights=np.ones((n,1))
weights_array=np.array([])
for i in range(maxCycles):
y=sigmoid(dataSet * weights)                 # 預測值
error=labelSet -y
weights=weights  alpha *dataSet.transpose()*error
weights_array=np.append(weights_array,weights)
weights_array = weights_array.reshape(maxCycles, n)
#print(weights_array)
return weights.getA(),weights_array  ##getA():將Mat轉化為ndarray,因為mat不能用index
"""
函式:隨機梯度上升演算法0.0
改進:每次用一個樣本來更新迴歸係數
"""
def stocGradAscent0(dataMat,labelMat):
m, n = np.shape(dataMat)  # m*n: m個樣本,n個特徵
alpha = 0.001  # 學習步長
maxCycles=500
weights = np.ones(n)
weights_array = np.array([])
for cycle in range(maxCycles):
for i in range(m):
y = sigmoid(sum(dataMat[i] * weights) )  # 預測值
error = labelMat[i] - y
weights = weights   alpha  * error* dataMat[i]
weights_array = np.append(weights_array, weights)
# print(type(weights))
weights_array = weights_array.reshape(maxCycles*m, n)
return weights,weights_array
"""
函式:改進的隨機梯度上升法1.0
改進:1.alpha隨著迭代次數不斷減小,但永遠不會減小到0
2.通過隨機選取樣本來更新迴歸係數
"""
def stocGradAscent1(dataMat,labelMat):
m, n = np.shape(dataMat)  # m*n: m個樣本,n個特徵
maxCycles = 500
weights = np.ones(n)
weights_array=np.array([])
for cycle in range(maxCycles):
dataIndex=list(range(m))
for i in range(m):
alpha = 4 / (1.0   cycle   i)   0.01                 # 學習步長
randIndex=int(np.random.uniform(0,len(dataIndex) ))  #隨機選取樣本
y = sigmoid(sum(dataMat[randIndex] * weights ))          # 預測值
error = labelMat[randIndex] - y
weights = weights   alpha  * error * dataMat[randIndex]
del(dataIndex[randIndex])
# print(type(weights))
weights_array = np.append(weights_array, weights)
weights_array = weights_array.reshape(maxCycles*m, n)
return weights,weights_array
"""
函式:視覺化迭代過程
"""
def visual(weights_array1,weights_array2,weights_array3):
import matplotlib as mpl
# 指定預設字型
mpl.rcParams['font.sans-serif'] = ['SimHei'] #顯示漢字
mpl.rcParams['axes.unicode_minus'] = False  # 能夠顯示符號(負號)
fig,axs=plt.subplots(3,3,figsize=(15,10))
x1=np.arange(0,len(weights_array1),1)
#繪製w0和迭代次數的關係
axs[0][0].plot(x1,weights_array1[:,0])
axs[0][0].set_title('梯度上升演算法:迴歸係數與迭代次數的關係')
axs[0][0].set_ylabel('w0',)
# 繪製w1和迭代次數的關係
axs[1][0].plot(x1, weights_array1[:, 1])
axs[1][0].set_ylabel('w1')
# 繪製w2和迭代次數的關係
axs[2][0].plot(x1, weights_array1[:, 2])
axs[2][0].set_ylabel('w2')
axs[2][0].set_xlabel('迭代次數')
x2 = np.arange(0, len(weights_array2), 1)
# 繪製w0和迭代次數的關係
axs[0][1].plot(x2, weights_array2[:, 0])
axs[0][1].set_title('隨機梯度上升演算法')
axs[0][1].set_ylabel('w0')
# 繪製w1和迭代次數的關係
axs[1][1].plot(x2, weights_array2[:, 1])
axs[1][1].set_ylabel('w1')
# 繪製w2和迭代次數的關係
axs[2][1].plot(x2, weights_array2[:, 2])
axs[2][1].set_ylabel('w2')
axs[2][1].set_xlabel('迭代次數')
x3 = np.arange(0, len(weights_array3), 1)
# 繪製w0和迭代次數的關係
axs[0][2].plot(x3, weights_array3[:, 0])
axs[0][2].set_title('改進的隨機梯度上升演算法')
axs[0][2].set_ylabel('w0')
# 繪製w1和迭代次數的關係
axs[1][2].plot(x3, weights_array3[:, 1])
axs[1][2].set_ylabel('w1')
# 繪製w2和迭代次數的關係
axs[2][2].plot(x3, weights_array3[:, 2])
axs[2][2].set_ylabel('w2')
axs[2][2].set_xlabel('迭代次數')
plt.show()
if __name__ == '__main__':
dataMat, labelMat = loadDataSet()
weights1, weights_array1 = gradAscent(dataMat, labelMat)
weights2, weights_array2 = stocGradAscent0(np.array(dataMat), labelMat)
weights3, weights_array3 = stocGradAscent1(np.array(dataMat), labelMat)
visual(weights_array1, weights_array2,weights_array3)

三.資料中的缺失項處理

        資料中的缺失值是一個非常棘手的問題,很多文獻都致力於解決這個問題。那麼,資料缺失究竟帶來了什麼問題?假設有100個樣本和20個特徵,這些資料都是機器收集回來的。若機器上的某個感測器損壞導致一個特徵無效時該怎麼辦?它們是否還可用?答案是肯定的。因為有時候資料相當昂貴,扔掉和重新獲取都是不可取的,所以必須採用一些方法來解決這個問題。下面給出了一些可選的做法:

  • 使用可用特徵的均值來填補缺失值;
  • 使用特殊值來填補缺失值,如-1;
  • 忽略有缺失值的樣本;
  • 使用相似樣本的均值添補缺失值;
  • 使用另外的機器學習演算法預測缺失值

四.logistic實現馬疝氣病預測

    1.預處理

    由於該資料集存在30%的缺失,那麼首先必須對資料集進行預處理。

    預處理資料做兩件事:

  • 如果測試集中一條資料的特徵值已經確實,那麼我們選擇實數0來替換所有缺失值,因為本文使用Logistic迴歸。因此這樣做不會影響迴歸係數的值。sigmoid(0)=0.5,即它對結果的預測不具有任何傾向性。
  • 如果測試集中一條資料的類別標籤已經缺失,那麼我們將該類別資料丟棄,因為類別標籤與特徵不同,很難確定採用某個合適的值來替換。

    所有處理後的資料集可見:(感謝Jack-Cui大神處理的資料集)

        https://github.com/LisaPig/ML_LogisticRegression

2.完整程式碼

# -*- coding:utf-8 -*-
"""
@author:Lisa
@file:horse_colic.py
@note:疝氣馬能否痊癒預測
@time:2018/7/12 0012下午 10:11
"""
import numpy as np
import logisticRegression as LR
"""
函式:sigmoid迴歸分類
"""
def classifyVector(dataIn,weights):
h=LR.sigmoid(sum(dataIn*weights))
if h>0.5:
return 1.0
else:
return 0.0
"""
函式:疝氣預測
"""
def colicTest():
trainData=open('data\horseColicTraining.txt')
testData = open('data\horseColicTest.txt')
trainSet=[]
trainLabel=[]
for line in trainData.readlines():
curLine=line.strip().split('\t')
lineArr=[]
for i in range(21):
lineArr.append(float (curLine[i]))
trainSet.append(lineArr)
trainLabel.append(float(curLine[21]))
trainWeights=LR.stocGradAscent1(np.array(trainSet),trainLabel)
errorCount=0
numTestVec=0
testSet = []
testLabel = []
for line in testData.readlines():
numTestVec =1
curLine=line.strip().split('\t')
lineArr=[]
for i in range(21):
lineArr.append(float (curLine[i]))
if int(classifyVector(np.array(lineArr),trainWeights))!=int(curLine[21]):
errorCount =1
errorRate=float(errorCount/numTestVec)
print("the error rate is: %f" % errorRate)
return errorRate
def multiTest():
numTests=10
errorSum=0
for k in range(numTests):
errorSum =colicTest()
print("after %d iterations the average error rate is: %f"
%(numTests,errorSum/float(numTests)))
if __name__ == '__main__':
multiTest()

執行結果:


———————————————————    end    ———————————————————————

參考:

1.《機器學習實戰》

2. https://blog.csdn.net/c406495762/article/details/77851973