机器学习实践之Logistic回归

Posted 2020-10-18 蓝色之旅

    关于本文说明，本人原博客地址位于http://blog.csdn.net/qq_37608890，本文来自笔者于2017年12月17日 19:18:31所撰写内容（http://blog.csdn.net/qq_37608890/article/details/78827013）。

 

本文根据最近学习机器学习书籍 网络文章的情况,特将一些学习思路做了归纳整理,详情如下.如有不当之处,请各位大拿多多指点,在此谢过。

         （今天发现第二部分 第4点中，部分代码不整齐，重新梳理了。2017.12.21记）

   本文将开展最优化算法的学习和梳理。在日常生活中，我们遇到过很多最优化的问题，例如你打算明天从北京出发去青岛，如何在最短的时间里到达青岛？如何投入最少的工作量而获得最大的报酬？等等。接下来，将要介绍一些最优化算法并利用它们训练出一个非线性函数用于分类。

一 概述

1  概念理解  

      假设现在有一些数据点，我们使用一条直线对这些点进行拟合（该线被称为最佳拟合直线），这个拟合过程称为回归。

       其实使用Logistic回归进行分类的主要思想是：根据现有数据对分类边界线建立回归公式，以此进行分类。

        其实“回归”一词来源于最佳拟合，表示找到最佳拟合参数。训练分类器的做法：寻找最佳拟合参数，使用的是最优化算法（梯度上升法、改进的随机梯度上升法）。

2   基于Logistic回归和Sigmoid函数的分类

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

            

      Sigmoid函数也可表示为hΘ(X)=g(ΘTX).

       下图给出了 Sigmoid 函数在不同坐标尺度下的两条曲线图。当 x 为 0 时，Sigmoid 函数值为 0.5 。随着 x 的增大，对应的 Sigmoid 值将逼近于 1 ; 而随着 x 的减小， Sigmoid 值将逼近于 0 。如果横坐标刻度足够大， Sigmoid 函数看起来很像一个阶跃函数。
      

 

  

 

    为了实现Logistic回归分类器，我们需要在每个特征上乘以一个回归系数，然后把所有结果值相加，将这个总和代入Sigmoid函数，进而得到0~1之间的数值。

   

  

  现在就可以对标签Y进行分类：

   

Cost function:

    线性回归的Cost Function依据最小二乘法使得观察值和估计值差的平方和最小，即：

   

  但这里对于logistic回归而言，我们却不能最小化观察值和估计值的差的平方和，因为这样就会发现J(θ)为非凸函数，就存在多个局部极值点，从而不能开展梯度迭代得到最优参数，所以，这里我们重新定义一种Cost Function。

 

通过以上两个函数的函数曲线，我们会发现当y=1，而估计值h=1或者当y=0，而估计值h=0,即预测准确了，此时的cost就为0，，但是当预测错误了cost就会无穷大，很明显满足cost function的定义。

可以将上面的分组函数写在一起：

从而总体损失函数（Cost Function）为：

  

为什么要采用上面的函数作为cost function？

       Andrew Ng给的解释：因为最小估计值和观察值的差平方和为非凸函数，通过函数曲线观察得到上面的cost function满足条件。

        还可以给出另外一种解释，那就是最大似然估计:

            由上面可以知道：hθ(x)≥0.5，此时y=1， 小于0.5，y=0。 那么我们就用hθ作为y=1发生的概率，那么当y=0时，hθ＜0.5，此时不能用hθ作为y=0的概率，否则由于最大似然的思想使已有的数据发生的概率最大化，小于0.5太小了，但我们可以用1-hθ作为y=0的概率，这样以来就可以作为y=0的概率了，接下来我们只需要最大化联合概率密度函数就可以了。

       联合概率密度函数可以写成：

    

 

     为了和上面的似然函数保持一致，这里再转换成对数似然函数：

         

 

3   基于最优化方法的回归系数确定

    Sigmoid 函数的输入记为 z ，由下面公式得到:

     

     如果采用向量的写法，上述公式可以写成z=w^Tx,  ，它表示将这两个数值向量对应元素相乘然后全部加起来即得到 z 值。其中的向量 x 是分类器的输入数据，向量 w 也就是我们要找到的最佳参数（系数），从而使得分类器尽可能地精确。为了寻找该最佳参数，需要用到最优化理论的一些知识。我们这里使用的是——梯度上升法。
        

