Logistic回归应用-预测马的死亡率

Logistic回归应用-预测马的死亡率

本文所有代码均来自《机器学习实战》，数据也是

本例中的数据有以下几个特征：

1. 部分指标比较主观、难以很好的定量测量，例如马的疼痛级别
2. 数据集中有30%的数据是缺失的，例如：

如何处理数据缺失情况

数据是很宝贵的，对于有缺失的数据我们不能直接丢弃，而应该想办法将其修补后再利用起来

常用的方法有：

1. 使用可用样本中该特征的均值来填补缺失样本中缺失特征的值
2. 使用特殊值来填补缺失值，如-1
3. 忽略有缺失值的样本
4. 使用相似样本的均值填补缺失值
5. 使用另外的机器学习算法来预测缺失值

但如果是Label信息丢失的话，一个简单的办法就是将该data丢弃，因为Label和attributes不一样，很难确定采用某个合适的值来替换缺失值。但是这种方法在类似KNN这样的、样本数量和样本完整度对结果有较为明显影响的算法中就不太好了。

---

回到本例，这里是采用了将缺失数据置换为0的操作。因为在Logistic回归的梯度下降法中，如果某个样本的一个attribute为0的话，那么就不使用它进行weights的更新了，这显然是一种很好的结果。另外，由于sigmoid(0)=0.5,既不偏向0也不偏向1，所以这个特征设置为0的话，该特征不会对结果的预测产生偏向性。

这是未处理的数据：

这是处理后的数据：

关于每个attribute代表什么意思，可以在这个网站上查到：http://archive.ics.uci.edu/ml/datasets/Horse+Colic。反正就是关于马的一些特征，由上图我们可以看到，这里我们并没有使用原数据中的全部特征，只选了一部分。

---

将Logistic回归用于分类

之前一篇文章：https://www.cnblogs.com/jiading/p/11660271.html。中，我们只是画出了Logistic回归生成的边界，但是并没有完善分类算法，这里我们完成这个工作,代码如下：

#用于返回分类结果（0或1）
def classifyVector(inX, weights):
    prob = sigmoid(sum(inX*weights))
    if prob > 0.5: return 1.0
    else: return 0.0

def colicTest():
    #训练集和测试集
    frTrain = open('horseColicTraining.txt'); frTest = open('horseColicTest.txt')
    trainingSet = []; trainingLabels = []
    #调用训练集进行训练
    for line in frTrain.readlines():
        currLine = line.strip().split('	')
        lineArr =[]
        #一共有21个特征
        for i in range(21):
            #lineArr以数组形式保存由字符转换为浮点数的一条数据
            lineArr.append(float(currLine[i]))
        trainingSet.append(lineArr)
        trainingLabels.append(float(currLine[21]))
    trainWeights = stocGradAscent1(array(trainingSet), trainingLabels, 1000)
    errorCount = 0; numTestVec = 0.0
    #调用测试集进行测试
    for line in frTest.readlines():
        numTestVec += 1.0
        currLine = line.strip().split('	')
        lineArr =[]
        for i in range(21):
            lineArr.append(float(currLine[i]))
        if int(classifyVector(array(lineArr), trainWeights))!= int(currLine[21]):
            errorCount += 1
    errorRate = (float(errorCount)/numTestVec)#返回误差率
    print("the error rate of this test is: %f" % errorRate)
    return errorRate

完整代码（实现+测试）如下：

'''
Created on Oct 27, 2010
Logistic Regression Working Module
@author: Peter
'''
'''
注意，本例中默认data在testSet.txt中且有两个attributes，
使用该模型的时候要注意根据实际使用的attributes个数修改代码
'''
from numpy import *

def loadDataSet():#加载数据
    dataMat = []; labelMat = []
    fr = open('testSet.txt')
    for line in fr.readlines():
        lineArr = line.strip().split()
        #Python strip() 方法用于移除字符串头尾指定的字符（默认为空格）或字符序列。
        #注意，在这里为data添加了1的一列
        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
        #数据有两个attributes
        labelMat.append(int(lineArr[2]))
    return dataMat,labelMat

def sigmoid(inX):
    return 1.0/(1+exp(-inX))
#梯度下降方法
def gradAscent(dataMatIn, classLabels):
    #转换为numpy矩阵类型，便于优化和并行计算
    dataMatrix = mat(dataMatIn)             #convert to NumPy matrix
    labelMat = mat(classLabels).transpose() #convert to NumPy matrix
    #将Label转置为列向量
    m,n = shape(dataMatrix)#m行n列
    alpha = 0.001#学习速度
    maxCycles = 500#最大学习步数
    weights = ones((n,1))#将来的输出
    for k in range(maxCycles):              #heavy on matrix operations
        h = sigmoid(dataMatrix*weights)     #matrix mult
        error = (labelMat - h)              #vector subtraction
        weights = weights + alpha * dataMatrix.transpose()* error #matrix mult
    return weights

