Logistic回归实战篇之预测病马死亡率(二)

作 者：崔家华
编 辑：李文臣

三、从疝气病症状预测病马的死亡率

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 np

import random

"""函数说明:sigmoid函数

Parameters:

   inX - 数据

Returns:

   sigmoid函数

Author:    Jack Cui

Blog:

   http://blog.csdn.net/c406495762

Zhihu:

   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 Cui

Blog:    http://blog.csdn.net/c406495762

Zhihu:    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 Cui

Blog:

   http://blog.csdn.net/c406495762

Zhihu:

   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 Cui

Blog:

   http://blog.csdn.net/c406495762

Zhihu:

   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 np

import random

"""函数说明:sigmoid函数

Parameters:

   inX - 数据

Returns:

   sigmoid函数

Author:

   Jack Cui

Blog:

   http://blog.csdn.net/c406495762

Zhihu:

   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 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.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 Cui

Blog:

   http://blog.csdn.net/c406495762

Zhihu:

   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 Cui

Blog:

   http://blog.csdn.net/c406495762

Zhihu:

   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回归实战篇之预测病马死亡率(三)

往期精彩文章推荐

支持向量机原理篇之手撕线性SVM 

章神的私房菜之数据预处理 

深入浅出--基于密度的聚类方法 

贝叶斯系列-贝叶斯与其它统计流派的区别和联系

扫描燕哥微信号，

拉你进机器学习大牛群。

福利满满，名额已不多…

群里目前包括：

清华张长水教授，

清华顾险峰教授，

北大黄铁军教授，

西安电子科技大学焦李成教授，

新加坡南洋理工大学黄广斌教授，

北交李清勇教授

等等……