3.1 梯度的介绍

         向量 = 值 + 方向
         梯度 = 向量
         梯度 = 梯度值 + 梯度方

3.2 梯度上升法的思想

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

    

 

    这个梯度意味着要沿 x 的方向移动 ，沿 y 的方向移动 。其中，函数f(x, y) 必须要在待计算的点上有定义并且可微。下图是一个具体的例子。

       

            上图展示的，梯度上升算法到达每个点后都会重新估计移动的方向。从 P0 开始，计算完该点的梯度，函数就根据梯度移动到下一点 P1。在 P1 点，梯度再次被重新计算，并沿着新的梯度方向移动到 P2 。如此循环迭代，直到满足停止条件。迭代过程中，梯度算子总是保证我们能选取到最佳的移动方向。

           上图中的梯度上升算法沿梯度方向移动了一步。可以看到，梯度算子总是指向函数值增长最快的方向。这里所说的是移动方向，而未提到移动量的大小。该量值称为步长，记作 α 。用向量来表示的话，梯度上升算法的迭代公式如下:

            

          该公式将一直被迭代执行，直至达到某个停止条件为止，比如迭代次数达到某个指定值或者算法达到某个可以允许的误差范围。

          介绍一下几个相关的概念：

           例如：y = w1x1 + w2x2 + ... + wnxn
         梯度：参考上图的例子，二维图像，x方向代表第一个系数，也就是 w1，y方向代表第二个系数也就是 w2，这样的向量就是梯度。
           α：上面的梯度算法的迭代公式中的阿尔法，这个代表的是移动步长。移动步长会影响最终结果的拟合程度，最好的方法就是随着迭代次数更改移动步长。
          步长通俗的理解，100米，如果我一步走10米，我需要走10步；如果一步走20米，我只需要走5步。这里的一步走多少米就是步长的意思。
          ▽f(w)：代表沿着梯度变化的方向。

             

       当然对于随机化的梯度迭代每次只使用一个样本进行参数更新，则如下所示（也是下面代码中公式的来源）：

        

         例如：data=[1,2,3;4,5,6;7,8,9;10,11,12]为4个样本点，3个特征的数据集，，此时标签为[1,0,0,0]，

       那么用梯度上升 表达的就是当j=0时，就是第一列[1,4,7,10]与标签差的乘积。

 

3.3 梯度下降法简要介绍

        平时常听到的是梯度下降算法，它与这里的梯度上升算法是一样的，只是公式中的加法需要变成减法。因此，对应的公式可以写成

       

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

3.4 关于局部最优现象

       如下图所示，参数 θ 与误差函数 J(θ) 的关系图，红色的部分是表示 J(θ) 有着比较高的取值，我们需要的是，能够让 J(θ) 的值尽量的低。也就是深蓝色的部分。θ0，θ1 表示 θ 向量的两个维度。

可能梯度下降的最终点并非是全局最小点，可能是一个局部最小点，如我们上图中的右边的梯度下降曲线，描述的是最终到达一个局部最小点，这是我们重新选择了一个初始点得到的。

看来我们这个算法将会在很大的程度上被初始点的选择影响而陷入局部最小点。

 

    

 

4 Logistic回归的一般过程

    (1)收集数据: 可以使用任何方法。

    (2)准备数据: 由于需要进行距离计算，因此要求数据类型为数值型。当然，结构化数据格式最佳。

    (3)分析数据: 采用任意方法对数据进行分析。

    (4)训练算法: 大部分时间将用于训练，训练的目的是为了找到最佳的分类回归系数。

    (5)测试算法: 一旦训练阶段完成，分类将很快完成。

    (6)使用算法: 首先，需要输入一些数据，并将其转换成对应的结构化数值；然后，基于训练好的回归系数就可以对这些数值进行简单的回归计算，判定它们属于哪个类别；最后，可以在输出的类别上做一些其他分析工作。

 

5 相关特性

      优点：计算代价不高，易于理解和实现。

      缺点： 容易欠拟合，分类精度可能不高。

      适用数据类型： 数值型和标称型。

 

二 项目案例1: 使用 Logistic 回归在简单数据集上的分类

1   项目概述

         在一个简单的数据集上，采用梯度上升法找到 Logistic 回归分类器在此数据集上的最佳回归系数。