def plotBestFit(weights):
    import matplotlib.pyplot as plt
    dataMat,labelMat=loadDataSet()#加载源数据，用于画图
    dataArr = array(dataMat)#将data向量化
    n = shape(dataArr)[0] #n为data的个数
    xcord1 = []; ycord1 = []
    xcord2 = []; ycord2 = []
    #给data分类
    for i in range(n):
        if int(labelMat[i])== 1:
            xcord1.append(dataArr[i,1]); ycord1.append(dataArr[i,2])
        else:
            xcord2.append(dataArr[i,1]); ycord2.append(dataArr[i,2])
    fig = plt.figure()
    ax = fig.add_subplot(111)
    #用不同的颜色画出来
    ax.scatter(xcord1, ycord1, s=30, c='red', marker='s')
    ax.scatter(xcord2, ycord2, s=30, c='green')
    #arange([start,] stop[, step,], dtype=None)
    #根据start与stop指定的范围以及step设定的步长，生成一个 ndarray。
    '''
    这一段是在干什么呢？其实就是想要画一条直线出来，但是只能画点，由点组成线
    '''
    x = arange(-3.0, 3.0, 0.1)
    y = (-weights[0]-weights[1]*x)/weights[2]
    ax.plot(x, y)
    plt.xlabel('X1'); plt.ylabel('X2');
    plt.show()
#随机梯度下降法，只对整个数据集迭代了一次，显然不够
def stocGradAscent0(dataMatrix, classLabels):
    m,n = shape(dataMatrix)
    alpha = 0.01
    weights = ones(n)   #initialize to all ones
    for i in range(m):
        h = sigmoid(sum(dataMatrix[i]*weights))
        error = classLabels[i] - h
        #这里其实不是真正的随机，样本是按顺序取的
        weights = weights + alpha * error * dataMatrix[i]
    return weights
#改进的随机梯度下降法
def stocGradAscent1(dataMatrix, classLabels, numIter=150):
    m,n = shape(dataMatrix)
    weights = ones(n)   #initialize to all ones
    for j in range(numIter):#对样本迭代150次
        dataIndex = list(range(m))#dataIndex保存了所有在本次迭代中还没有使用过的data
        for i in range(m):
            '''
            改进1：学习速度在每次迭代后都会调整，越往后学习速度越小，这会缓解在最优点附近的波动.
            由于我们在alpha之后加了一个常数，这使得alpha即使再迭代也不会减小到0.
            在降低alpha的函数中，j是迭代次数，i是data下标，这样可以避免参数的严格下降
            '''
            alpha = 4/(1.0+j+i)+0.0001    #apha decreases with iteration, does not
            #改进2：随机选择用来计算的data
            randIndex = int(random.uniform(0,len(dataIndex)))#go to 0 because of the constant
            h = sigmoid(sum(dataMatrix[randIndex]*weights))
            error = classLabels[randIndex] - h
            weights = weights + alpha * error * dataMatrix[randIndex]
            del(dataIndex[randIndex])
    return weights
#用于返回分类结果（0或1）
def classifyVector(inX, weights):
    prob = sigmoid(sum(inX*weights))
    if prob > 0.5: return 1.0
    else: return 0.0

def colicTest():
    #训练集和测试集
    frTrain = open('horseColicTraining.txt'); frTest = open('horseColicTest.txt')
    trainingSet = []; trainingLabels = []
    #调用训练集进行训练
    for line in frTrain.readlines():
        currLine = line.strip().split('	')
        lineArr =[]
        #一共有21个特征
        for i in range(21):
            #lineArr以数组形式保存由字符转换为浮点数的一条数据
            lineArr.append(float(currLine[i]))
        trainingSet.append(lineArr)
        trainingLabels.append(float(currLine[21]))
    trainWeights = stocGradAscent1(array(trainingSet), trainingLabels, 1000)
    errorCount = 0; numTestVec = 0.0
    #调用测试集进行测试
    for line in frTest.readlines():
        numTestVec += 1.0
        currLine = line.strip().split('	')
        lineArr =[]
        for i in range(21):
            lineArr.append(float(currLine[i]))
        if int(classifyVector(array(lineArr), trainWeights))!= int(currLine[21]):
            errorCount += 1
    errorRate = (float(errorCount)/numTestVec)#返回误差率
    print("the error rate of this test is: %f" % errorRate)
    return errorRate
#马类问题的测试接口
def multiTest():
    #重复执行10次以获取模型分类的平均误差
    numTests = 10; errorSum=0.0
    for k in range(numTests):
        errorSum += colicTest()
    print ("after %d iterations the average error rate is: %f" % (numTests, errorSum/float(numTests)))
        
dataArr,labelMat=loadDataSet()
weights=gradAscent(dataArr,labelMat)
print("梯度下降法：")
#getA()是numpy的一个函数，作用是将矩阵转成一个ndarray，getA()函数和mat()函数的功能相反，是将一个矩阵转化为数组。
plotBestFit(weights.getA())
print("随机梯度下降法")
weights2=stocGradAscent0(array(dataArr),labelMat)
plotBestFit(weights2)
print("一次随机梯度下降很可能是得不到最优结果的")
print("改进的随机梯度下降法")
weights3=stocGradAscent1(array(dataArr),labelMat)
plotBestFit(weights3)
print("马类问题")
multiTest()