啃书-《机器学习实战》：Logistic回归

假设现在有一些数据点，我们用一条直线对这些点进行拟合（该线称为佳拟合直线），这个拟合过程就称作回归。利用Logistic回归进行分类的主要思想是：根据现有数据对分类边界线建立回归公式，以此进行分类。这里的“回归”一词源于最佳拟合，表示要找到佳拟合参数集，其背后的数学分析将在下一部分介绍。训练分类器时的做法就是寻找最佳拟合参数，使用的是优化算法。接下来介绍这个二值型输出分类器的数学原理。

1. 简介

Logistic回归的一般过程：

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

Logistic回归优缺点：

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

2. sigmoid函数

我们想要的函数应该是，能接受所有的输入然后预测出类别。例如，在两个类的情况下，上述函数输出0或1。或许你之前接触过具有这种性质的函数，该函数称为海维塞德阶跃函数（Heaviside step function），或者直接称为单位阶跃函数。然而，海维塞德阶跃函数的问题在于：该函数在跳跃点上从0瞬间跳跃到1，这个瞬间跳跃过程有时很难处理。幸好，另一个函数也有类似的性质，且数学上更易处理，这就是Sigmoid函数。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回归也可以被看成是一种概率估计。
确定了分类器的函数形式之后，现在的问题变成了：最佳回归系数是多少？ 如何确定它们的大小？这些问题将在下一节解答。

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

Sigmoid函数的输入记为z，由下面公式得出：
z = w₀x₀ + w₁x₁ + w₂x₂ + … +w_nx_n
如果采用向量的写法，上述公式可以写成z = w^Tx，它表示将这两个数值向量对应元素相乘然后全部加起来即得到z值。其中的向量x是分类器的输入数据，向量w也就是我们要找到的最佳参数 （系数），从而使得分类器尽可能地精确。为了寻找该佳参数，需要用到优化理论的一些知识。
下面首先介绍梯度上升的优化方法，我们将学习到如何使用该方法求得数据集的佳参数。接下来，展示如何绘制梯度上升法产生的决策边界图，该图能将梯度上升法的分类效果可视化地呈现出来。后我们将学习随机梯度上升算法，以及如何对其进行修改以获得更好的结果。

3.1 梯度上升法

我们介绍的第一个优化算法叫做梯度上升法。梯度上升法基于的思想是：要找到某函数的大值，好的方法是沿着该函数的梯度方向探寻。如果梯度记为∇，则函数f(x,y)的梯度由下式表示：

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

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

该公式将一直被迭代执行，直至达到某个停止条件为止，比如迭代次数达到某个指定值或算法达到某个可以允许的误差范围。
我们最经常听到的应该是梯度下降算法，它与这里的梯度上升算法是一样的，只是公式中的加法需要变成减法。他是沿着下降最快的方向更新参数。

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

图中有100个样本点，每个点包含两个数值型特征：X1和X2。在此数据集上，我们将通过使用梯度上升法找到佳回归系数，也就是拟合出Logistic回归模型的佳参数。

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

注意：本书中省略了Logistic回归的损失函数和梯度计算推导过程，这里直接给出结论：
损失函数为：

梯度为真实值减去预测值的差乘以特征值：

具体的推导过程可以参考：机器学习实战教程（六）：Logistic回归基础篇之梯度上升算法
下面用代码实现梯度上升算法。
首先看一下提供的数据集：

# 看下数据情况
df =  pd.read_table("testSet.txt")
print(df.shape)
print(df.head())
# 结果如下
(99, 3)
   -0.017612  14.053064  0
0  -1.395634   4.662541  1
1  -0.752157   6.538620  0
2  -1.322371   7.152853  0
3   0.423363  11.054677  0
4   0.406704   7.067335  1
# 前面两列是特征，最后一列是分类结果

下面定义一个载入并处理数据的函数，它返回特征列表和标签列表：

def loadDataSet():
    '''
    将文本数据集转为列表
    '''
    dataMat = []
    labelMat = []
    fr = open('testSet.txt')
    for line in fr.readlines():
        lineArr = line.strip().split()
        # 第一列和第二列为特征
        # x0设为1，相当于z=w0+w1x1 + ... ,多了个常数项，有些是用b表示常数项
        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
        # 第三列为分类标签
        labelMat.append(int(lineArr[2]))
    return dataMat,labelMat

定义sigmoid函数：

def sigmoid(inX):
    '''
    定义sigmoid函数
    '''
    return 1.0/(1 + np.exp(-inX))

梯度上升算法，注意这里用到矩阵的运算：

def gradAscent(dataMatIn, classLabels):
    '''
    利用梯度上升法更新参数w，找到最优参数
    dataMatIn:特征的二维列表
    classLabels:标签列表
    '''
    # 把数据列表转为数组
