【Machine Learning in Action --5】逻辑回归（LogisticRegression）

1、概述

　　Logistic regression（逻辑回归）是当前业界比较常用的机器学习方法，用于估计某种事物的可能性。

　　在经典之作《数学之美》中也看到了它用于广告预测，也就是根据某广告被用 户点击的可能性，把最可能被用户点击的广告摆在用户能看到的地方，然后叫他“你点我啊！”用户点了，你就有钱收了。这就是为什么我们的电脑现在广告泛滥的 原因。还有类似的某用户购买某商品的可能性，某病人患有某种疾病的可能性啊等等。这个世界是随机的（当然了，人为的确定性系统除外，但也有可能有噪声或产生错误的结果，只是这个错误发生的可能性太小了，小到千万年不遇，小到忽略不计而已），所以万物的发生都可以用可能性或者几率（Odds）来表达。“几率”指的是某事物发生的可能性与不发生的可能性的比值。

       Logistic regression可以用来回归，也可以用来分类，主要是二分类。

2、基本理论

2.1Logistic regression和Sigmoid函数

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

　　Sigmoid函数具体的计算公式如下：

　　　　

　　        z=w₀x₀+w₁x₁+w₂x₂+...+w_nx_n， z=w^Tx 　其中w是我们要找的最佳参数（系数），x是分类器的输入数据特征。

　　当x为0时，Sigmoid函数值为0.5，随着x的增大，对应的Sigmoid值将逼近于1；而随着x的减小，Sigmoid值将逼近于0。如果横坐标刻度足够大（如下图所示），Sigmoid函数看起来很像一个阶跃函数。

    

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

2.2最优化理论

 　　由上述问题得到，我们现在的问题变成了：最佳回归系数时多少？

　　  z=w₀x₀+w₁x₁+w₂x₂+...+w_nx_n， z=w^Tx 

　　向量x是分类器的输入数据，向量w是我们要找的最佳参数（系数），从而使得分类器尽可能地精确，为了寻找最佳参数，需要用到最优化理论的一些知识。

　　下面首先介绍梯度上升的最优化方法，我们将学习到如何使用该方法求得数据集的最佳参数。接下来，展示如何绘制梯度上升法产生的决策边界图，该图能将梯度上升法的分类效果可视化地呈现出来。最后我们将学习随机梯度上升法，以及如何对其进行修改以获得更好的结果。

2.2.1梯度上升法

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

    

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

　　　　　　

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

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

　　

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

　　事例（用梯度上升找到最佳参数）

　　该数据中有100个样本点，每个点包含两个数值型特征：X1和X2。在此数据集上，我们将通过使用梯度上升法找到最佳回归系数，也就是说拟合出Logistic回归模型的最佳参数。在以下的数据集中，每行的前两个值分别是X1和X2，它们是数据特征，第三个值是数据对应的类别标签，为了方便计算，该函数还将X0的值设为1.0。

　　testSet.txt数据集如下：

　　

　　梯度上升算法代码如下，建立一个logRegres.py的文件，添加如下代码：

#!/usr/bin/python
# -*- coding: utf-8 -*-
from numpy import *
#Logistic回归梯度上升优化算法
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，X1，X2的值
        labelMat.append(int(lineArr[2]))    #获取最后一列的类别标签
    return dataMat,labelMat

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

def gradAscent(dataMatIn, classLabels):     #dataMatIn是一个2维Numpy数组
    dataMatrix = mat(dataMatIn)             #convert to NumPy matrix
    labelMat = mat(classLabels).transpose() #convert to NumPy matrix
    m,n = shape(dataMatrix)
    alpha = 0.001
    maxCycles = 500                         #迭代次数为500
    weights = ones((n,1))                   #回归系数初始化为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

 　　在python提示符下，写下面代码：

>>>import logRegres
>>> dataArr,labelMat=logRegres.loadDataSet()
>>> logRegres.gradAscent(dataArr,labelMat)
matrix([[ 4.12414349],
        [ 0.48007329],
        [-0.6168482 ]])

解释：

>>> dataMat
matrix([[  1.00000000e+00,  -1.76120000e-02,   1.40530640e+01],
        [  1.00000000e+00,  -1.39563400e+00,   4.66254100e+00],
           ...
        [  1.00000000e+00,   1.38861000e+00,   9.34199700e+00],
        [  1.00000000e+00,   3.17029000e-01,   1.47390250e+01]])
>>> labelMat
matrix([[0],
        [1],
        ...
        [0],
        [0]])

　　上面解出了一组回归系数，它确定了不同类别数据之间的分隔线。那么怎样画出该分隔线，从而使得优化过程便于理解，打开logRegres.py添加如下代码：

#画出数据即和Logisitic回归最佳拟合直线的函数
def plotBestFit(weights):
    import matplotlib.pyplot as plt
    dataMat,labelMat=loadDataSet()
    dataArr = array(dataMat)
    n = shape(dataArr)[0] 
    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]-weights[1]*x)/weights[2]  #W0X0+W1X1+W2X2=0(X0=1),X2=(-W0-W1X1)/W2
    ax.plot(x, y)
    plt.xlabel('X1'); plt.ylabel('X2');
    plt.show()

 　　运行程序清单，在python提示符下输入：

