机器学习实战笔记-Logistic回归

假设现在有一些数据点，我们用一条直线对这些点进行拟合（该线称为最佳拟合直线），这个拟合过程就称作回归。
利用Logistic回归进行分类的主要思想是：根据现有数据对分类边界线建立回归公式，以此进行分类。这里的“ 回归” 一词源于最佳拟合，表示要找到最佳拟合参数集
训练分类器时的做法就是寻找最佳拟合参数，使用的是最优化算法。

Logistic回归的一般过程
(1)收集数据：采用任意方法收集数据。
(2)准备数据：由于需要进行距离计算，因此要求数据类型为数值型。另外，结构化数据
格式则最佳。
(3)分析数据：采用任意方法对数据进行分析。
(4)训练算法：大部分时间将用于训练，训练的目的是为了找到最佳的分类回归系数。
(5)测试算法：一旦训练步驟完成，分类将会很快。
(6)使用算法：首先，我们需要输入一些数据，并将其转换成对应的结构化数值；
接着，基于训练好的回归系数就可以对这些数值进行简单的回归计算，判定它们属于
哪个类别.，在这之后，我们就可以夺输出的类别上做一些其他分析工作。

5.1基于Logistic回归和Sigmoid函数的分类
Logistic回归
优点：计算代价不高，易于理解和实现。
缺点：容易欠拟合，分类精度可能不高。 .
适用数据类型：数值型和标称型数据。

我们想要的函数应该是，能接受所有的输人然后预测出类别。例如 ，在两个类的情况下，上述函数输出0或 1。或许你之前接触过具有这种性质的函数，该函数称为海维塞德阶跃函数(Heaviside step function) ,或者直接称为单位阶跃函数。然而，海维塞德阶跃函数的问题在于：该函数在跳跃点上从0瞬间跳跃到1，这个瞬间跳跃过程有时很难处理。幸好，另一个函数也有类似的性质，且数学上更易处理，这就是Sigmoid函数。Sigmoid函数具体的计算公式如下：
f(z)=11+e−z

图5-1给出了Sigmoid函数在不同坐标尺度下的两条曲线图。当x为0时,Sigmoid函数值为0.5。随着x的增大，对应的Sigmoid值将逼近于1; 而随着x的减小，Sigmoid值将逼近于0。如果横坐标刻度足够大(图5-1下图),Sigmoid函数看起来很像一个阶跃函数。
因此，为了实现Logistic回归分类器，我们可以在每个特征上都乘以一个回归系数，然后把所有的结果值相加，将这个总和代人Sigmoid函数中，进而得到一个范围在0〜1之间的数值。任何大于0.5的数据被分人1类 ，小于0.5即被归人0类 。所 以,Logistic回归也可以被看成是一种概率估计。

5.2 基于最优化方法的最佳回归系数确定
Sigmoid函数的输人记为z，由下面公式得出:
z=w0x0+w1x1+w2x2+...wnxn
如果采用向量的写法，上述公式可以写成z=wTx, 它表示将这两个数值向量对应元素相乘然后全部加起来即得到z值。其中的向量x是分类器的输人数据，向量w就是我们要找到的最佳参数(系数) , 从而使得分类器尽可能地精确。为了寻找该最佳参数，需要用到最优化理论的一些知识。
下面首先介绍梯度上升的最优化方法，我们将学习到如何使用该方法求得数据集的最佳参数。接下来，展示如何绘制梯度上升法产生的决策边界图，该图能将梯度上升法的分类效果可视化地呈现出来。最后我们将学习随机梯度上升算法，以及如何对其进行修改以获得更好的结果。

5.2.1梯度上升法
梯度上升法基于的思想是：要找到某函数的最大值，最好的方法是沿着该函数的梯度方向探寻。如果梯度记为∇ , 则函数f(x,y)的梯度由下式表示:

图5 - 2 梯度上升算法到达每个点后都会重新估计移动的方向。从P0开始，计算完该点的梯度，函数就根据梯度移动到下一点P1。在P1点，梯度再次被重新计算，并沿新的梯度方向移动到P2。如此循环迭代，直到满足停止条件。迭代的过程中，梯度算子总是保证我们能选取到最佳的移动方向
图5-2中的梯度上升算法沿梯度方向移动了一步。可以看到，梯度算子总是指向函数值增长最快的方向。这里所说的是移动方向，而未提到移动量的大小。该量值称为步长，记 做 《。用向量来表示的话，梯度算法的迭代公式如下：
w=w+α∇wf(w)
该公式将一直被迭代执行，直至达到某个停止条件为止，比如迭代次数达到某个指定值或算法达到某个可以允许的误差范围。

梯度下降算法
公式:w=w−α∇wf(w)

梯度上升算法用来求函数的最大值，而梯度下降算法用来求函数的最小值。

