机器学习：Logistic回归

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

1.1 Logistic回归的一般过程

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

1.2 Sigmoid函数

我们想要的函数应该是，能接受所有的输入然后预测出类别。例如，在两个类的情况下，上述函数输出0或1，具有这种性质的函数称为海维赛德阶跃函数，或者直接称为单位阶跃函数。然而这个函数的问题在于：该函数在跳跃点上从0瞬间跳跃到1，这个瞬间跳跃有时很难处理。而另一个函数也有类似性质，且数学上更易处理，这就是Sigmoid函数。该函数具体的计算公式如下：

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

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

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

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

2.1 梯度上升算法

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


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

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

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

2.2 梯度下降算法

该算法与上升算法一样，只是公式中的加法需要变成减法。公式为：

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

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

2.3 使用梯度上升找到最佳参数

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

from numpy import *
def loadDataSet():
    dataMat = []; labelMat = []
    fr = open('testSet.txt') #打开文本
    # 逐行读取，每行前两个值x1和x2，第三个值是数据对应的类别标签，将x0的值设置为1.0
    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

# 定义Sigmoid函数
def sigmoid(inX):
    return 1.0/(1+exp(-inX))

# 梯度上升算法
# dataMathIn参数是一个二维Numpy数组，每列分别代表每个不同的特征，每行则代表每个训练样本
# dataMathIn存放100x3的矩阵
# classLabels参数是类别标签，是一个1x100的行向量
def gradAscent(dataMatIn, classLabels): 
    dataMatrix = mat(dataMatIn)  # 原始数据转化为矩阵          
    # 行向量转换为列向量(转置)赋值给labelMat
    labelMat = mat(classLabels).transpose() 
    # 得到矩阵大小
    m,n = shape(dataMatrix)
    alpha = 0.001 # 向目标移动的步长
    maxCycles = 500 #迭代次数500
    weights = ones((n,1)) #初始化参数为1
    # 矩阵运算
    # 计算真实类别与预测类别的差值，然后按照该差值方向调整回归系数
    for k in range(maxCycles):              
        h = sigmoid(dataMatrix*weights)     # h为列向量，元素个数等于样本个数
        error = (labelMat - h)              # 矢量减法
        weights = weights + alpha * dataMatrix.transpose()* error 
    return weights

2.4 分析数据画出决策边界

画出数据集和Logistic回归最佳拟合直线的函数，代码实现如下：

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] # 此处设置了Sigmoid函数为0
    ax.plot(x, y)
    plt.xlabel('X1'); plt.ylabel('X2');
    plt.show()

之前提到，0是两个类别的分界处。因此，我们设定0 = w0x0 + w1x1 + w2x2，然后解出X2和X1的关系式（即分离线的方程，注意X0 = 1）
尽管数据集小，整个方法却需要大量计算(300次乘法)，接下来我们进行改进。

三、训练算法：随机梯度上升算法

3.1 随机梯度上升算法

梯度上升算法在每次更新回归系数时都需要遍历整个数据集，导致该方法在处理数十亿样本和成千上万特征时计算的复杂度就太高了。一种改进方法是一次仅用一个样本点来更新回归系数，该方法称为随机梯度上升算法。由于可以在新样本到来时对分类器进行增量式更新，因而随机梯度上升算法是一个在线学习算法，所以一次处理所有数据被称作批处理。

随机梯度上升算法的伪代码如下：
所有回归系数初始化为1
对数据集中的每个样本
计算该样本梯度
使用alpha x gradient更新回归系数值
返回回归系数值

随机梯度上升算法代码实现：

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

可以看到，随机梯度上升算法和梯度上升算法在代码上相似，区别是：
1.后者的变量h和误差error都是向量，而前者是数值；
2.前者没有矩阵转换过程，所有变量的数据类型都是Numpy数组。

3.2 随机梯度上升算法改进

一个判断优化算法优劣的可靠方法是看它是否收敛，也就是说参数是否达到了稳定值，是否还会不断变化？运行随机梯度上升算法，在数据集的一次遍历中回归系数与迭代次数的关系图显示，回归系数经过大量迭代才能达到稳定值，并且仍然有局部波动现象，产生这种现象的原因是存在一些不能正确分类的样本点(数据集并非线性可分)，在每次迭代时会引发系数的剧烈改变。而我们期望算法能避免来回波动，从而收敛到某个值，而且收敛速度也需要加快。优化代码如下：

def stocGradAscent1(dataMatrix, classLabels, numIter=150):
    m,n = shape(dataMatrix)
    weights = ones(n)   #initialize to all ones
    for j in range(numIter):
        dataIndex = list(range(m))
        for i in range(m):
            alpha = 4/(1.0+j+i)+0.0001    # 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

程序改进两处代码：
1.alpha = 4/(1.0+j+i)+0.0001，也就是说alpha每次迭代时需要调整，这将缓解数据波动或高频波动。虽然alpha会随着迭代次数不断减小，但永远不会减到0，因为公式中有常数项。这样做的原因是为了保证在多次迭代之后新数据仍然有一定的影响。如果要处理的问题是动态变化的，可以适当加大常数项的值来确保新的值获得更大的回归系数。

2.另一处是通过随机选取样本来更新回归系数，这将减少周期性波动。这种方法每次随机从列表中选一个值，然后从列表中删掉该值再进行下一次迭代。

此外，改进算法还增加了一个迭代次数作为第三参数，若没有给定该参数算法默认迭代150次。

四、从疝气病症预测病马的死亡率

4.1 准备数据：处理数据中的缺失值

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

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

2.如果数据集中一条数据的类标签已经缺失，简单做法是将该条数据丢弃。因为类别标签与特征不同，很难确定采用某个合适的值来替换。

原始数据集经过预处理后保存成两个文件：horseColicTest.txt和horseColic-Training.txt。

4.2 测试算法：用Logistic回归进行分类

使用Logistic回归方法进行分类所需要做的只是把测试集上每个特征向量乘以最优化方法得来的回归系数，再将该乘积结果求和，最后输入到Sigmoid函数中即可。如果对应的Sigmoid值大于0.5就预测为1，否则为0。

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('\t')
        lineArr =[]
        for i in range(21):
            lineArr.append(float(currLine[i]))
        trainingSet.append(lineArr)
        trainingLabels.append(float(currLine[21]))
    trainWeights = stocGradAscent1(array(trainingSet), trainingLabels, 1000)
    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("此测试的错误率为: %f" % errorRate)
    return errorRate

# 调用函数colicTest()函数10次求结果的平均值
def multiTest():
    numTests = 10; errorSum=0.0
    for k in range(numTests):
        errorSum += colicTest()
    print("在%d次迭代之后，平均错误率为：%f" % (numTests, errorSum/float(numTests)))

运行结果如图：

分析：从结果可以得到，10次迭代之后的平均错误率为34%。由于有30%的数据缺失，所以这个结果还是不错的，如果调整colicTest()中的迭代次数和stochGradAscent1()中的步长，平均错误率可以降到20%左右。

4.3 小结

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

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