作 者：崔家華 編 輯：李文臣

三、從疝氣病症狀預測病馬的死亡率

1、實戰背景

本次實戰內容，將使用Logistic迴歸來預測患疝氣病的馬的存活問題。原始資料集下載地址：http://archive.ics.uci.edu/ml/datasets/Horse+Colic

這裡的資料包含了368個樣本和28個特徵。這種病不一定源自馬的腸胃問題，其他問題也可能引發馬疝病。該資料集中包含了醫院檢測馬疝病的一些指標，有的指標比較主觀，有的指標難以測量，例如馬的疼痛級別。另外需要說明的是，除了部分指標主觀和難以測量外，該資料還存在一個問題，資料集中有30%的值是缺失的。下面將首先介紹如何處理資料集中的資料缺失問題，然後再利用Logistic迴歸和隨機梯度上升演算法來預測病馬的生死。

2、準備資料

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

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

預處理資料做兩件事：

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

原始的資料集經過處理，儲存為兩個檔案：horseColicTest.txt和horseColicTraining.txt。已經處理好的“乾淨”可用的資料集下載地址：

• https://github.com/Jack-Cherish/Machine-Learning/blob/master/Logistic/horseColicTraining.txt
• https://github.com/Jack-Cherish/Machine-Learning/blob/master/Logistic/horseColicTest.txt

有了這些資料集，我們只需要一個Logistic分類器，就可以利用該分類器來預測病馬的生死問題了。

3、使用Python構建Logistic迴歸分類器

在使用Sklearn構建Logistic迴歸分類器之前，我們先用自己寫的改進的隨機梯度上升演算法進行預測，先熱熱身。使用Logistic迴歸方法進行分類並不需要做很多工作，所需做的只是把測試集上每個特徵向量乘以最優化方法得來的迴歸係數，再將乘積結果求和，最後輸入到Sigmoid函式中即可。如果對應的Sigmoid值大於0.5就預測類別標籤為1，否則為0。

