機器學習python實戰----羅吉斯迴歸

標籤：電腦   就是   提取   pytho   ada   梯度下降法   選擇   ipy   數字   

　　當看到這部分內容的時候我是激動的，因為它終於能跟我之前學習的理論內容聯絡起來了，這部分內容就是對之前羅吉斯迴歸理論部分的代碼實現，所以如果有不甚理解的內容可以返回對照著理論部分來理解，下面我們進入主題----logistic regression

一、sigmoid函數

　　在之前的理論部分我們知道，如果我們需要對某物進行二分類，那麼我們希望輸出函數的值在區間[0,1]，於是我們引入了sigmoid函數。函數的形式為。

曲線圖

　　

根據函數運算式，我們可以用代碼來表示

def sigmoid(Inx):    return 1.0/(1+exp(-Inx))

二、梯度上升法

　　這裡採用的方法為梯度上升法。其實梯度上升法跟梯度下降法的原理是一樣的，不過是換種形式表達。梯度上升法的公式表示為，梯度下降法的公式表示為。這裡雖然看上去只有一個運算子不同，事實上這裡的f(x,θ)也是有些小區別，在後面的代碼中我會解釋一下。

　　這裡首先還是要從文字檔中擷取特徵資料集和標籤，方法在之前已經闡述過幾次了，這部分就直接給出代碼。得到了特徵資料集和標籤，我們就可以用梯度上升法和梯度下降法來尋找最佳參數

def loaddataSet(filename):    dataMat = []    labelsVec = []    file = open(r‘logRegres\testSet.txt‘)    for line in file.readlines():        lineArr = line.strip().split()        dataMat.append([1.0,float(lineArr[0]),float(lineArr[1])])        labelsVec.append(int(lineArr[-1]))    return dataMat,labelsVecdef GradAscent(dataSetIn,classLabels):  #梯度上升法    dataMat = mat(dataSetIn)    labelsVec = mat(classLabels).transpose()    m,n = shape(dataMat)    sigma = ones((n,1))    alpha = 0.001    maxCycles = 500    for i in range(maxCycles):        h = sigmoid(dataMat*sigma)        error = (labelsVec-h)        sigma = sigma + alpha*dataMat.transpose()*error    return sigmadef GradAscent1(dataSetIn,classLabels):   #梯度下降法    dataMat = mat(dataSetIn)    labelsVec = mat(classLabels).transpose()    m,n = shape(dataMat)    sigma = zeros((n,1))    alpha = 0.001    maxCycles = 500    for i in range(maxCycles):        h = sigmoid(dataMat*sigma)        error = (h-labelsVec)        sigma = sigma - alpha*dataMat.transpose()*error    return sigma

這裡的.transpose()方法是對矩陣進行轉置，maxCycles為迭代次數。上面說的梯度下降法與梯度上升法的f(x,θ)的區別在誤差error，兩種方法的error互為相反數。其實只要把梯度下降法的error代入sigma = sigma - alpha*dataMat.transpose()*error就可以發現，將error的負號提出來得到的就跟梯度上升法一樣的式子了。這裡的參數sigma的更新我們省略了推導過程，不熟悉的可以去之前理論部分看看。得到了最佳參數之後，我們就可以畫出決策邊界來看看是否能完美的契合資料集。

def plotBestFit(sigma): #sigma為numpy數組的形式傳入    dataMat,labelsVec = loaddataSet(r‘logRegres\testSet.txt‘)    ArrdataMat = array(dataMat)    n = shape(ArrdataMat)[0]    #sigmaT = GradAscent(dataMat,labelsVec)    #sigma = sigmaT.getA()    xcord1=[];ycord1=[]    xcord2=[];ycord2=[]    for i in range(n):        if labelsVec[i]==1:            xcord1.append(ArrdataMat[i,1])            ycord1.append(ArrdataMat[i,2])        else:            xcord2.append(ArrdataMat[i,1])            ycord2.append(ArrdataMat[i,2])    fig = plt.figure()    ax = plt.subplot(111)    ax.scatter(xcord1,ycord1,s=30,c = ‘red‘,marker = ‘s‘)    ax.scatter(xcord2,ycord2,s=30,c = ‘blue‘)    x = arange(-3.0,3.0,0.1)    y=(-sigma[0]-sigma[1]*x)/sigma[2]    ax.plot(x,y)    plt.show()

這個函數中調用了之前的loaddataSet()函數，這裡需要修改為自己電腦上文字檔所在路徑。然後我們把資料集中的兩個類別0、1分開顯示以便看得更清晰，之後就是繪製決策邊界了。y=(-sigma[0]-sigma[1]*x)/sigma[2]這個運算式是怎麼來的呢？我們的特徵向量X = [x0,x1,x2]，x0=1，運算式中的x相當於x1，y相當於x2，因此我們就用這個式子來表示決策邊界。結果顯示如下