重点内容
5.2.2 训练算法：使用梯度上升找到最佳参数
图5-3中有100个样本点，每个点包含两个数值型特征：X1和X2 。在此数据集上，我们将通过使用梯度上升法找到最佳回归系数，也就是拟合出logistic回归模型的最佳参数。
梯度上升法的伪代码如下：
每个回归系数初始化为1
重复R次：
计算整个数据集的梯度
使用alpha * gradient更新回归系数的向量
返回回归系数

Logistic回归梯度上升优化算法,代码:

from numpy import *
def loadDataSet():
    dataMat = []; labelMat = []
    fr = open('testSet.txt')
    for line in fr.readlines():
        lineArr = line.strip().split()
        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
        labelMat.append(int(lineArr[2]))
    return dataMat,labelMat

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

def gradAscent(dataMatIn, classLabels):
    #将二维矩阵转换为矩阵，之后有矩阵运算
    dataMatrix = mat(dataMatIn)
    #将labelMat转换为列向量
    labelMat = mat(classLabels).transpose()
    m,n = shape(dataMatrix)
    alpha = 0.001
    maxCycles = 500
    weights = ones((n,1))
    for k in range(maxCycles):
        #可以感性地将这两行代码理解为计算推测值与实际值的差距
        h = sigmoid(dataMatrix*weights)
        error = (labelMat - h)
        #其中的dataMatrix.transpose()*error对应梯度
        #在我转载的一篇博客中有详细的梯度公式推导
        weights = weights + alpha * dataMatrix.transpose()*error
    return weights

5.2.3 分析数据：画出决策边界
画出数据集和Logistic回归最佳拟合直线的函数,代码如下:

def plotBestFit(weights):
    import matplotlib.pyplot as plt
    dataMat,labelMat = loadDataSet()
    dataArr = array(dataMat)
    n = shape(dataArr)[0]
    #将矩阵转化为二维数组
    weights = weights.getA()
    #分别将类别为0,1的点存储在对应的数组中，方便之后绘点
    xcord1 = [];ycord1 = []
    xcord2 = [];ycord2 = []
    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')
    x = arange(-3.0,3.0,0.1)
    y = (-weights[0][0] - weights[1][0] * x)/weights[2][0]
    ax.plot(x,y)
    plt.xlabel('X1');plt.ylabel('X2')
    plt.show()

测试截图如下:

5.2.4 训练算法：随机梯度上升
梯度上升算法在每次更新回归系数时都需要遍历整个数据集, 该方法在处理100个左右的数据集时尚可，但如果有数十亿样本和成千上万的特征，那么该方法的计算复杂度就太高了。一种改进方法是一次仅用一个样本点来更新回归系数，该方法称为随机梯度上升算法。由于可以在新样本到来时对分类器进行增量式更新，因而随机梯度上升算法是一个在线学习算法。与 “在线学习”相对应，一次处理所有数据被称作是“批处理” 。
随机梯度上升算法可以写成如下的伪代码：
所有回归系数初始化为1
对数据集中每个样本
计算该样本的梯度
使用alpha∗gradient更新回归系数值
返回回归系数值

随机梯度上升算法 代码:

def stocGradAscent0(dataMatrix, classLabels):
    m,n = shape(dataMatrix)
    alpha = 0.01
    weights = ones(n)
    for i in range(m):
        h = sigmoid(sum(dataMatrix[i] * weights))
        error = classLabels[i] - h
        weights = weights + alpha * error * dataMatrix[i]
    return weights

测试截图:

可以看到，随机梯度上升算法与梯度上升算法在代码上很相似，但也有一些区别：第一，后者的变量h和误差error都是向量，而前者则全是数值；第二，前者没有矩阵的转换过程，所有变量的数据类型都是Numpy数组。
执行完毕后将得到图5-5所示的最佳拟合直线图，该图与图5-4有一些相似之处。可以看到，拟合出来的直线效果还不错，但并不像图5-4那样完美。这里的分类器错分了三分之一的样本。
直接比较程序清单5-3和程序清单5-1的代码结果是不公平的，后者的结果是在整个数据集上迭代了500次才得到的。一个判断优化算法优劣的可靠方法是看它是否收敛，也就是说参数是否达到了稳定值，是否还会不断地变化？对此 ，我们在程序清单5-3中随机梯度上升算法上做了些修改，使其在整个数据集上运行200次。最终绘制的三个回归系数的变化情况如图5-6所示。

图5-6展示了随机梯度上升算法在200次迭代过程中回归系数的变化情况。其中的系数2 ,也就是图5-5中的X2只经过了50次迭代就达到了稳定值，但系数1和0则需要更多次的迭代。另外值得注意的是，在大的波动停止后，还有一些小的周期性波动。不难理解，产生这种现象的原因是存在一些不能正确分类的样本点（数据集并非线性可分)，在每次迭代时会引发系数的剧烈改变。我们期望算法能避免来回波动，从而收敛到某个值。另外，收敛速度也需要加快。

改进的随机梯度上升算法,代码如下:

def stocGradAscent1(dataMatrix, classLabels, numIter = 150):
    m,n = shape(dataMatrix)
    weights = ones(n)
    for j in range(numIter):
        dataIndex = list(range(m))
        for i in range(m):
            alpha = 4/(1.0+j+i)+0.01
            randIndex = int(random.uniform(0, len(dataIndex)))
            h = sigmoid(sum(dataMatrix[randIndex] * weights))
            error = classLabels[randIndex] - h
            weights = weights + alpha * error * dataMatrix[randIndex]
            del(dataIndex[randIndex])
    return weights

测试截图如下:

书中的代码解读:
程序清单5-4与程序清单5-3类似，但增加了两处代码来进行改进。一方面,alpha在每次迭代的时候都会调整，这会缓解图5-6上的数据波动或者高频波动。另外，虽然alpha会随着迭代次数不断减小，但永远不会减小到0，这是因为alpha计算中还存在一个常数项。必须这样做的原因是为了保证在多次迭代之后新数据仍然具有一定的影响。如果要处理的问题是动态变化的，那么可以适当加大上述常数项，来确保新的值获得更大的回归系数。另一点值得注意的是，在降低alpha的函数中，alpha每次减少lj+i>，其中j是迭代次数，i是样本点的下标。这样,当j<<max(i)时，alpha就不是严格下降的。避免参数的严格下降也常见于模拟退火算法等其他优化算法中。
第二个改进的地方在,通过随机选取样本来更新回归系数。这种方法将减少周期性的波动（如图5-6中的波动)。具体实现方法与第3章类似，这种方法每次随机从列表中选出一个值，然后从列表中删掉该值（再进行下一次迭代)。此外，改进算法还增加了一个迭代次数作为第3个参数。如果该参数没有给定的话，算法将默认迭代150次。如果给定，那么算法将按照新的参数值进行迭代。

5.3 示例：从疝气病症预测病马的死亡率
示例:使用Logistic回归估计马疝病的死亡率
(1)收集数据：给定数据文件。
(2)准备数据：用Python解析文本文件并填充缺失值。
(3)分析数据：可视化并观察数据。
(4)训练算法：使用优化算法，找到最佳的系数。
(5)测试算法：为了量化回归的敢果，需要观察错误率。根据错误率决定是否回退到训练
阶段，通过改变迭代的次数和步长等参数来得到更好的回归系数。
(6)使用算法：实现一个简单的命令行程序来收集马的症状并输出预测结果并非难事，这
可以做为留给读者的一道习题。

5.3.1 准备数据：处理被据中的缺失值
数据中的缺失值是个非常棘手的问题，有很多文献都致力于解决这个问题。那么，数据缺失究竟带来了什么问题？假设有100个样本和20个特征，这些数据都是机器收集回来的。若机器上的某个传感器损坏导致一个特征无效时该怎么办？此时是否要扔掉整个数据？这种情况下，另外19个特征怎么办？它们是否还可用？答案是肯定的。因为有时候数据相当昂贵，扔掉和重新获取都是不可取的，所以必须采用一些方法来解决这个问题。
下面给出了一些可选的做法：
□使用可用特征的均值来填补缺失值；
□使用特殊值来±真补缺失值，如-1;
□忽略有缺失值的样本；
□使用相似样本的均值添补缺失值；
□使用另外的机器学习算法预测缺失值。

现在，我们对下一节要用的数据集进行预处理，使其可以顺利地使用分类算法。在预处理阶段需要做两件事：第一，所有的缺失值必须用一个实数值来替换，因为我们使用的 Numpy数据类型不允许包含缺失值。这里选择实数0来替换所有缺失值，恰好能适用于Logistic回归。这样做的直觉在于，我们需要的是一个在更新时不会影响系数的值。回归系数的更新公式如下：
weights=weights+alpha∗error∗dataMatrix[randIndex]
如果dataMatrix某特征对应值为0，那么该特征的系数将不做更新，即：
w e i g h t s = w e i g h t s
另外，由于sigmoid(0)=12 ，即它对结果的预测不具有任何倾向性，因此上述做法也不会对误差项造成任何影响。基于上述原因，将缺失值用0代替既可以保留现有数据，也不需要对优化算法进行修改。此外，该数据集中的特征取值一般不为0 , 因此在某种意义上说它也满足“特殊值”这个要求。

预处理中做的第二件事是，如果在测试数据集中发现了一条数据的类别标签已经缺失，那么我们的简单做法是将该条数据丢弃。这是因为类别标签与特征不同，很难确定采用某个合适的值来替换。采用Logistic回归进行分类时这种做法是合理的，而如果采用类似KNN的方法就可能不太可行

5.3.2 测试算法：用Logistic回归进行分类
代码:

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 = []
        for i in range(21):
            lineArr.append(float(currLine[i]))
        trainingSet.append(lineArr)
        trainingLabels.append(float(currLine[21]))
    trainWeights = stocGradAscent1(array(trainingSet),trainingLabels,500)
    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 rror rate of thits test is: %f" % errorRate)
    return errorRate

def multiTest():
    numTests = 10;errorSum=0.0
    for k in range(numTests):
        errorSum += colicTest()
    print("after %d iterations the average err

测试截图: