Python《机器学习实战》读书笔记（五）——Logistic回归

第五章 Logistic回归

前言

本章内容

• Sigmoid函数和Logistic回归分类器
• 最优化理论初步
• 梯度下降最优化算法
• 数据中的缺失项处理

假设现在有一些数据点，我们用一条直线对这些点进行拟合（该线称为最佳拟合直线），这个拟合过程就称为回归。利用Logistic回归进行分类的主要思想是：根据现有的数据对分类边界线建立回归公式，以此进行分类。

Logistic回归的一般过程：

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

5-1 基于Logistic回归和Sigmoid函数的分类

Logistic回归

优点：计算代价不高，易于理解和实现。
缺点：容易欠拟合，分类精度可能不高。
适用数据类型：数值型和标成型数据。

想得到的函数是能接受所有的输入然后预测出类别。例如：在两个类的情况下，该函数输出0或1。Sigmoid函数具有能扣输出0或1的这种性质。其函数具体的计算公式如下：

$$σ(z) = \dfrac{1}{1+e^{-z}}$$

当z=0时，Sigmoid函数值为0.5。随着z的增大，对应的Sigmoid值将逼近于1；随着z的减小，Sigmoid值将逼近于0。如下图所示是在不同坐标尺度下的两条曲线图：

从图中可以看出如果横坐标足够大，Sigmoid函数看起来很像一个阶跃函数。

为了实现Logistic回归分类器，可以在每个特征上都乘以一个回归系数，然后把所有的结果值相加，将这个总和带入Sigmoid函数中，进而得到一个范围在0~1之间的数值。任何大于0.5的数据被分为1类，小于0.5即被归入0类。所以Logistic回归也可以被看成一种概率估计。

5-2 基于最优化方法的最佳回归系数确定

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

$$z = ω_0x_0 + ω_1x_1 + ω_2x_2 + \cdots + ω_nx_n$$

如果采用向量的写法，上述公式可以写成：

$$z = w^Tx$$

它表示将两个数值向量对应元素相乘然后全部加起来即得到z值。

结合前面的Sigmoid函数：

$$σ(z) = \dfrac{1}{1+e^{-z}}$$

可以得到如下公式：

$$σ(w^Tx) = \dfrac{1}{1+e^{-w^Tx}}$$

其中向量 $x$ 是分类器的输入数据(样本列向量)，向量 $w$ 是要找到的最佳参数列向量(需要求解的系数)， $σ(w^Tx)$ 是实现了任意实数到 $[0,1]$ 的映射，这样，数据集 $[x_0,x_2,\cdots,x_n]$ 不管大于1还是小于0，都可以映射到 $[0,1]$ 区间进行分类。 $σ(w^Tx)$ 实际输出的是给定数据（给定样本）为1的概率，若概率大于0.5，则该样本分类为1，反之，样本分类为0，前提是这类分类属于两类分类问题。

如果我们找到合适的参数列向量$w(w_0,w_1,\cdots,w_n)$,同时，一直样本列向量$x[x_0,x_2,\cdots,x_n]$，那么我们就可以利用上面的公式计算出$σ(w^Tx)$概率，从而对样本X进行分类了。

不难看出公式中$x$我们容易得到，但怎么找到合适的参数向量$w$呢？

根据Sigmoid函数的特性，我们可以得到如下的公式：

$$P(y = 1|x;w) = σ_w(x)$$
$$P(y = 0|x;w) = 1 - σ_w(x)$$

上式即为在已知样本 $x$ 和参数 $w$ 的情况下，样本$x$属于类别1($y = 1$)和类别0($y = 0$)的条件概率。

将上述两个公式合二为一：

$$Cost(σ_w(x),y) = σ_w(x)^y(1-σ_w(x))^{1-y}$$

这儿的$Cost$函数成为代价函数，当 $y = 1$ 时，函数为

$$Cost(σ_w(x),1) = σ_w(x)$$

当当 $y = 0$ 时，函数为

$$Cost(σ_w(x),0) = (1-σ_w(x))$$

为了简化计算，对$Cost(σ_w(x),y)= σ_w(x)^y(1-σ_w(x))^{1-y}$ 两边求对数得到：

$$Cost(σ_w(x),y) = y logσ_w(x)+(1-y)log(1-σ_w(x))$$

上式的代价函数，如果给定一个样本，则可以通过它求出样本所属类别的概率，概率越大越好，所以便是求解这个代价函数的最大值。

