机器学习实战学习笔记（四）Logistic回归

PS：该系列数据都可以在图灵社区（点击此链接）中随书下载中下载（如下）

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

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

  我们想要的函数应该是，能接受所有的输入然后预测出类别。例如在两个类的情况下，函数输出0或1，这种性质的函数称为海维赛德阶跃函数（Heaviside step function），或者直接称为单位阶跃函数。然而，海维赛德阶跃函数的问题在于：该函数在跳跃点上从0瞬间跳跃到1，这个瞬间跳跃过程有时很难处理。找一个类似性质的函数，且数学上更易处理，这就是Sigmoid函数。下面给出Sigmoid函数具体的计算公式和图像：
$$\sigma (z)=\frac{1}{1+e^{-z}}$$

  由上图可知，当x为0时，Sigmoid函数值为0.5。随着x的增大，对应的Sigmoid值将逼近于1；而随着x的减小，Sigmoid值将逼近于0。因此，为了实现Logistic回归分类器，我们可以在每个特征上都乘以一个回归系数，然后把所有的结果值相加，将这个总和带入Sigmoid函数中，进而得到一个范围在0~1之间的数值。任何大于0.5的数据被分入1类，小于0.5即被归入0类。所以Logistic回归也可以被看成是一种概率估计。

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

  Sigmoid函数的输入记为$z$，$z=w_0{x}_0+w_1{x}_1+w_2{x}_2+...+w_n{x}_n$，如果采用向量的写法，上述公式可以写成$z=w^{T}x$，他表示将这两个数值向量对应元素相乘然后全部加起来得到$z$值。其中向量$x$是分类器的输入数据，向量$w$也就是我们要找到的最佳参数（系数）。

2.1 梯度上升法

  梯度上升法基于的思想是：要找到某函数的最大值，最好的方法是沿着该函数的梯度方向探寻。如果梯度记为$\nabla$，则函数$f(x,y)$的梯度为：
$$\nabla f(x,y)=\left(\begin{array}{l}\frac{\partial f(x,y)}{\partial x}\\ \frac{\partial f(x,y)}{\partial y}\end{array}\right)$$
  这个梯度意味着要沿$x$的方向移动$\frac{\partial f(x,y)}{\partial x}$，沿$y$的方向移动$\frac{\partial f(x,y)}{\partial y}$。其中函数$f(x,y)$必须要在待计算的点上有定义并且可微。梯度算子总是指向函数值增长最快的方向，而移动量的大小称为步长，记做$\alpha$。用向量来表示的话，梯度算法的迭代公式为：$w:=w+\alpha {\nabla }_{w}f(w)$。该公式将一直被迭代执行，直到达到某个停止条件为止，比如迭代次数达到某个指定值或算法达到某个可以允许的误差范围。
  梯度下降算法与梯度上升算法时一样的，公式变为：$w:=w-\alpha {\nabla }_{w}f(w)$前者用来求最小值，后者用来求最大值。

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

#伪代码：
每个回归系数初始化为1
重复R次：
	计算整个数据集的梯度
	使用alpha×gradient
	返回回归系数

  创建一个名为logRegres.py的文件，并下载数据集，将其中的testSet.txt文件复制到logRegres.py所在的文件夹，输入如下代码：

import numpy as np

def  loadDataSet():
    '''创建数据集'''
    dataMat = []; labelMat = []
    with open('testSet.txt', 'r') as fileObject:
        for line in fileObject.readlines():
            lineArr = line.strip().split()
            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 + np.exp(-inX))

def gradAscent(dataMatIn, classLabels):
    '''梯度上升'''
    #转换为NumPy矩阵数据类型
    dataMatrix = np.mat(dataMatIn)
    labelMat = np.mat(classLabels).transpose()
    m, n = np.shape(dataMatrix)
    alpha = 0.001
    maxCycles = 500
    weights = np.ones((n, 1))
    for k in range(maxCycles):
        #梯度更新，定性地说是计算真实类别与预测类别的差值
        h = sigmoid(dataMatrix * weights)
        error = (labelMat - h)
        weights = weights + alpha * dataMatrix.transpose() * error
    return weights

  便利函数loadDataSet()，它的主要功能是打开文本文件testSet.txt文件并逐行读取。每行前两个值分别是X1和X2，第三个值是数据对应的类别标签，为了方便计算，该函数还将X0的值设为1.0。我们利用NumPy矩阵对数据进行处理，这里需要详细说明的是梯度更新的部分的两行更新公式。下面进入数学推导：

