Logistic回归算法及实践

参考书籍《机器学习实战》

我们的生活中遇到很多最优化问题，比如如何在最短时间从A到达B？如果投入最少的工作量获得最大效益？等等...现在，我们假设有一些数据点，我们用一条直线对这些点进行拟合，这个拟合的过程就称为回归。利用Logistic回归进行分类的主要思想是：根据现有数据对分类边界线建立回归公式，依次进行分类。

Logistic回归的一般过程

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

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

Logistic回归

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

我们想要的函数应该是，能接受所有的输入然后预测出类别。例如在两个类的情况下输出0或者1。 Heaviside step funtion就有这种功能，但它的问题在于该函数在跳跃点上直接由0跳跃到1，这个瞬间的跳跃过程很难处理。幸好，另一个函数也有类似的性质，也更容易处理，它就是Sigmoid函数。

图5-1给出了Sigmoid函数在不同坐标尺度下的两条曲线。当x逼近于0时，Sigmoid函数值为0.5，随着x增大，y逼近1；x减小，y逼近0.如果横坐标刻度足够大，Sigmoid函数很像阶跃函数。

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

确定了分类器的函数形式后，问题就变成了：最佳回归系数是多少？如何确定？

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

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

采用向量的写法，上面公式可以写成z=wx，向量w就是我们要找的最佳参数。接下来我们将学习到如何使用方法求得数据集的最佳参数。

梯度上升法

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

参考图5-2

图5-2中梯度上升算法沿梯度方向移动了一步。可以看到，梯度算子总是指向函数值增长最快的方向。该量值称为步长，记做α。用向量来表示的话，梯度上升算法的迭代公式如下：

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

---

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

图5-3中有100个样本点，每个点包含两个数值型特征：X1和X2。

梯度上升的伪代码如下：

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

创建文件logRegres.py，输入以下代码：

#coding=utf-8
from numpy import *
def loadDataSet():
    dataMat = []
    labelMat = []
    fr = open('testSet.txt')
    for line in fr.readlines():
        lineArr = line.strip().split()
        #回归系数，坐标1，坐标2
        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 = mat(classLabels).transpose()#训练集标签
    m,n = shape(dataMatrix)#m为训练样本数，n为每个样本的属性数
    alpha = 0.001
    maxCycles = 500
    weights = ones((n,1))
    for k in range(maxCycles):
        h = sigmoid(dataMatrix*weights)
        error = (labelMat - h)
        weights = weights + alpha * dataMatrix.transpose() * error
    return weights

上述函数主要功能为：
loadDataSet：从TreeSet.txt读取样本
sigmoid：阶跃函数
gradAscent：将得到的样本通过梯度上升算法，反复迭代后得到对样本每个属性的权重（回归系数）

为了验证实际效果，输入代码：

import logRegres
dataArr,labelMat = logRegres.loadDataSet()
print(logRegres.gradAscent(dataArr,labelMat))

分析数据：画出决策边界

上面接解出了回归系数，他确定了不同类别数据之间的分隔线。我们应该如何画出分割线。从而使得优化的过程便于理解呢？
在logRegres.py添加如下代码：

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]
    ax.plot(x,y)
    plt.xlabel('X1')
    plt.ylabel('X2')
    plt.show()

在这个案例中，0是两个分类的分解处，因此我们设定0=w0x0+w1x1+w2x2，然后解出X2与X1的关系式。
输入测试代码：

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

可以看到分类的结果相当不错。

训练算法：随机梯度上升

虽然分类结果不错，但是可以发现，在数据集很小的情况下，需要大量计算（300次乘法）。如果数据上升到成千上万的特征。该方法计算复杂度就太高了。一种改进方法是一次仅用一个样本点来更新回归系数，该方法称为随机梯度上升算法。
随机梯度算法可以写成如下伪代码：
所有回归系数初始化为1
对数据集中每个样本：
计算该样本的梯度
使用alpha x gradient更新回归系数值
返回回归系数值
在logRegres.py中添加如下代码：

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

测试代码：

from numpy import *
import logRegres
dataArr,labelMat = logRegres.loadDataSet()
#weights = logRegres.gradAscent(dataArr,labelMat)
#logRegres.plotBestFit(weights.getA())
weights = logRegres.stocGradAscent0(array(dataArr),labelMat)
logRegres.plotBestFit(weights)

可以发现，分类器分错了三分之一的样本。当然直接比较两个分类结果是不公平的，前者的结果是在整个数据集上迭代500次得到。所以我们对随机梯度算法做一些改进，使其在数据集上运行多次，最终会址的三个回归系数变化情况如图所示

系数2，也就是X2只经过了50次迭代就达到了稳定值，但系数1和0则需要更多次迭代。另外，在大的波动结束后，依然有一些小的周期性波动。这是因为存在一些不能正确分类的样本点，每次迭代时会引发系数的剧烈改变。我们期望算法能回避来回波动，从而收敛到某个值。另外，收敛的速度也要加快。
在logRegres.py中添加如下代码：

def stocGradAscent1(dataMatrix,classLabels,numIter=150):
    m,n = shape(dataMatrix)
    weights = ones(n)
    for j in range(numIter):
        dataIndex = 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

代码改进两处，alpha在每次迭代的时候都会调整，这会缓解图5-6的数据流动或高频波动。另外alpha会随着迭代次数不断变小，但永不会减小到0，这样保证在多次迭代后新数据依然会有影响。如果要处理的问题是动态变化的，可以适当加大常数项。另外一点要注意的是，在降低alpha的函数中，alpha每次减少1/(j+i)，其中j是迭代次数，i是样本点的下标，这样当j<<max(i)时，alpha就不是严格下降的。第二个改进地方在于随机选取更新，它将减少周期性的波动。

输入测试代码：

from numpy import *
import logRegres
dataArr,labelMat = logRegres.loadDataSet()
#weights = logRegres.gradAscent(dataArr,labelMat)
#logRegres.plotBestFit(weights.getA())
weights = logRegres.stocGradAscent1(array(dataArr),labelMat)
logRegres.plotBestFit(weights)

效果已经很接近之前的算法，同时减少了计算复杂度和迭代次数。