输出 $y$ 发生的概率与回归参数$w$ 有关，因此我们可以对 $w$ 进行最大似然估计，使得y发生的概率最大，此时的$w$便是最优的回归系数

假设给定多个样本，样本与样本之间相互独立，则整个样本集生成的概率即为所有样本生成概率的乘积，将得到的公式对数化，将乘积转换成相加：

$$J(w) =\sum_{i=1}^n{y^{(i)} logσ_w(x^{(i)})+(1-y^{(i)})log(1-σ_w(x^{(i)}))}$$

其中，$n$ 为样本的总数，$y^{(i)}$为第$i$ 个样本的类别，$x^{(i)}$为第$i$个样本。$w$和$x^{(i)}$均为多维向量。

现在要求满足$J(w)$最大的$w$值，这儿使用梯度上升法进行求解。

5-2-1 梯度上升法

梯度上升法基于的思想是：要找到某函数的最大值，最好的方法是沿着该函数的梯度方向探寻。

函数$f(x,y)$的梯度由下式表示：

$$\nabla{f(x,y)} = \binom{\dfrac{\partial{f(x,y)}}{\partial{x}}}{\dfrac{\partial{f(x,y)}}{\partial{y}}}$$

这个梯度意味着沿$x$的方向移动$\dfrac{\partial{f(x,y)}}{\partial{x}}$,沿$y$的方向移动$\dfrac{\partial{f(x,y)}}{\partial{y}}$。

其中，$f(x,y)$必须要在待计算的点上有定义并且可微。

如下图所示，是一个具体的函数例子。

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

如上图所示，梯度算子总是指向函数增值增长最快的方向。这里所说的是移动方向，并没有提到移动量的大小。该量值称为步长，记做$α$。用向量来表示的话，梯度算法的迭代公式如下：

$$w:= w+α\nabla_wf(w)$$

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

梯度下降算法对应的公式如下：

$$w:= w-α\nabla_wf(w)$$

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

这儿的$\nabla_wf(w)$对应$\nabla_wJ(w)$，求导过程如下：

$J(w)$函数为$J(w) =\sum_{i=1}^n{y^{(i)} logσ_w(x^{(i)})+(1-y^{(i)})log(1-σ_w(x^{(i)}))}$

$σ_w(x)$函数为$σ_w(x) = g(w^Tx) =\dfrac{1}{1+e^{-w^Tx}}$

只要求出$J(w)$的偏导数，就可以利用梯度上升法求出$J(w)$极大值对应的$w$值。

$$\dfrac{\partial}{w_j}J(w) = \dfrac{\partial{J(w)}}{\partial{g(w^Tx)}}{*}{\dfrac{\partial{g(w^Tx)}}{\partial{w^Tx}}}{*}{\dfrac{\partial{w^Tx}}{\partial{w_j}}}$$

其中

$$\dfrac{\partial{J(w)}}{\partial{g(w^Tx)}} = y*\dfrac{1}{g(w^Tx)} + (y-1)*\dfrac{1}{1-g(w^Tx)}$$

$$\dfrac{\partial{g(w^Tx)}}{\partial{w^Tx}} = \dfrac{d}{dz}{\dfrac{1}{1+e^{-z}}} = \dfrac{1}{(1+e^{-z})^2}(e^{-z}) = \dfrac{1}{1+e^{-z}}(1-\dfrac{1}{1+e^{-z}}) = σ(z)(1-σ(z))$$

则

$$\dfrac{\partial{g(w^Tx)}}{\partial{w^Tx}} = σ(w^Tx)(1-σ(w^Tx))$$

$$\dfrac{\partial{w^Tx}}{w_j} = \dfrac{\partial(w_1x_1 + w_2x_2+\cdots+w_nx_n )}{\partial(w_j)} = x_j$$

综上所述：

$$\dfrac{\partial}{w_j}J(w) = (y-σ_w(x))x_j$$

所以：

$$w_j := w_j + α\sum_{i=1}^n({y^{(i)}}-σ_w(x^{(i)}))x_j^{(i)})$$

其中，$y^{(i)}$是第$i$个样本的类型；$x^{(i)}$第$i$个样本；$x_j^{(i)}$是第$i$个样本的第$j$个特征值。

5-2-2 训练算法：使用梯度上升找到最佳参数