2  Logistic回归工作原理

           每个回归系数初始化为 1
           重复 R 次:
                    计算整个数据集的梯度
                    使用 步长 x 梯度 更新回归系数的向量
         返回回归系数   

3  简单数据集estSet.txt中数据存储格式

         -0.346811 -1.678730 1
          -2.124484 2.672471 1
           1.217916 9.597015 0
          -0.733928 9.098687 0
          -3.642001 -1.618087 1
            0.315985 3.523953 1

 

4 训练样本

          100个样本点，每个点包含两个数值型特征：x1和x2。

    

    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])])#为方便计算，将x0值设为1.0  
            labelMat.append(int(lineArr[2]))  
        return dataMat, labelMat  
      
    def sigmoid(inX):  
        return 1.0/(1+exp(-inX))  
      
    def gradAscent(dataMatIn, classLabels):#梯度上升：dataMatIn:2维NumPy数组,100*3矩阵;classLabels:类别标签，1*100行向量  
        dataMatrix = mat(dataMatIn)#特征矩阵  
        labelMat = mat(classLabels).transpose()#类标签矩阵：100*1列向量  
        m,n = shape(dataMatrix)  
        alpha = 0.001#向目标移动的步长  
        maxCycles = 500#迭代次数  
        weights = ones((n,1))#n*1列向量：3行1列  
        for k in range(maxCycles):  
            h = sigmoid(dataMatrix*weights)#100*3*3*1=100*1,dataMatrix * weights代表不止一次乘积计算，事实上包含了300次乘积  
            error = (labelMat - h)#真实类别与预测类别的差值  
            weights = weights + alpha * dataMatrix.transpose()* error#w:=w+α▽wf(w)  
        return weights  
      
    #画出数据集和Logistic回归最佳拟合直线的函数  
    def plotBestFit(weights):  
        import matplotlib.pyplot as plt  
        dataMat, labelMat = loadDataSet()  
        dataArr = array(dataMat)  
        n = shape(dataArr)[0]#n=100  
        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)#arange创建等差数组,-3.0起始点,3.0终止点(不包含3.0),间隔为0.1  
        y = (-weights[0] - weights[1] * x)/weights[2]#最佳拟合直线,设置sigmoid函数为0,0是两个分类(类别1和类别0)的分界处。因此设定0=w0x0+w1x1+w2x2,解出x1和x2关系(即分割线的方程，x0=1)。  
        ax.plot(x, y)  
        plt.xlabel(‘X1‘);plt.ylabel(‘X2‘);  
        plt.show()  

 

 执行如下代码：

   

dataArr,labelMat = loadDataSet()weights =gradAscent(dataArr,labelMat)plotBestFit(weights.getA())

  程序运行之后能看到类似于下图的结果图：

 

         其实，这个分类结果是很不错，虽然这里是简单且数据集很小，但是这个方法却需要大量的计算（300次乘法）。

         下面会对该算法稍作改进，从而使它可以用在其他真是数据上。

5 训练算法

       梯度上升算法：在每次更新回归系数时都需要遍历整个数据集，在数十亿样本上该算法复杂度太高。

       改进方法：随机梯度上升算法：一次仅用一个样本点更新回归系数。

       由于可以在新样本到来时对分类器进行增量式更新，因此随机梯度上升算法是一个在线学习算法。与“在线学习”相对应，一次处理所有数据被称作“批处理”。

 

    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  

 

 执行 

    dataArr,labelMat = loadDataSet()  
    weights =stocGradAscent0(array(dataArr),labelMat)  
    plotBestFit(weights)  

可以得到

  

   

           拟合效果没有梯度上升算法完美。这里的分类器错分了三分之一的样本。梯度上升算法的结果是在整个数据集上迭代了500次才得到的。对此，在程序中随机梯度上升做些修改，使其在整个数据集上运行200次。最终绘制的三个回归系数变化情况如下图：

判断优化算法优劣的可靠方法：看它是否收敛，也就是说参数是否达到稳定值，是否不断地变化。

     

        可以看到，随机梯度上升算法与梯度上升算法在代码上很相似，但也有一些区别: 第一，后者的变量 h 和误差 error 都是向量，而前者则全是数值；第二，前者没有矩阵的转换过程，所有变量的数据类型都是 NumPy 数组。

        判断优化算法优劣的可靠方法是看它是否收敛，也就是说参数是否达到了稳定值，是否还会不断地变化？下图展示了随机梯度上升算法在 200 次迭代过程中回归系数的变化情况。其中的系数2，也就是 X2 只经过了 50 次迭代就达到了稳定值，但系数 1 和 0 则需要更多次的迭代。如下图所示:

     

 