#    dataMatrix = np.mat(dataMatIn)   #书中是转为矩阵matrix了
    dataMatrix = np.array(dataMatIn)   #维度(100, 3)
    # 得到特征数据的行和列
    m,n = dataMatrix.shape[0],dataMatrix.shape[1]
    #标签列表转为数组
#    labelMat = np.mat(classLabels).transpose()   #书中是转为矩阵matrix了
    labelMat = np.array(classLabels).T.reshape((m,1))   
    #特别注意要转换维度，一开始是(100,)，转换之后是(100,1),两个是不一样的,使用(100,)计算erro时会出现维度错误

    alpha = 0.001   # 初始化步长
    maxCycles = 500   # 设置最大更新次数
    # 初始化所有参数为1
    weights = np.ones((n,1))  # (3, 1)
    # 循环更新w
    for k in range(maxCycles):              
        h = sigmoid(np.dot(dataMatrix,weights))    #100x1
        error = (labelMat - h)   #100x1  
        weights = weights + alpha*np.dot(dataMatrix.T, error)  # 3x1
    return weights   #(3,1)

运行下看下实际效果：

# test1
dataArr,labelMat = loadDataSet()
weights = gradAscent(dataArr, labelMat)
print(weights)
# 结果如下
[[ 4.12414349]
 [ 0.48007329]
 [-0.6168482 ]]

3.3 分析数据：画出决策边界

上面已经解出了一组回归系数，它确定了不同类别数据之间的分隔线。那么怎样画出该分隔线，从而使得优化的过程便于理解呢？下面将解决这个问题。
画图思路：我们知道当sigmoid函数为0.5时，是分界线，此时，z则为0，即w₀ + w₁x₁ + w₂x₂=0。就可以以x₁为x轴，x₂为y轴画图了。