這裡只給出了梯度上升法的結果圖，梯度下降法得到的結果跟這個一模一樣，有興趣的可以自己繪製出來試試。圖中可以看出分類的效果還是蠻不錯的。但是，這種方法存在一個問題，就是計算複雜度較高，對於數百個樣本可以，但是如果有數百萬的樣本呢？下面我們會對上面的演算法進行最佳化。

三、隨機梯度上升法

　　隨機梯度上升法的原理是：一次只用一個樣本點來更新迴歸係數。代碼如下

def stocGradAscent0(dataMat,labelsVec):    ArrdataMat = array(dataMat)    m,n = shape(ArrdataMat)    sigma = ones(n)    alpha = 0.01    for i in range(m):        h = sigmoid(sum(ArrdataMat[i]*sigma))        error = float(labelsVec[i] - h)        sigma = sigma+alpha*error*ArrdataMat[i]    return sigma    

我們可以看出，代入sigmoid函數的不再是整個特徵資料集，而只是一個樣本點，這樣計算量就比之前的方法小了很多了。我們來看下這樣做的效果

紅色部分有將近一半的資料劃分產生的錯誤。那這種方法是不是不好呢？當然不是，這裡我們只運算了一次，還沒有開始迭代，而最佳化之前的方法我們迭代了500次，所以這並不能看出來最佳化後的方法的效果。下面我們會修改代碼，進行150次迭代再來看看效果

def stocGradAscent1(dataMat,labelsVec,numIter = 150):    ArrdataMat = array(dataMat)    m,n = shape(ArrdataMat)    sigma = ones(n)    for i in range(numIter):        for j in range(m):            alpha = 1.0/(1+i+j)+0.01            randIndex = int(random.uniform(0,m))            h = sigmoid(sum(ArrdataMat[randIndex]*sigma))            error = labelsVec[randIndex]-h            sigma = sigma+alpha*error*ArrdataMat[randIndex]    return sigma

這裡的第三個預設參數numIter是迭代次數，預設是150，我們可以自己指定。代碼中的alpha是變化的，會隨著迭代次數不斷減少，但不會減到0，因為有一個常數項，這樣做是為了保證在多次迭代後對新資料仍有一定的影響。這裡的樣本點也需要電腦隨機播放，以減少周期性的波動。我們來看下結果吧

跟最佳化前的結果進行比較，效果還是不錯的，但是這裡的計算複雜度就小了很多很多了。

四、實戰演練

　　這裡給出一些文字檔資料集，有30%的資料缺失。經過資料的預先處理，特徵的缺失我們用實數0來代替，而標籤的缺失我們選擇丟棄該條資料。我們需要從訓練集資料中提取特徵資料集和標籤對模型進行訓練，然後用測試資料集測試訓練效果。當然給出的這裡的測試我們用測試資料特徵集與訓練模型找到的最佳迴歸係數相乘並代入sigmoid函數，我們之前的理論中說過，只要sigmoid函數值大於0.5，我們可以認為其分類是1，否則認為分類為0.

 

def classifyVector(Inx,sigma):    prob = sigmoid(sum(Inx*sigma))    if prob>0.5:        return 1.0    else:        return 0.0    def colicTest():    Trainfile = open(r‘D:\ipython\logRegres\horseColicTraining.txt‘)    Testfile = open(r‘D:\ipython\logRegres\horseColicTest.txt‘)    TrainSet = [];TrainLabels = []    for line in Trainfile.readlines():        line_s = line.strip().split(‘\t‘)        lineArr = []        for i in range(21):            lineArr.append(float(line_s[i]))        TrainSet.append(lineArr)        TrainLabels.append(float(line_s[-1]))    sigma = stocGradAscent1(TrainSet,TrainLabels,500)    error_cnt=0.0;numTestVec = 0    for line1 in Testfile.readlines():        numTestVec+=1        line_s1 = line1.strip().split(‘\t‘)        lineArr1 = []        for j in range(21):           lineArr1.append(float(line_s1[j]))        if int(classifyVector(array(lineArr1),sigma))!=int(line_s1[-1]):            error_cnt+=1    error_rate = float(error_cnt)/numTestVec    print(‘the num of error is %d,the error rate is %f\n ‘%(error_cnt,error_rate))    return error_rate    def multiTest():    numTests = 10    error_sum = 0.0    for i in range(numTests):        error_sum+=colicTest()    print(‘%d iterations the average error rate is %f\n‘%(numTests,float(error_sum)/numTests))

最後一個函數是進行10次測試並取平均值。這裡一定要注意的一個問題是：利用.split()方法分隔得到的是str，也就是說是‘2.0‘,‘1.0‘的形式，這並不是數字，當然在測試時預測的分類(數字)跟標籤中(str)也就是不相等的。所以這裡一定要注意，如果你也產生的error_rate=1的結果，可以查看是不是犯了跟我一樣的錯誤！

這裡的預測結果中有34.3%的錯誤率，這個結果其實並不差，因為我們使用的資料集中有30%的資料缺失。

 

資料集及代碼下載 http://pan.baidu.com/s/1qYHT3i8

 

機器學習python實戰----羅吉斯迴歸