针对这个问题，我们改进了之前的随机梯度上升算法，如下：

 

def stocGradAscent1(dataMatrix, classLabels, numIter=150):  
    m,n = shape(dataMatrix)  
    weights = ones(n)     
      
    for j in range(numIter):  
        # [0, 1, 2 .. m-1]  
        dataIndex = list(range(m))  
        for i in range(m):  
            
            alpha = 4/(1.0+j+i)+0.0001     
             
            randIndex = int(random.uniform(0,len(dataIndex)))  
            # sum(dataMatrix[i]*weights)为了求 f(x)的值， f(x)=a1*x1+b2*x2+..+nn*xn  
            h = sigmoid(sum(dataMatrix[randIndex]*weights))  
            error = classLabels[randIndex] - h  
            # print weights, ‘__h=%s‘ % h, ‘__‘*20, alpha, ‘__‘*20, error, ‘__‘*20, dataMatrix[randIndex]  
            weights = weights + alpha * error * dataMatrix[randIndex]  
            del(dataIndex[randIndex])  
    return weights  

 

 

上面的改进版随机梯度上升算法，我们修改了两处代码。

第一处改进为 alpha 的值。alpha 在每次迭代的时候都会调整，这回缓解上面波动图的数据波动或者高频波动。另外，虽然 alpha 会随着迭代次数不断减少，但永远不会减小到 0，因为我们在计算公式中添加了一个常数项。

第二处修改为 randIndex 更新，这里通过随机选取样本拉来更新回归系数。这种方法将减少周期性的波动。这种方法每次随机从列表中选出一个值，然后从列表中删掉该值（再进行下一次迭代）。

程序运行之后能看到类似于下图的结果图。

         

        该方法比采用固定的alpha收敛速度更快。主要归功于：1.stocGradAscent1()的样本随机机制避免周期性波动；2.stocGradAscent1()收敛更快。这次仅对数据集做了20次遍历，而之前的方法是500次。

    

 

 

三  项目案例2: 从疝气病症预测病马的死亡率

1 项目概述

      使用 Logistic 回归来预测患有疝病的马的存活问题。疝病是描述马胃肠痛的术语。然而，这种病不一定源自马的胃肠问题，其他问题也可能引发马疝病。这个数据集中包含了医院检测马疝病的一些指标，有的指标比较主观，有的指标难以测量，例如马的疼痛级别。