如下图所示，有100个样本点，每一个点包含两个数值型特征：X1和X2。基于这些数据，使用梯度上升法找到最佳回归系数，也就是拟合出Logistic回归模型的最佳参数。

梯度上升法的伪代码如下：

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

创建名为logRegre.py的文件，并写入下列程序。

"""
    Function:
        读取数据
    Parameters:
        None
    Return：
        dataMat——数据集
        labelMat——数据集中每组数据对应的类别标签集
    Modify：
        2018-01-29    
"""
def loadDataSet():
    dataMat = []; labelMat = []
    #创建列表dataMat与labelMat分别存储数据集和数据集对应的类别标签
    fr = open(r'E:\Python_Files\CodeofMe\Chapter5\testSet.txt')
    #打开数据集文件
    for line in fr.readlines():
    #读取文件中的每一行
        lineArr = line.strip().split()
        #将每一行的每一个数据作为元素放在列表lineArr中
        # S.strip（[chars]）—— 返回字符串S的副本，删除前导和尾随空白
        #S.split（[sep [，maxsplit]]） - >字符串列表
        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
        #得到数据dataMat列表，每个元素是一个列表，每个列表中X0初始化为1(为了计算方便)，X1为第一个数据，X2为第二个数据
        labelMat.append(int(lineArr[2]))
        #得到类别标签的向量列表
    return dataMat,labelMat

"""
    Function:
        Sigmoid函数
    Parameters:
        inX——样本数据
    Return：
        1.0/(1+exp(-inX))——Sigmoid函数值
    Modify：
        2018-01-29    
"""
def sigmoid(inX):
    return 1.0/(1+exp(-inX))
    #sigmoid函数公式，返回函数值

"""
    Function:
        梯度上升算法
    Parameters:
        dataMatIn——准备好的数据集，二维Numpy数组
        classLabels——数据集对应的类别标签集
    Return：
        weights——进行词条切分后的文档集合
    Modify：
        2018-01-29    
"""
def gradAscent(dataMatIn, classLabels):
    dataMatrix = mat(dataMatIn)  
    #将列表dataMatIn转换成100x3的dataMatrix矩阵形式
    labelMat = mat(classLabels).transpose() 
    #将列表classLabels转换成1x100的矩阵，然后再将矩阵转置得到100x1的矩阵labelMat
    m,n = shape(dataMatrix)
    #得到矩阵dataMatrix的行数和列数
    alpha = 0.001
    #设置向目标移动的步长
    maxCycles = 500
    #设置梯度上升算法的迭代次数
    weights = ones((n,1))
    #初始化向量中每个特征Xi的最优回归系数为1
    for k in range(maxCycles):              
    #矩阵运算次数k
        h = sigmoid(dataMatrix*weights)     
        #矩阵计算，计算sigmoid的函数值，h是一个100x1的矩阵（一个样本一个值，100个样本）
        error = (labelMat - h)  
        #矢量作差（这儿是公式的推导出来的，用于回归最佳系数计算的迭代）
        weights = weights + alpha * dataMatrix.transpose()* error 
        #矩阵迭代，根据公式计算最佳回归系数
    return weights

需要注意的是，程序中对dataMat中的每个元素里的X0初始化为1，是为了方便计算。

写入上述程序前加入下列程序，以实现上述程序中的矩阵的运算：

from numpy import *

最后，在Python提示符下，敲入下面的代码：

>>> sys.path.append('E:\Python_Files\CodeofMe\Chapter5')
>>> import logRegres
>>> dataArr,labelMat = logRegres.loadDataSet()
>>> logRegres.gradAscent(dataArr,labelMat)
matrix([[ 4.12414349],
        [ 0.48007329],
        [-0.6168482 ]])

如上程序输出结果所示，得到的matrix是对应的回归系数。

5-2-3 分析数据：画出决策边界

得到一组回归系数后，可以确定不同类别数据之间的分割线，画出该分割线，使得优化的过程便于理解。

在logRegres.py添加下列代码，画出数据集和Logistic回归最佳拟合直线的函数。