def plotBestFit(weights):
    '''
    可视化最优分类线
    '''
    import matplotlib.pyplot as plt
    dataMat,labelMat=loadDataSet()
    dataArr = np.array(dataMat)
    n = dataArr.shape[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 = np.arange(-3.0, 3.0, 0.1)   # (60,)
    y = (-weights[0]-weights[1]*x)/weights[2]
    ax.plot(x, y)
    plt.xlabel('X1'); plt.ylabel('X2');
    plt.show()

运行看一下效果：

# test2
dataArr,labelMat = loadDataSet()
weights =  gradAscent(dataArr, labelMat)
plotBestFit(weights)

这个分类结果相当不错。但是，尽管例子简单且数据集很小， 这个方法却需要大量的计算（300次乘法）。因此下一节将对该算法稍作改进，从而使它可以用在真实数据集上

3.4 训练算法：随机梯度上升

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

以下是随机梯度上升算法的实现代码。

    def stocGradAscent0(dataMatrix, classLabels):
    '''
    随机梯度下降法求最优解
    '''
    dataMatrix = np.array(dataMatrix)   #维度(100, 3)
    m,n = dataMatrix.shape[0],dataMatrix.shape[1]
    labelMat = np.array(classLabels).T.reshape((m,1))
    alpha = 0.01
    weights = np.ones(n)   #(3,)
    for i in range(m):
        h = sigmoid(np.dot(dataMatrix[i], weights))   #(1,3)*(3,)得到h为numpy.float64类型
        print(h.shape)
        error = labelMat[i] - h
        weights = weights + alpha * error * dataMatrix[i]
    return weights   #(3,)

可以看到，随机梯度上升算法与梯度上升算法在代码上很相似，但也有一些区别：后者的变量h和误差error都是矩阵，而前者则全是数值。
运行代码测试效果：

#test3
dataArr,labelMat = loadDataSet()
weights =  stocGradAscent0(dataArr, labelMat)
plotBestFit(weights)

可以看到，拟合出来的直线效果不太好。这里的分类器错分了三分之一的样本。
下图是随机梯度上升算法在200次迭代过程中回归系数的变化情况。运行随机梯度上升算法，在数据集的一次遍历中回归系数与迭代次数的关系图。回归系数经过大量迭代才能达到稳定值，并且仍然有局部的波动现象。

其中的系数2，也就是图中的X2只经过了50次迭代就达到了稳定值，但系数1和0则需要更多次的迭代。另外值得注意的是，在大的波动停止后，还有一些小的周期性波动。不难理解，产生这种现象的原因是存在一些不能正确分类的样本点（数据集并非线性可分），在每次迭代时会引发系数的剧烈改变。我们期望算法能避免来回波动，从而收敛到某个值。另外，收敛速度也需要加快。
为此做了改进，改进的随机梯度上升算法如下：

def stocGradAscent1(dataMatrix, classLabels, numIter=150):
    '''
    改进的随机梯度下降法
    numIter:默认的迭代次数
    '''
    dataMatrix = np.array(dataMatrix)
    m,n = dataMatrix.shape[0],dataMatrix.shape[1]
    labelMat = np.array(classLabels).T.reshape((m,1))
    weights = np.ones(n)
    for j in range(numIter):
        dataIndex = list(range(m))  #python3改动
        for i in range(m):
            alpha = 4/(1.0+j+i)+0.0001 
            randIndex = int(np.random.uniform(0,len(dataIndex)))
            h = sigmoid(sum(dataMatrix[randIndex]*weights))
            error = labelMat[randIndex] - h
            weights = weights + alpha * error * dataMatrix[randIndex]
            del(dataIndex[randIndex])
    return weights   #(3,)

第一处改进是：alpha在每次迭代的时候都会调整，这会缓解图上图的数据波动或者高频波动。另外，虽然alpha会随着迭代次数不断减小，但永远不会减小到0，这是因为中还存在一个常数项。必须这样做的原因是为了保证在多次迭代之后新数据仍然具有一定的影响。 如果要处理的问题是动态变化的，那么可以适当加大上述常数项，来确保新的值获得更大的回归系数。另一点值得注意的是，在降低alpha的函数中，alpha每次减少1/(j+i) ，其中j是迭代次数， i是样本点的下标。这样当j<<max(i)时，alpha就不是严格下降的。避免参数的严格下降也常见于模拟退火算法等其他优化算法中。
第二个改进是：通过随机选取样本来更新回归系数。这种方法将减少周期性的波动。
此外，改进算法还增加了一个迭代次数作为第3个参数。如果该参数没有给定的话，算法将默认迭代150次。如果给定，那么算法将按照新的参数值进行迭代。

上图是使用样本随机选择和alpha动态减少机制的随机梯度上升算法stocGradAscent1() 所生成的系数收敛示意图。该方法比采用固定alpha的方法收敛速度更快 。
下面看一下拟合效果：

#test4
dataArr,labelMat = loadDataSet()
weights =  stocGradAscent1(dataArr, labelMat)
plotBestFit(weights)  

4. 示例：从疝气病症预测病马的死亡率

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

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

数据中的缺失值是个非常棘手的问题，有很多文献都致力于解决这个问题。那么，数据缺失究竟带来了什么问题？假设有100个样本和20个特征，这些数据都是机器收集回来的。若机器上的某个传感器损坏导致一个特征无效时该怎么办？此时是否要扔掉整个数据？这种情况下，另外19个特征怎么办？它们是否还可用？答案是肯定的。因为有时候数据相当昂贵，扔掉和重新获取都是不可取的，所以必须采用一些方法来解决这个问题。
下面给出了一些可选的做法：

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

书中原始的数据集经过预处理之后保存成两个文件：horseColicTest.txt和horseColic- Training.txt。
看一下我们的数据集：

# 示例：从疝气病症预测马的死亡率
df =  pd.read_table("horseColicTraining.txt")
print(df.head())

得到数据为[5 rows x 22 columns]，可知实际只有有21个特征。

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

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

def classifyVector(inX, weights):
    '''
    计算输入的一个样本，sigmoid的值来预测是0还是1
    inX：一个特征向量
    weights：回归系数
    '''
    prob = sigmoid(sum(np.dot(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(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]))   # 提取特征
        # 如果预测值不等于真实值，错误数量加1
        if int(classifyVector(np.array(lineArr), trainWeights))!= int(currLine[21]):  
            errorCount += 1
    errorRate = (float(errorCount)/numTestVec)
    print ("the error rate of this test is: %f" % errorRate)
    return errorRate

为了多次拟合看下效果，我们定义一个函数：

def multiTest():
    '''
    多次分类预测
    '''
    numTests = 100   #预测次数
    errorSum=0.0   #错误率总和
    for k in range(numTests):
        errorSum += colicTest()
    print ("after %d iterations the average error rate is: %f" % (numTests, errorSum/float(numTests)))

运行看下效果：

#test5
multiTest()
# 运行结果如下
'after 100 iterations the average error rate is: 0.351045'

从上面的结果可以看到，10次迭代之后的平均错误率为35%。事实上，这个结果并不差，因为有30%的数据缺失。当然，如果调整colicTest()中的迭代次数和stochGradAscent1()中的步长，平均错误率可以降到20%左右。

总结

Logistic回归的目的是寻找一个非线性函数Sigmoid的佳拟合参数，求解过程可以由优化算法来完成。在优化算法中，常用的就是梯度上升算法，而梯度上升算法又可以简化为随机梯度上升算法。
随机梯度上升算法与梯度上升算法的效果相当，但占用更少的计算资源。此外，随机梯度上升是一个在线算法，它可以在新数据到来时就完成参数更新，而不需要重新读取整个数据集来进行批处理运算。
机器学习的一个重要问题就是如何处理缺失数据。这个问题没有标准答案，取决于实际应用中的需求。现有一些解决方案，每种方案都各有优缺点。