2   病马的训练数据已经给出来了，如下形式存储在文本文件中(这里打数据集样式重新选择打，不影响模型设计和校验，与本人在csdn博客上有些不同，但都是来源于同一数据集，即horseColicTraining.txt）

 

2.000000	9.000000	38.300000	90.000000	0.000000	1.000000	0.000000	1.000000	1.000000	5.000000	3.000000	1.000000	2.000000	1.000000	0.000000	3.000000	0.000000	40.000000	6.200000	1.000000	2.200000	1.000000
1.000000	1.000000	38.100000	66.000000	12.000000	3.000000	3.000000	5.000000	1.000000	3.000000	3.000000	1.000000	2.000000	1.000000	3.000000	2.000000	5.000000	44.000000	6.000000	2.000000	3.600000	1.000000
2.000000	1.000000	39.100000	72.000000	52.000000	2.000000	0.000000	2.000000	1.000000	2.000000	1.000000	2.000000	1.000000	1.000000	0.000000	4.000000	4.000000	50.000000	7.800000	0.000000	0.000000	1.000000
1.000000	1.000000	37.200000	42.000000	12.000000	2.000000	1.000000	1.000000	1.000000	3.000000	3.000000	3.000000	3.000000	1.000000	0.000000	4.000000	5.000000	0.000000	7.000000	0.000000	0.000000	1.000000
2.000000	9.000000	38.000000	92.000000	28.000000	1.000000	1.000000	2.000000	1.000000	1.000000	3.000000	2.000000	3.000000	0.000000	7.200000	1.000000	1.000000	37.000000	6.100000	1.000000	0.000000	0.000000
1.000000	1.000000	38.200000	76.000000	28.000000	3.000000	1.000000	1.000000	1.000000	3.000000	4.000000	1.000000	2.000000	2.000000	0.000000	4.000000	4.000000	46.000000	81.000000	1.000000	2.000000	1.000000
1.000000	1.000000	37.600000	96.000000	48.000000	3.000000	1.000000	4.000000	1.000000	5.000000	3.000000	3.000000	2.000000	3.000000	4.500000	4.000000	0.000000	45.000000	6.800000	0.000000	0.000000	0.000000
1.000000	9.000000	0.000000	128.000000	36.000000	3.000000	3.000000	4.000000	2.000000	4.000000	4.000000	3.000000	3.000000	0.000000	0.000000	4.000000	5.000000	53.000000	7.800000	3.000000	4.700000	0.000000
2.000000	1.000000	37.500000	48.000000	24.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	1.000000
1.000000	1.000000	37.600000	64.000000	21.000000	1.000000	1.000000	2.000000	1.000000	2.000000	3.000000	1.000000

  

3   处理数据中的缺失值

         假设有100个样本和20个特征，这些数据均由机器收集获得。如果机器上的某个传感器损坏引起一个特征无效，此时我们该怎么处理？是否要扔掉整个数据？这种情况下，另外19个特征怎么办？ 它们是否还可以用？答案是肯定的。因为有时候数据相当昂贵，扔掉和重新获取都是不可取的，所以必须采用一些方法来解决这个问题。

         下面给出了一些可选的做法：                  

• 使用可用特征的均值来填补缺失值；

• 使用特殊值来填补缺失值，如 -1；

• 忽略有缺失值的样本；

• 使用有相似样本的均值添补缺失值；

• 使用另外的机器学习算法预测缺失值。

       现在，我们对下一节要用的数据集进行预处理，使其可以顺利地使用分类算法。在预处理需要做两件事:

   （1）    所有的缺失值必须用一个实数值来替换，因为我们使用的 NumPy 数据类型不允许包含缺失值。我们这里选择实数 0 来替换所有缺失值，恰好能适用于 Logistic 回归。这样做的直觉在于，我们需要的是一个在更新时不会影响系数的值。回归系数的更新公式如下:

                   weights = weights + alpha * error * dataMatrix[randIndex]

                   如果 dataMatrix 的某个特征对应值为 0，那么该特征的系数将不做更新，即:

                  weights = weights

                  另外，由于 Sigmoid(0) = 0.5 ，即它对结果的预测不具有任何倾向性，因此我们上述做法也不会对误差造成任何影响。基于上述原因，将缺失值用 0 代替既可以保留现有数据，也不需要对优化算法进行修改。此外，该数据集中的特征取值一般不为 0，因此在某种意义上说它也满足 “特殊值” 这个要求。

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

原始的数据集经过预处理后，保存成两个文件: horseColicTest.txt 和 horseColicTraining.txt 。

4 分析数据

          可视化并观察数据

     将数据使用 MatPlotlib 打印出来，观察数据是否是我们想要的格式。

5  训练算法

          使用优化算法，找到最佳的系数

      原始的梯度上升算法，随机梯度上升算法，改进版随机梯度上升算法的代码如下:

 

# 正常的处理方案  
# 两个参数：第一个参数==> dataMatIn 是一个2维NumPy数组，每列分别代表每个不同的特征，每行则代表每个训练样本。  
# 第二个参数==> classLabels 是类别标签，它是一个 1*100 的行向量。为了便于矩阵计算，需要将该行向量转换为列向量，做法是将原向量转置，再将它赋值给labelMat。  
def gradAscent(dataMatIn, classLabels):  
    dataMatrix = mat(dataMatIn)             # 转换为 NumPy 矩阵  
    labelMat = mat(classLabels).transpose() # 首先将数组转换为 NumPy 矩阵，然后再将行向量转置为列向量  
    m,n = shape(dataMatrix)  
    alpha = 0.001  
    maxCycles = 500  
    weights = ones((n,1))  
    for k in range(maxCycles):              #heavy on matrix operations  
        h = sigmoid(dataMatrix*weights)     # 矩阵乘法  
        error = (labelMat - h)              # 向量相减  
        weights = weights + alpha * dataMatrix.transpose() * error # 矩阵乘法，最后得到回归系数  
    return array(weights)  
  
  
# 随机梯度上升  
# 梯度上升优化算法在每次更新数据集时都需要遍历整个数据集，计算复杂都较高  
# 随机梯度上升一次只用一个样本点来更新回归系数  
def stocGradAscent0(dataMatrix, classLabels):  
    m,n = shape(dataMatrix)  
    alpha = 0.01  
    weights = ones(n)   # 初始化长度为n的数组，元素全部为 1  
    for i in range(m):  
        h = sigmoid(sum(dataMatrix[i]*weights))  
        error = classLabels[i] - h  
        print(weights, "*"*10 , dataMatrix[i], "*"*10 , error)   
        weights = weights + alpha * error * dataMatrix[i]  
    return weights  
  
  
# 随机梯度上升算法（随机化）  
def stocGradAscent1(dataMatrix, classLabels, numIter=150):  
    m,n = shape(dataMatrix)  
    weights = ones(n)   # 创建与列数相同的矩阵的系数矩阵，所有的元素都是1  
    for j in range(numIter):  
        dataIndex = range(m)  
        for i in range(m):  
            alpha = 4/(1.0+j+i)+0.0001    # alpha 会随着迭代不断减小，但永远不会减小到0，因为后边还有一个常数项0.0001  
            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 

 

6  测算算法

     为了量化回归的效果，需要观察错误率。根据错误率决定是否回退到训练阶段，通过改变迭代的次数和步长的参数来得到更好的回归系数。

        

    # 分类函数，根据回归系数和特征向量来计算 Sigmoid的值  
    def classifyVector(inX, weights):  
        prob = sigmoid(sum(inX*weights))  
        if prob > 0.5: return 1.0  
        else: return 0.0  
      
    # 打开测试集和训练集,并对数据进行格式化处理  
    def colicTest():  
        frTrain = open(‘input/5.Logistic/horseColicTraining.txt‘)  
        frTest = open(‘input/5.Logistic/horseColicTest.txt‘)  
        trainingSet = []  
        trainingLabels = []  
        for line in frTrain.readlines():  
            currLine = line.strip().split(‘\t‘)  
            lineArr = []  
            for i in range(21):  
                lineArr.append(float(currLine[i]))  
            trainingSet.append(lineArr)  
            trainingLabels.append(float(currLine[21]))  
        # 使用 改进后的 随机梯度下降算法 求得在此数据集上的最佳回归系数 trainWeights  
        trainWeights = stocGradAscent1(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(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  
      
      
    # 调用 colicTest() 10次并求结果的平均值  
    def multiTest():  
        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)) )   

 

 执行  