"""
    Function:
        画出数据集和Logistic回归最佳拟合直线
    Parameters:
        weights—回归系数矩阵集
    Return：
        None
    Modify：
        2018-01-29    
"""
def plotBestFit(weights):
    import matplotlib.pyplot as plt
    #导入相应画图的库
    dataMat,labelMat=loadDataSet()
    #读取数据和数据对应的类别标签
    dataArr = array(dataMat)
    #将dataMat列表数据转换成矩阵
    n = shape(dataArr)[0] 
    #得到矩阵dataArr的行数
    xcord1 = []; ycord1 = []
    #初始化列表xcord1和ycord1，将属于类型1的两个数据分别作为可视化坐标上的x和y存在xcord1和ycord1
    xcord2 = []; ycord2 = []
    #初始化列表xcord2和ycord2，将属于类型0的两个数据分别作为可视化坐标上的x和y存在xcord2和ycord2    
    for i in range(n):
    #读取类型标签集的每一个标签
        if int(labelMat[i])== 1:
            xcord1.append(dataArr[i,1]); ycord1.append(dataArr[i,2])
            #标签为1，将每个样本对应的两个数据存放在xcord1和ycord1里
        else:
            xcord2.append(dataArr[i,1]); ycord2.append(dataArr[i,2])
            #标签为0，将每个样本对应的两个数据存放在xcord2和ycord2里
    fig = plt.figure()
    #建立可视化图形
    ax = fig.add_subplot(111)
    #建立一个窗口
    ax.scatter(xcord1, ycord1, s=30, c='red', marker='s')
    #将标签为1的每个样本画在可视化图像上
    ax.scatter(xcord2, ycord2, s=30, c='green')
    #将标签为0的每个样本画在可视化图形上
    x = arange(-3.0, 3.0, 0.1)
    #以0.1为步长，从-3.0到3.0显示x的刻度
    y = (-weights[0]-weights[1]*x)/weights[2]
    #该处设置了sigmoid函数为0，0是两个分类（类别1和类别0的分界处）
    #因此设定0=w0x0+w1x1+w2x2,解出如上式X2和X1的关系式，此处的X0前面的程序已经初始化为1。    
    ax.plot(x, y)
    #将决策边界线画在可视化图像上
    plt.xlabel('X1'); plt.ylabel('X2');
    #标注横轴为X1,纵轴为X2
    plt.show()
    #可视化图形显示

这段程序需要注意的是Logistic回归最佳拟合直线对应的函数。函数中设置Sigmoid函数为0。0是两个分类（类别1和类别0）的分界处。所以设定$z = 0 = w_0x_0+w_1x_1+w_2x_2$，然后解出$x_2$和$x_1$的关系式，其中$x_0 = 1$（为了方便计算，程序中初始化赋值为1）。

在Python提示符下输入：