# -*- coding:UTF-8 -*-import numpy as npimport random"""函式說明:sigmoid函式Parameters:    inX - 資料Returns:    sigmoid函式Author:    Jack CuiBlog:    http://blog.csdn.net/c406495762Zhihu:    https://www.zhihu.com/people/Jack--Cui/Modify:    2017-09-05"""def sigmoid(inX):    return 1.0 / (1 + np.exp(-inX))"""函式說明:改進的隨機梯度上升演算法Parameters:    dataMatrix - 資料陣列    classLabels - 資料標籤    numIter - 迭代次數Returns:    weights - 求得的迴歸係數陣列(最優引數)Author:    Jack CuiBlog:    http://blog.csdn.net/c406495762Zhihu:    https://www.zhihu.com/people/Jack--Cui/Modify:    2017-09-05"""def stocGradAscent1(dataMatrix, classLabels, numIter=150):
     #返回dataMatrix的大小。m為行數,n為列數。    weights = np.ones(n)   #儲存每次更新的迴歸係數
    for j in range(numIter):                                           
        dataIndex = list(range(m))
        for i in range(m):           
            alpha = 4/(1.0+j+i)+0.01 
    #降低alpha的大小，每次減小1/(j+i)。
            randIndex = int(random.uniform(0,len(dataIndex)))
     #隨機選取樣本
            h = sigmoid(sum(dataMatrix[randIndex]*weights)) 
     #選擇隨機選取的一個樣本，計算h
            error = classLabels[randIndex] - h 
      #計算誤差
            weights = weights + alpha
                  * error * dataMatrix[randIndex] 
                           #更新迴歸係數
            del(dataIndex[randIndex]) 
                        #刪除已經使用的樣本    return weights       #返回"""函式說明:使用Python寫的Logistic分類器做預測Parameters:    無Returns:    無Author:    Jack CuiBlog:    http://blog.csdn.net/c406495762Zhihu:    https://www.zhihu.com/people/Jack--Cui/Modify:    2017-09-05"""def colicTest():    frTrain = open('horseColicTraining.txt')
      #開啟訓練集
    frTest = open('horseColicTest.txt')
      #開啟測試集
    trainingSet = []; trainingLabels = []
    for line in frTrain.readlines():
        currLine = line.strip().split('t')
        lineArr = []
        for i in range(len(currLine)-1):
            lineArr.append(float(currLine[i]))
        trainingSet.append(lineArr)
        trainingLabels.append(float(currLine[-1]))
    trainWeights =
    stocGradAscent1(np.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(len(currLine)-1):
            lineArr.append(float(currLine[i]))
        if int(classifyVector
        (np.array(lineArr), trainWeights))!= int(currLine[-1]):
            errorCount += 1
    errorRate = (float(errorCount)/numTestVec) * 100  
     #錯誤率計算    print("測試集錯誤率為: %.2f%%" % errorRate)
"""函式說明:分類函式Parameters:    inX - 特徵向量    weights - 迴歸係數Returns:    分類結果Author:    Jack CuiBlog:    http://blog.csdn.net/c406495762Zhihu:    https://www.zhihu.com/people/Jack--Cui/Modify:    2017-09-05"""def classifyVector(inX, weights):    prob = sigmoid(sum(inX*weights))
    if prob > 0.5: return 1.0
    else: return 0.0if __name__ == '__main__':
    colicTest()

執行結果如下：

錯誤率還是蠻高的，而且耗時1.9s，並且每次執行的錯誤率也是不同的，錯誤率高的時候可能達到40%多。為啥這樣？首先，因為資料集本身有30%的資料缺失，這個是不能避免的。另一個主要原因是，我們使用的是改進的隨機梯度上升演算法，因為資料集本身就很小，就幾百的資料量。用改進的隨機梯度上升演算法顯然不合適。讓我們再試試梯度上升演算法，看看它的效果如何？

# -*- coding:UTF-8 -*-import numpy as npimport random"""函式說明:sigmoid函式Parameters:    inX - 資料Returns:    sigmoid函式Author:    Jack CuiBlog:    http://blog.csdn.net/c406495762Zhihu:    https://www.zhihu.com/people/Jack--Cui/Modify:    2017-09-05"""def sigmoid(inX):    return 1.0 / (1 + np.exp(-inX))
"""函式說明:梯度上升演算法Parameters:    dataMatIn - 資料集    classLabels - 資料標籤Returns:    weights.getA() - 求得的權重陣列(最優引數)Author:    Jack CuiBlog:    http://blog.csdn.net/c406495762Zhihu:    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.01 
    #移動步長,也就是學習速率,控制更新的幅度。
    maxCycles = 500 
    #最大迭代次數
    weights = np.ones((n,1))
       #梯度上升向量化公式
        error = labelMat - h
        weights = weights + alpha * 
                         dataMatrix.transpose() * error    return weights.getA()  
    #將矩陣轉換為陣列，並返回
"""函式說明:使用Python寫的Logistic分類器做預測Parameters:    無Returns:    無Author:    Jack CuiBlog:    http://blog.csdn.net/c406495762Zhihu:    https://www.zhihu.com/people/Jack--Cui/Modify:    2017-09-05"""def colicTest():    frTrain = open('horseColicTraining.txt')#開啟訓練集
    frTest = open('horseColicTest.txt')#開啟測試集
    trainingSet = []; trainingLabels = []
    for line in frTrain.readlines():
        currLine = line.strip().split('t')
        lineArr = []
        for i in range(len(currLine)-1):
            lineArr.append(float(currLine[i]))
        trainingSet.append(lineArr)
        trainingLabels.append(float(currLine[-1]))
    trainWeights = gradAscent(np.array(trainingSet), trainingLabels)
      #使用改進的隨即上升梯度訓練
    errorCount = 0; numTestVec = 0.0
    for line in frTest.readlines():
        numTestVec += 1.0
        currLine = line.strip().split('t')
        lineArr =[]
        for i in range(len(currLine)-1):
            lineArr.append(float(currLine[i]))
        if int(classifyVector
        (np.array(lineArr), trainWeights[:,0]))!= int(currLine[-1]):
            errorCount += 1
    errorRate = (float(errorCount)/numTestVec) * 100

#錯誤率計算    print("測試集錯誤率為: %.2f%%" % errorRate)
"""函式說明:分類函式Parameters:    inX - 特徵向量    weights - 迴歸係數Returns:     分類結果Author:    Jack CuiBlog:    http://blog.csdn.net/c406495762Zhihu:    https://www.zhihu.com/people/Jack--Cui/Modify:    2017-09-05"""def classifyVector(inX, weights):    prob = sigmoid(sum(inX*weights))
    if prob > 0.5: return 1.0
    else: return 0.0if __name__ == '__main__':
    colicTest()

執行結果如下：

可以看到演算法耗時減少了，錯誤率穩定且較低。很顯然，使用隨機梯度上升演算法，反而得不償失了。所以可以得到如下結論：

• 當資料集較小時，我們使用梯度上升演算法
• 當資料集較大時，我們使用改進的隨機梯度上升演算法

對應的，在Sklearn中，我們就可以根據資料情況選擇優化演算法，比如資料較小的時候，我們使用liblinear，資料較大時，我們使用sag和saga。

本系列篇章：

Logistic迴歸實戰篇之預測病馬死亡率(一)

Logistic迴歸實戰篇之預測病馬死亡率(二)

Logistic迴歸實戰篇之預測病馬死亡率(三)