multiTest()  

得到 

   

    /home/ibm_x/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:32: RuntimeWarning: overflow encountered in exp  
      
    the error rate of this test is: 0.358209  
    the error rate of this test is: 0.328358  
    the error rate of this test is: 0.402985  
    the error rate of this test is: 0.298507  
    the error rate of this test is: 0.328358  
    the error rate of this test is: 0.343284  
    the error rate of this test is: 0.417910  
    the error rate of this test is: 0.328358  
    the error rate of this test is: 0.343284  
    the error rate of this test is: 0.388060  
    after 10 iterations the average error rate is: 0.353731  

 

 
  7  结果分析

       从上面的结果来看，经过10次迭代之后的平均错误率为35%.实事求是地讲，这个结果并不差，因为有30%的数据缺失。当然，如果调整colicTest()中的迭代次数和stocGradAscent（）中的步长，平均错误率可以下降到20%左右。以后有机会再做进一步的优化。

 

  四  小结

          Logistic回归的目的应该是找到一个非线性函数sigmoid的最佳拟合参数，求解过程可以由最优化算法来完成。在使用最优化算法过程中，最常用的是梯度上升算法，而它又可以简化为随机梯度上升算法。

        随机梯度上升算法相比于梯度上升算法，从实际效果来看其实二者相当，但随机梯度上升算法占用的计算资源更少，且它是一个在线算法，可以在新数据到来时就完成参数更新，这就不需要重新读取整个数据集来进行批处理运算。

        机器学习的一个重要问题就是如何处理缺失数据。当然，这个问题没有标准答案，一切取决于实际应用中的需求。现成的一些解决方案，都有自己的优缺点，选择哪一种，具体问题具体分析。

         本文引出了回归分析的知识点，后续有机会的话，想单独梳理下回归分析背后的各种原理性的理论知识。

       如有不足之处，望各位网友多提宝贵意见，谢谢。