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

一、sigmoid函數

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

曲線圖

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

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

二、梯度上升法

　　這裏采用的方法為梯度上升法。其實梯度上升法跟梯度下降法的原理是一樣的，不過是換種形式表達。梯度上升法的公式表示為

　　這裏首先還是要從文本文件中獲取特征數據集和標簽，方法在之前已經闡述過幾次了，這部分就直接給出代碼。得到了特征數據集和標簽，我們就可以用梯度上升法和梯度下降法來尋找最佳參數

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,labelsVec

def 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 sigma

def 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實戰----邏輯回歸