2.2.1 Logistic回归梯度更新公式推导

  我们记Sigmoid函数为$$g(z)=\frac{1}{1+e^{-z}}$$
  对Sigmoid函数求导很容易得到$$g^{\prime }(z)=g(z)(1-g(z))$$
  逻辑回归的假设函数形式如下：$$h_{\theta }(x)=g({\theta }^{T}x),g(z)=\frac{1}{1+e^{-z}}$$
  所以:$$h_{\theta }(x)=\frac{1}{1+e^{-{\theta }^{T}x}}$$
  其中$x$是我们的输入，$\theta$为我们要求取的参数。一个机器学习的模型，实际上是把决策函数限定在某一组条件下，这组条件就决定了模型的假设空间。而逻辑回归模型所做的假设是：$$P(y=1|x;\theta )=g({\theta }^{T}x)=\frac{1}{1+e^{-{\theta }^{T}x}}=h_{\theta }(x)$$
  这个函数的意思就是在给定$x$和$\theta$的条件下$y=1$的概率，基于假设显然我们可以得到：$$P(y=0|x;\theta )=1-h_{\theta }x,P(y|x;\theta )=(h_{\theta }(x){)}^y(1-h_{\theta }(x){)}^{1-y}$$
  假设样本与样本之间相互独立，那么整个样本集生成的概率即为所有样本生成概率的乘积：
$$\begin{array}{rl}L(\theta )& =p(y|X;\theta )\\ & =\prod _{i=1}^{m}p\left(y^{(i)}|x^{(i)};\theta \right)\\ & =\prod _{i=1}^m{\left(h_{\theta }\left(x^{(i)}\right)\right)}^{y^{(i)}}{\left(1-h_{\theta }\left(x^{(i)}\right)\right)}^{1-y^{(i)}}\end{array}$$
  通过log转换以便于求导：
$$\begin{array}{rl}\ell (\theta )& =\log L(\theta )\\ & =\sum _{i=1}^m{y}^{(i)}\log h\left(x^{(i)}\right)+\left(1-y^{(i)}\right)\log \left(1-h\left(x^{(i)}\right)\right)\end{array}$$
  对上述式子求导：
$$\begin{array}{rl}\frac{\partial }{\partial {\theta }_j}(\ell (\theta ))& =\frac{\partial }{\partial {\theta }_j}\left(\sum _{i=1}^m{y}^{(i)}\log \left(h\left(x^{(i)}\right)\right)+\left(1-y^{(i)}\right)\log \left(1-h\left(x^{(i)}\right)\right)\right)\\ & =\left(\frac{y^{(i)}}{h\left(x^{(i)}\right)}-\left(1-y^{(i)}\right)\frac{1}{1-h\left(x^{(i)}\right)}\right)\frac{\partial }{\partial {\theta }_j}\left(h\left(x^{(i)}\right)\right)\\ & =\left(\frac{y^{(i)}}{g\left({\theta }^T{x}^{(i)}\right)}-(1-y^{(i)})\frac{1}{1-g\left({\theta }^T{x}^{(i)}\right)}\right)\frac{\partial }{\partial {\theta }_j}\left(g\left({\theta }^T{x}^{(i)}\right)\right)\\ & =\left(\frac{y^{(i)}}{g\left({\theta }^T{x}^{(i)}\right)}-(1-y^{(i)})\frac{1}{1-g\left({\theta }^T{x}^{(i)}\right)}\right)g\left({\theta }^T{x}^{(i)}\right)\left(1-g\left({\theta }^T{x}^{(i)}\right)\right)\frac{\partial {\theta }^T{x}^{(i)}}{\partial {\theta }_j}\\ & =\left(y^{(i)}\left(1-g\left({\theta }^T{x}^{(i)}\right)-\left(1-y^{(i)}\right)g\left({\theta }^T{x}^{(i)}\right)\right)\right)x_j\\ & =\left(y^{(i)}-h_{\theta }\left(x^{(i)}\right)\right)x_j\end{array}$$
  所以我们得到梯度更新的公式为$$w:=w+\alpha \ast gradient,gradient=(y_i-sigmoid(w^{T}x))\ast x$$
  现在回过头去看梯度更新中的那两句更新公式代码就可以看懂了。然后我们在Python命令行进行测试：

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