>>> from numpy import *
>>> reload(logRegres）     
>>> reload(logRegres)
<module 'logRegres' from 'E:\Python_Files\CodeofMe\Chapter5\logRegres.py'>
>>> logRegres.plotBestFit(weights.getA())

getA()是将输入转换成Numpy数组。

输出结果如下：

5-2-4 训练算法：随机梯度上升

梯度上升算法存在的缺点：每次更新回归系数时都需要遍历整个数据集。

假设有数十亿样本和成千上万的特征，那该方法的计算复杂程度相当高。

改进方法是采用随机梯度上升算法。

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

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

随机梯度上升算法的伪代码：

所有回归系数初始化为1
对数据集中每个样本   
      计算该样本的梯度
      使用alpha x gradient更新回归系数值
返回回归系数值

随机梯度上升算法的代码如下：

"""
    Function:
        随机梯度上升算法计算回归系数
    Parameters:
        dataMatrix——数据集（此处应该输入数据为Numpy数组）
        classLabels——数据集对应的标签列表
    Return：
        weights——回归系数
    Modify：
        2018-01-30    
"""
def stocGradAscent0(dataMatrix, classLabels):
    m,n = shape(dataMatrix)
    #获取dataMatrix数据数组的行列数
    alpha = 0.01
    #设置向目标移动的步长
    weights = ones(n)   
    #将所有回归系数初始化为1
    for i in range(m):
    #对每一行数据进行循环，一行数据一个样本
        h = sigmoid(sum(dataMatrix[i]*weights))
        #仅用计算一个样本的Sigmoid函数值
        error = classLabels[i] - h
        #推导出的计算回归系数所需的计算因式（一个数值，非矩阵）
        weights = weights + alpha * error * dataMatrix[i]
        #根据公式计算回归系数
    return weights

随机梯度上升算法与梯度上升算法代码上很相似，但也有区别：

1. 随机梯度上升算法全是数值，而梯度上升算法的变量h和误差error都是向量。
2. 随机梯度上升算法没有矩阵的转换过程，所以变量的数据类型都是Numpy数组。

在Python提示符下输入下列程序，验证该方法的结果：

>>> from numpy import *
>>> import sys
>>> sys.path.append('E:\Python_Files\CodeofMe\Chapter5')
>>> import logRegres
>>> reload(logRegres)
<module 'logRegres' from 'E:\Python_Files\CodeofMe\Chapter5\logRegres.pyc'>
>>> dataArr,labelMat = logRegres.loadDataSet()
>>> weights = logRegres.stocGradAscent0(array(dataArr),labelMat)
>>> weights
array([ 1.01702007,  0.85914348, -0.36579921])
>>> logRegres.plotBestFit(weights)

输出结果如图所示：

不难看出，按照上面的程序，得到的回归系数只经过一次迭代得到的，相比于梯度上升算法的500次遍及所有数据的迭代，只有一次迭代的随机梯度上升算法分类器错分的结果比梯度上升算法多很多。

我们将stocGradAscent0(dataMatrix, classLabels)函数稍作修改，看一下增加迭代次数与回归系数的关系如何，以及对应分类效果如何。

修改后的程序如下：

"""
    Function:
        随机梯度上升算法计算回归系数
    Parameters:
        dataMatrix——数据集（此处应该输入数据为Numpy数组）
        classLabels——数据集对应的标签列表
        IterationNum——计算回归系数迭代的次数
    Return：
        weights——回归系数
    Modify：
        2018-01-30    
"""
def stocGradAscent0(dataMatrix, classLabels, IterationNum):
    import matplotlib.pyplot as plt
    m,n = shape(dataMatrix)
    #获取dataMatrix数据数组的行列数
    alpha = 0.01
    #设置向目标移动的步长
    weights = ones(n)   
    #将所有回归系数初始化为1
    x0cord1 = [];x1cord1 = [];x2cord1 = []
    for j in range(IterationNum):
        for i in range(m):
        #对每一行数据进行循环，一行数据一个样本
            h = sigmoid(sum(dataMatrix[i]*weights))
            #仅用计算一个样本的Sigmoid函数值
            error = classLabels[i] - h
            #推导出的计算回归系数所需的计算因式（一个数值，非矩阵）
            weights = weights + alpha * error * dataMatrix[i]
            #根据公式计算回归系数
            x0cord1.append(weights[0])
            x1cord1.append(weights[1])
            x2cord1.append(weights[2])
    fig = plt.figure()
    #建立可视化图形
    ax0 = fig.add_subplot(311)
    ax1 = fig.add_subplot(312)
    ax2 = fig.add_subplot(313)
    #建立3x1的窗口
    #ax0.scatter(range(500*m), x0cord1, s=0.0001, c='blue')
    #ax1.scatter(range(500*m), x1cord1, s=0.0001, c='blue')
    #ax2.scatter(range(500*m), x2cord1, s=0.0001, c='blue')
    #在可视化窗口画出散点图
    ax0.plot(range(IterationNum*m), x0cord1, marker='.' ,markersize = 0.000001)
    ax1.plot(range(IterationNum*m), x1cord1, marker='.' ,markersize = 0.000001)
    ax2.plot(range(IterationNum*m), x2cord1, marker='.' ,markersize = 0.000001)
    #在可视化窗口画出折线图
    plt.show()
    #图像显示
    return weights

在Python命令提示符下，我们输入下列代码，看一下迭代次数为1次，100次，200次对应的回归系数与迭代次数的关系，以及对应的分类的效果：

>>> reload(logRegres)
<module 'logRegres' from 'E:\Python_Files\CodeofMe\Chapter5\logRegres.py'>
>>> dataMat,labelMat = logRegres.loadDataSet()
>>> weights = logRegres.stocGradAscent0(array(dataMat),labelMat,1)
>>> logRegres.plotBestFit(weights)
>>> weights = logRegres.stocGradAscent0(array(dataMat),labelMat,100)
>>> logRegres.plotBestFit(weights)
>>> weights = logRegres.stocGradAscent0(array(dataMat),labelMat,200)
>>> logRegres.plotBestFit(weights)

一次迭代：

100次迭代：

200次迭代：

回归系数与迭代次数的关系图中，横坐标是迭代的次数（因为样本是100，所以刻度点是迭代次数的100倍），纵坐标是回归系数的值。

从图像的变化情况来看，迭代次数到达200次，回归系数差不到收敛于某个数。所以对应的分类效果图也比迭代一次（回归系数还未达到收敛）的好很多。

通过迭代次数与回归系数的图像我们看到，即使最后回归系数有收敛的趋势，但是在收敛范围内，依然有小的周期性波动。产生这种现象的原因是存在一些不能正确分类的样本点（数据集并非线性可分），在每次迭代时会引发系数的剧烈改变。

为了让算法避免来回波动，从而收敛到某个值，同时提高收敛速度，对随机梯度上升算法进行修改：

"""
    Function:
        改进随机梯度上升算法计算回归系数
    Parameters:
        dataMatrix——数据集（此处应该输入数据为Numpy数组）
        classLabels——数据集对应的标签列表
        numIter——计算回归系数迭代的次数
    Return：
        weights——回归系数
    Modify：
        2018-01-30    
"""
def stocGradAscent1(dataMatrix, classLabels, numIter):
    m,n = shape(dataMatrix)
    #获取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    
            #apha 每次迭代会进行调整 缓解回归系数的波动或者高频波动
            randIndex = int(random.uniform(0,len(dataIndex)))
            #通过随机选取样本来更新回归系数以减少周期性的波动
            h = sigmoid(sum(dataMatrix[randIndex]*weights))
            #仅用计算一个样本的Sigmoid函数值
            error = classLabels[randIndex] - h
            #推导出的计算回归系数所需的计算因式（一个数值，非矩阵）
            weights = weights + alpha * error * dataMatrix[randIndex]
            #根据公式计算回归系数
            del(dataIndex[randIndex])
            #删掉已迭代过的样本
    return weights

算法改进有三处：

1.alpha在每次迭代的时候都会调整，这会缓解原来算法上回归系数的波动或者高频波动。虽然alpha会随着迭代次数不断减少，但永远不会减小到0，因为

alpha = 4/(1.0+j+i)+0.0001

中还存在一个常数项。必须这样做的原因是为了保证在多次迭代之后新数据仍然具有一定的影响。如果要处理的问题是动态变化的，可以适当加大上述常数项，来确保新的值获得更大的回归系数。
2. 通过随机选取样本来更新回归系数。这种方法可以减少周期性的波动，这种方法每次随机从列表中选出一个值，然后从列表中删除该值（再进行下一次迭代）。
3. 增加了一个迭代次数的输入参数。

在Python的命令提示符下输入下列代码，与没改进的随机上升梯度算法输出结果进行对比。

>>> reload(logRegres)
<module 'logRegres' from 'E:\Python_Files\CodeofMe\Chapter5\logRegres.pyc'>
>>> dataArr,labelMat = logRegres.loadDataSet()
>>> weights = logRegres.stocGradAscent1(array(dataArr),labelMat,50)
>>> logRegres.plotBestFit(weights)

输出结果如下：

我们再看看没有改进的算法50次迭代的效果：

5-3 示例：从疝气病症预测病马的死亡率

使用Logistic回归来预测患有疝病的马的存货问题。本例数据包含368个样本和28个特征。

示例：使用Logistic回归估计马疝病的死亡率

1. 收集数据：给定数据文件。
2. 准备数据：用Python解析文本文件并填充缺失值
3. 分析数据：可视化并观察数据。
4. 训练算法：使用优化算法，找到最佳的系数。
5. 测试算法：为了量化回归的效果，需要观察错误率。根据错误率决定是否会退到训练阶段，通过改变迭代的次数和步长等参数来得到更好的回归系数。
6. 使用算法：实现一个简单的命令行程序来收集马的症状并输出预测结果并非难事。

5-3-1 准备数据：处理数据中的缺失值

首先来理解一下数据缺失这个问题：假设有100个样本和20个特征，这些数据都是机器收集回来的。若机器上的某个传感器损坏导致一个特征无效时该怎么办？此时是否要扔掉整个数据？这种情况下，另外19个特征怎么办？他们是否还可用？答案是肯定的。因为有时候数据相当昂贵，扔掉和重新获取都是不可取的，所以可采取下列方法应对数据缺失的问题：

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

对现数据集进行预处理，要做两件事：

第一，所有的缺失值必须用一个实数值必须用一个实数值来替换，因为我们使用的Numpy数据类型不允许包含缺失值。这里选择实数0来替代所有缺失值，因为sigmoid(0)=0.5时，他对结果的预测不具有任何倾向性，他是一个在更新时不会影响系数的值的选择。满足“特殊值的要求”。

第二，如果在测试集中发现了亿条数据的类别标签已经缺失，那么简单的做法是将该条数据丢弃。

5-3-2 测试算法：用Logistic回归进行分类

使用Logistic回归方法进行分类并不需要做很多工作，所需要做的只是：

将测试集上每个特征向量乘以最优化方法得来的回归系数，再将该乘积求和，最后输入到Sigmoid函数中，得到的Sigmoid函数值大于0.5就预测类标签为1，反之为0。

在logRegres.py文件里添加下列函数：

"""
    Function:
        分类函数
    Parameters:
        inX——测试集样本特征向量
        weights——训练集训练出的最优回归系数
    Return：
        0/1——测试集类别标签
    Modify：
        2018-01-30    
"""
def classifyVector(inX, weights):
    prob = sigmoid(sum(inX*weights))
    #计算测试集样本的Sigmoid函数值
    if prob > 0.5: return 1.0
    #函数值大于0.5，返回标签1，否则返回0
    else: return 0.0

"""
    Function:
        病马预测函数
    Parameters:
        None
    Return：
        errorRate——预测的错误率
    Modify：
        2018-01-30    
"""
def colicTest():
    frTrain = open(r'E:\Python_Files\CodeofMe\Chapter5\horseColicTraining.txt')
    #读取训练集数据
    frTest = open(r'E:\Python_Files\CodeofMe\Chapter5\horseColicTest.txt')
    #读取测试集数据
    trainingSet = []; trainingLabels = []
    #创建训练集数据列表和数据对应的类别标签列表
    for line in frTrain.readlines():
    #读取训练集的每一行数据
        currLine = line.strip().split('\t')
        #将数据集的每行数据每个字符串类型的特征+类别转换成列表存放（去除了空格）
        lineArr =[]
        #创建存放每行数据的列表
        for i in range(21):
        #疝病有20特征数据+1个类别标签数据
            lineArr.append(float(currLine[i]))
            #将每行数据转换成float类型存放
        trainingSet.append(lineArr)
        #以float类型存放每个样本的特征和类别标签
        trainingLabels.append(float(currLine[21]))
        #得到float类型的类别标签集
    trainWeights = stocGradAscent1(array(trainingSet), trainingLabels, 1000)
    #调用优化的随机梯度上升算法计算最优回归系数，迭代次数为1000次
    errorCount = 0; numTestVec = 0.0
    #初始化测试错误计数参数和
    for line in frTest.readlines():
    #读取测试集数据的每一行
        numTestVec += 1.0
        #测试个数自动加1
        currLine = line.strip().split('\t')
        #将测试集每一行数据转换成字符串列表存储
        lineArr =[]
        #建立存放每行数据的列表
        for i in range(21):
        #疝病有20特征数据+1个类别标签数据
            lineArr.append(float(currLine[i]))
            #将每行数据转换成float类型存放
        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

"""
    Function:
        调用colicTest()10次求结果的平均值进行评估
    Parameters:
        None
    Return：
        None
    Modify：
        2018-01-30    
"""
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))

Python命令提示符下输入下列代码进行测试：

>>> reload(logRegres)
<module 'logRegres' from 'E:\Python_Files\CodeofMe\Chapter5\logRegres.py'>
>>> logRegres.multiTest()
the error rate of this test is: 0.313433
the error rate of this test is: 0.402985
the error rate of this test is: 0.417910
the error rate of this test is: 0.313433
the error rate of this test is: 0.358209
the error rate of this test is: 0.283582
the error rate of this test is: 0.313433
the error rate of this test is: 0.343284
the error rate of this test is: 0.313433
the error rate of this test is: 0.432836
after 10 iterations the average error rate is: 0.349254

从上面的结果来看，10次迭代之后的平均错误率为34.9254%。

5-4 本章小结

Logistic回归目的：寻找一个非线性函数Sigmoid的最佳拟合参数。

求解最佳拟合参数常用梯度上升算法，该算法可以简化为随机梯度上升算法。

随机梯度上升算法与梯度上升算法效果相当，但前者占用的计算资源更少。同时，前者是一个在线算法可以在新数据到来时就完成参数更新，不需要重新读取整个数据集来进行批处理运算。

5-5 参考文献

《机器学习实战》