>>> from numpy import *
>>> reload(logRegres)
<module 'logRegres' from 'logRegres.pyc'>
>>> logRegres.plotBestFit(wei)
'''

　　这个分类结果相当不错，从图上看只错分了两到四个点。但是，尽管例子简单且数据集很小，但是这个方法确需要大量的计算（300次乘法），下一节将对该算法稍作改进，从而使它能用在真实数据集上。

2.2.2随机梯度上升

　　梯度上升算法在每次更新回归系数时都需要遍历整个数据集，该方法在处理1001个左右的数据集时尚克，但如果有数十亿样本和成千上万的特征，那么该方法的计算复杂度就太高了。

　　一种改进方法是一次仅用一个样本点来更新回归系数，该方法称为随机梯度上升算法。由于可以在信仰本到来时对分类器进行增量式更新，因而随机梯度上升算法是一个在线学习算法。

　　随机梯度上升算法代码如下，将下面代码添加到logRegres.py中：

#随机梯度上升函数
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

 　　为了验证该方法的结果，我们在python提示符中输入以下命令：

>>> from numpy import *
>>> reload(logRegres)
<module 'logRegres' from 'logRegres.py'>
>>> dataArr,labelMat=logRegres.loadDataSet()
>>> weights=logRegres.stocGradAscent0(array(dataArr),labelMat)
>>> logRegres.plotBestFit(weights)

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

　　梯度上升算法：变量h和误差error都是向量，有矩阵的转换过程

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

　　由上述图得到的最佳拟合直线图，但拟合的直线没有那么完美，这里的分类器错分了三分之一的样本。一个判断优化算法优劣的可靠方法是看它是否收敛，也就是说参数是否达到了稳定值，是否还会不断地变化。

　　对上述随机梯度上升算法上做了些修改，使其在整个数据集上运行200次，最终绘制的三个回归系数的变化情况如图所示：

　　上图展示了随机梯度上升算法在200次迭代过程中（我们的数据库有100个样本，每个样本都对系数调整一次，所以共有200*100=20000次调整）回归系数的变化情况。由图中看到X2只经过了50次迭代就达到了稳定值，但X0和X1则需要更多次的迭代。另外值得注意的是，在大的波动停止后，还有一些小的波动。产生这种现象的原因是存在一些不能正确分类的样本点（数据集并非线性可分），在每次迭代时会引发系数的剧烈改变。我们期望算法能避免来回波动，从而收敛到某个值。另外，收敛速度也要加快。

2.2.3改进的随机梯度上升算法

　　对上述随机梯度上升算法，我们做两处改进来避免上述的波动问题：

1） 在每次迭代时，调整更新步长alpha的值。随着迭代的进行，alpha越来越小，这会缓解系数的高频波动（也就是每次迭代系数改变得太大，跳的跨度太 大）。当然了，为了避免alpha随着迭代不断减小到接近于0（这时候，系数几乎没有调整，那么迭代也没有意义了），我们约束alpha一定大于一个稍微大点的常数项，具体见代码。

2）每次迭代，改变样本的优化顺序。也就是随机选择样本来更新回归系数。这样做可以减少周期性的波动，因为样本顺序的改变，使得每次迭代不再形成周期性。

　　将下述代码添加到logRegres.py中：

#改进的随机梯度上升函数
def stocGradAscent1(dataMatrix, classLabels, numIter=150):
    m,n = shape(dataMatrix)
    weights = ones(n)   #initialize to all ones
    for j in range(numIter):
        dataIndex = range(m)
        for i in range(m):
            alpha = 4/(1.0+j+i)+0.01    #j是迭代次数，i是样本点的下标，alpha每次迭代时需要调整 
            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

 　　为了验证该方法的结果，我们在python提示符中输入以下命令：

>>> reload(logRegres)
<module 'logRegres' from 'logRegres.py'>
>>> dataArr,labelMat=logRegres.loadDataSet()
>>> weights=logRegres.stocGradAscent1(array(dataArr),labelMat)
>>> logRegres.plotBestFit(weights)

　　

　　改进算法增加了一个迭代次数作为第三个参数，如果该参数没有给定的话，算法将默认迭代150次，如果给定，那么算法将按照新的参数值进行迭代。与随机梯度上升算法中的回归系数类似，改进的随机梯度上升算法中的各个回归系数的变化情况如下：

　　比较随机梯度上升和改进后的梯度上升，可以看到两点不同：

　　1）系数不再出现周期性波动。

　　2）系数可以很快的稳定下来，也就是快速收敛。这里只迭代了20次就收敛了。而上面的随机梯度下降需要迭代200次才能稳定。

3、Logistic回归的一般过程

　　（1）收集数据：采用任何方法收集数据

　　（2）准备数据：数据类型是数值型的，另外，结构化数据格式则最佳

　　（3）分析数据：采用任意方法对数据进行分析

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

　　（5）测试算法：一旦训练步骤完成，分类将会很快

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

4、Logistic回归优缺点

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

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

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