def plotBestFit(wei):
    '''绘制数据集和Logistic回归最佳拟合直线'''
    import matplotlib.pyplot as plt
    weights = np.array(wei)
    dataMat, labelMat = loadDataSet()
    dataArr = np.array(dataMat)
    n = np.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 = np.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()

  代码里设置了sigmoid函数参数为0,0是两类的分解除，所以设定$0=w_0{x}_0+w_1{x}_1+w_2{x}_2$
  测试运行结果：


  从图上看分类结果只分错了两到四个点，但是该方法需要大量的计算。

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

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

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

def stocGradAscent0(dataMatrix, classLabels):
    '''随机梯度上升算法'''
    dataArr = np.array(dataMatrix)
    m, n = np.shape(dataArr)
    alpha = 0.01
    weights = np.ones(n)
    for i in range(m):
        h = sigmoid(sum(dataArr[i] * weights))
        error = classLabels[i] - h
        weights = weights + alpha * error * dataArr[i]
    return weights

  显然随机梯度算法在这里错分了三分之一左右的样本，所以我们需要优化随机梯度算法。

def stocGradAscent1(dataMatrix, classLabels, numIter=150):
    '''改进的随机梯度上升算法'''
    dataArr = np.array(dataMatrix)
    m, n = np.shape(dataArr)
    weights = np.ones(n)
    for j in range(numIter):
        dataIndex = range(m)
        for i in range(m):
            #alpha每次迭代需要调整
            alpha = 4.0 / (1.0 + j + i) + 0.01
            #随机选取更新
            randIndex = int(np.random.uniform(0, len(dataIndex)))
            h = sigmoid(sum(dataArr[randIndex] * weights))
            error = classLabels[randIndex] - h
            weights = weights + alpha * error * dataArr[randIndex]
            del(list(dataIndex)[randIndex])
    return weights

  (1)alpha在每次迭代的时候都会调整，以缓解数据波动或高频波动；（2）随机选取样本来更新回归系数，以减少周期性的波动。对结果进行测试发现达到了更好的效果，当然还可以通过修改迭代次数步长等参数来达到更好的分类效果。

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

3.1 处理数据：处理数据中的缺失值

• 使用可用特征的均值来填补缺失值
• 使用特殊值来填补缺失值，如-1
• 忽略有缺失值的样本
• 使用相似样本的均值填补缺失值
• 使用另外的机器学习算法预测缺失值
    上述是可选的处理数据中的缺失值的方法，而我们在数据预处理阶段使用一个实数值0来替换缺失值，因为回归系数更新公式：$weights=weights+alpha\ast error\ast dataMatrix[randIndex]$，我们选择0填补某个特征值，那么该特征值的系数将不会更新，即$weights=weights$；另外由于$sigmoid(0)=0.5$，即它对结果的预测不具有任何倾向性，因此上述做法也不会对误差项造成任何影响。
    如果一条数据的类别标签已经确实，那么我们简单的做法是丢弃该条数据。处理玩的数据保存在文件：horseColicTest.txt和horseColicTraining.txt中，并将其复制到logRegres.py所在文件夹。

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

def classifyVector(inX, weights):
    '''Logistic分类函数'''
    prob = sigmoid(sum(inX * weights))
    if prob > 0.5:
        return 1.0
    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, 500)
    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(np.array(lineArr), trainWeights)) != int(currLine[21]):
            errorCount += 1
    errorRate = float(errorCount) / numTestVec
    print("the error rate of this test is: ", errorRate)
    return errorRate

def multiTest():
    '''多组测试'''
    numTests = 10; errorSum = 0.0
    for k in range(numTests):
        errorSum += colicTest()
    print("after ", numTests, " iterations the average error rate is: ", errorSum / float(numTests))

  测试结果如下：

  从结果中可以看到10次迭代的平均错误率为35%，这是由于数据缺失导致的，当然可以通过调整一次分类的迭代次数和步长降低错误率。

4 小结

  Logistic回归的目的是寻找一个非线性函数Sigmoid的最佳拟合参数，求解过程可以由最优化算法来完成。在最优化算法中，最常用的就是梯度上升算法，而梯度上升算法又可以简化为随机梯度上升算法。
  随机梯度上升算法与梯度上升算法的效果相当，但占用更少的计算资源。此外，随机梯度算法是一个在线算法，它可以在新数据到来时就完成参数更新，而不需要重新读取整个数据集来进行批处理运算。