一、Logistic回归介绍

  回归算法有很多，比如：线性回归、Logistic回归、多项式回归、逐步回归、岭回归、Lasso回归等，Logistic回归是众多回归算法中的一员。我们常用Logistic回归模型做预测。通常，Logistic回归用于二分类问题，例如预测明天是否会下雨。
  1、Logistic回归
  假设现在有一些数据点，我们用一条直线对这些点进行拟合（该线称为最佳拟合直线），这个拟合过程就称作回归。利用Logistic回归进行分类的主要思想是：根据现有数据对分类边界线建立回归公式，以此进行分类。
  要想了解Logistic回归，我们必须先看一看Sigmoid函数 ，也可以称它为Logistic函数，它的值在[0,1]这个区间。它的公式如下：

  看下Sigmoid函数在不同坐标尺度下的两条曲线图。

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

  上两式整合成一个公式：

  z是一个矩阵，w是参数列向量(要求解的)，x是样本列向量(给定的数据集)，采用向量的写法可以写成 z = w^Tx ，它表示将这两个数值向量对应元素相乘然后全部加起来即得到z值。向量x是分类器的输入数据，向量w也就是我们要找到的最佳参数（系数），从而使得分类器尽可能地精确。
  如何得到合适的参数向量z？根据sigmoid函数的特性，我们可以做出如下的假设：

  式即为在已知样本x和参数w的情况下，样本x属性正样本(y=1)和负样本(y=0)的条件概率。理想状态下，根据上述公式，求出各个点的概率均为1，也就是完全分类都正确。但是考虑到实际情况，样本点的概率越接近于1，其分类效果越好。比如一个样本属于正样本的概率为0.51，那么我们就可以说明这个样本属于正样本。另一个样本属于正样本的概率为0.99，那么我们也可以说明这个样本属于正样本。但是显然，第二个样本概率更高，更具说服力。我们可以把上述两个概率公式合二为一，合并出来的Cost称之为代价函数(Cost Function)：

  当y等于1时，(1-y)项(第二项)为0；当y等于0时，y项(第一项)为0。为了简化问题，对整个Cost表达式求对数：

  这个代价函数，是对于一个样本而言的。给定一个样本，我们就可以通过这个代价函数求出，样本所属类别的概率，而这个概率越大越好，所以也就是求解这个代价函数的最大值。对于整体样本，假定各个样本与样本之间相互独立，那么根据极大似然估计法，整个样本集生成的概率即为所有样本生成概率的乘积，再将公式对数化，便可得到如下公式：

  m为样本的总数，y^{i表示第i个样本的类别，x}i表示第i个样本，w是多维向量，x(i)也是多维向量。满足J(w)的最大的w值即是我们需要求解的模型。
  求解使J(w)最大的w值可以使用梯度上升算法。如果面对的问题是求解使J(w)最小的w值，那么我们就需要使用梯度下降算法。面对我们这个问题，如果使J(θ) := -J(θ)，那么问题就从求极大值转换成求极小值了，使用的算法就从梯度上升算法变成了梯度下降算法，它们的思想都是相同的。本文使用梯度上升算法进行求解。
  2、梯度上升法
 &emsp梯度上升法基于的思想是：要找到某函数的最大值，最好的方法是沿着该函数的梯度方向探寻。如果梯度记为∇，则函数 f(x,y) 的梯度由下式表示：

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

  （该公式中的加法需要变成减法即为梯度下降算法）该公式将一直被迭代执行，直至达到某个停止条件为止，比如迭代次数达到某个指定值或算法达到某个可以允许的误差范围。

  同理，J(w)这个函数的极值，也可以这么求解。公式可以写为：

  由上面代价函数J(w)，Sigmoid函数：

  那么，只要求出J(θ)的偏导，就可以利用梯度上升算法求解J(θ)的极大值了。求解J(θ)对θ的偏导，求解如下：

  其中，为了简化显示，这里省略x的上标i，所以：

  因此，

  因此，梯度上升迭代公式为：

二、一个Logistic回归分类器的简单例子

贴上全部代码：

# -*- coding: utf-8 -*-
from cmath import exp
from matplotlib.font_manager import FontProperties
from numpy import *
import matplotlib.pyplot as plt
from sklearn.linear_model import LogisticRegression


def loadDataSet():
    """
    Function:
        加载数据集
    Parameters:
        无
    Returns:
        dataMat - 数据列表
        labelMat - 标签列表
    Modify:
        2018-08-27
    """
    dataMat = []
    labelMat = []
    fr = open('../Machine-Learning/machinelearninginaction/Ch05/testSet.txt')
    for line in fr.readlines():
        lineArr = line.strip().split()
        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
        labelMat.append(int(lineArr[2]))
    fr.close()
    return dataMat, labelMat


def plot_data_set():
    """
    Function:
        绘制数据集
    Parameters:
        无
    Returns:
        无
    Modify:
        2018-08-27
    """
    data_mat, label_mat = loadDataSet()
    data_arr = array(data_mat)
    n = shape(data_mat)[0]
    xcord_1 = []
    ycord_1 = []
    xcord_2 = []
    ycord_2 = []
    for i in range(n):
        if int(label_mat[i]) == 1:
            xcord_1.append(data_arr[i, 1])
            ycord_1.append(data_arr[i, 2])
        else:
            xcord_2.append(data_arr[i, 1])
            ycord_2.append(data_arr[i, 2])
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(xcord_1, ycord_1, s=20, c='red', marker='s', alpha=0.5)
    ax.scatter(xcord_2, ycord_2, s=20, c='green', marker='s', alpha=0.5)
    plt.title('data_set')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.show()


def sigmoid(inX):
    """
    Function:
        Sigmoid函数
    Parameters:
        inX - 数据
    Returns:
        Sigmoid函数
    Modify:
        2018-08-27
    """
    return 1.0 / (1 + exp(-inX))


def gradAscent(dataMatIn, classLabels):
    """
    Function:
        梯度上升算法
    Parameters:
        dataMatIn - 数据集
        classLabels - 类别标签
    Returns:
        weights - 求得的权重，最优参数
    Modify:
        2018-08-27
    """
    # dataMatrix 100 * 3
    dataMatrix = mat(dataMatIn)
    # labMatrix 100 * 1
    labMatrix = mat(classLabels).transpose()
    m, n = shape(dataMatrix)
    alpha = 0.01
    maxCycles = 500
    weights = ones((n, 1))
    weights_array = array([])
    for k in range(maxCycles):
        # h 100 * 1
        h = sigmoid(dataMatrix * weights)
        error = (labMatrix - h)
        # weights 3 * 1
        weights = weights + alpha * dataMatrix.transpose() * error
        weights_array = append(weights_array, weights)
    weights_array = weights_array.reshape(maxCycles, n)
    # return weights.getA()
    return weights.getA(), weights_array


def poltBestFit(weights):
    """
    Function:
        画出数据集和Logistic回归最佳拟合直线函数
    Parameters:
        weights - 权重，最优参数
    Returns:
        无
    Modify:
        2018-08-28
    """
    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(-4.0, 4.0, 0.1)
    y = (-weights[0] - weights[1] * x) / weights[2]
    ax.plot(x, y)
    plt.xlabel('X1');
    plt.ylabel('X2')
    plt.show()


def stocGradAscent0(dataMatrix, classLabel):
    """
    Function:
        随机梯度上升算法
    Parameters:
        dataMatrix - 数据集
        classLabels - 类别标签
    Returns:
        weights - 最优参数
    Modify:
        2018-08-27
    """
    m, n = shape(dataMatrix)
    alpha = 0.01
    weights = ones(n)
    for i in range(m):
        # h 1 * 1
        h = sigmoid(sum(dataMatrix[i] * weights))
        error = classLabel[i] - h
        # weights 3 * 1
        weights = weights + alpha * error * dataMatrix[i]
    return weights


def stocGradAscent1(dataMat, classLabel, numIter=150):
    """
    Function:
        改进的随机梯度上升算法
    Parameters:
        dataMatrix - 数据集
        classLabels - 类别标签
        numIter - 迭代次数
    Returns:
        weights - 最优参数
    Modify:
        2018-08-28
    """
    m, n = shape(dataMat)
    weights = ones(n)
    weights_array = array([])
    for i in range(numIter):
        dataIndex = list(range(m))
        for j in range(m):
            alpha = 4 / (1.0 + j + i) + 0.01
            randIndex = int(random.uniform(0, len(dataIndex)))
            h = sigmoid(sum(dataMat[randIndex] * weights))
            error = classLabel[randIndex] - h
            weights = weights + alpha * error * dataMat[randIndex]
            weights_array = append(weights_array, weights, axis=0)
            del (dataIndex[randIndex])
    weights_array = weights_array.reshape(numIter * m, n)
    # return weights
    return weights, weights_array


def plotWeights(weights_array1, weights_array2):
    """
    Function:
        绘制回归系数与迭代次数关系
    Parameters:
        weights_array1 - 梯度上升算法回归系数
        weights_array2 - 改进的随机梯度上升算法回归系数
    Returns:
        无
    Modify:
        2018-08-28
    """
    # 设置汉字格式
    font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=14)
    # fig画布被分为六个区域,axs[0][0]表示第一行第一列
    fig, axs = plt.subplots(nrows=3, ncols=2, sharex=False, sharey=False, figsize=(20, 10))
    x1 = arange(0, len(weights_array1), 1)
    # 绘制w0与迭代次数的关系
    axs[0][0].plot(x1, weights_array1[:, 0])
    axs0_title_text = axs[0][0].set_title(u'梯度上升算法：回归系数与迭代次数关系', FontProperties=font)
    axs0_ylabel_text = axs[0][0].set_ylabel(u'W0', FontProperties=font)
    plt.setp(axs0_title_text, size=20, weight='bold', color='black')
    plt.setp(axs0_ylabel_text, size=20, weight='bold', color='black')
    # 绘制w1与迭代次数的关系
    axs[1][0].plot(x1, weights_array1[:, 1])
    axs1_ylabel_text = axs[1][0].set_ylabel(u'W1', FontProperties=font)
    plt.setp(axs1_ylabel_text, size=20, weight='bold', color='black')
    # 绘制w2与迭代次数的关系
    axs[2][0].plot(x1, weights_array1[:, 2])
    axs2_xlabel_text = axs[2][0].set_xlabel(u'迭代次数', FontProperties=font)
    axs2_ylabel_text = axs[2][0].set_ylabel(u'W2', FontProperties=font)
    plt.setp(axs2_xlabel_text, size=20, weight='bold', color='black')
    plt.setp(axs2_ylabel_text, size=20, weight='bold', color='black')

    x2 = arange(0, len(weights_array2), 1)
    # 绘制w0与迭代次数的关系
    axs[0][1].plot(x2, weights_array2[:, 0])
    axs0_title_text = axs[0][1].set_title(u'改进的随机梯度上升算法：回归系数与迭代次数关系', FontProperties=font)
    axs0_ylabel_text = axs[0][1].set_ylabel(u'W0', FontProperties=font)
    plt.setp(axs0_title_text, size=20, weight='bold', color='black')
    plt.setp(axs0_ylabel_text, size=20, weight='bold', color='black')
    # 绘制w1与迭代次数的关系
    axs[1][1].plot(x2, weights_array2[:, 1])
    axs1_ylabel_text = axs[1][1].set_ylabel(u'W1', FontProperties=font)
    plt.setp(axs1_ylabel_text, size=20, weight='bold', color='black')
    # 绘制w2与迭代次数的关系
    axs[2][1].plot(x2, weights_array2[:, 2])
    axs2_xlabel_text = axs[2][1].set_xlabel(u'迭代次数', FontProperties=font)
    axs2_ylabel_text = axs[2][1].set_ylabel(u'W1', FontProperties=font)
    plt.setp(axs2_xlabel_text, size=20, weight='bold', color='black')
    plt.setp(axs2_ylabel_text, size=20, weight='bold', color='black')

    plt.show()


def classifyVector(inX, weights):
    """
    Function:
        Logistic回归分类函数
    Parameters:
        inX - 特征向量
        weights - 回归系数
    Returns:
        分类结果
    Modify:
        2018-08-29
    """
    prob = sigmoid(sum(inX * weights))
    if prob > 0.5:
        return 1.0
    else:
        return 0.0


def colicTest():
    """
    Function:
        Logistic回归分类测试函数
    Parameters:
        无
    Returns:
        分类结果错误率
    Modify:
        2018-08-29
    """
    frTrain = open('./machinelearninginaction/Ch05/horseColicTraining.txt')
    frTest = open('./machinelearninginaction/Ch05/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, 500)
    # 使用梯度上升算法
    trainWeights = gradAscent(array(trainingSet), trainingLabels)
    errorCount = 0.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]):
        # 使用梯度上升算法
        if int(classifyVector(array(lineArr), trainWeights[:, 0])) != int(currLine[-1]):
            errorCount += 1
    errorRate = (float(errorCount) / numTestVec)
    print('the error rate of this test is: %f' % errorRate)
    return errorRate


def multiTest():
    """
    Function:
        多次调用colicTest()函数，求结果的平均值
    Parameters:
        无
    Returns:
        分类结果平均错误率
    Modify:
        2018-08-29
    """
    numTests = 10;
    errorSum = 0.0
    for k in range(numTests):
        errorSum += colicTest()
    print("after %d iterations the average error rate is: %f" % (numTests, errorSum / float(numTests)))


def colic_sklearn():
    fr_train = open('./machinelearninginaction/Ch05/horseColicTraining.txt')
    fr_test = open('./machinelearninginaction/Ch05/horseColicTest.txt')
    training_set = []
    training_labels = []
    test_set = []
    test_labels = []
    for line in fr_train.readlines():
        curr_line = line.strip().split('\t')
        line_arr = []
        for i in range(len(curr_line) - 1):
            line_arr.append(float(curr_line[i]))
        training_set.append(line_arr)
        training_labels.append(float(curr_line[-1]))
    for line in fr_test.readlines():
        curr_line = line.strip().split('\t')
        line_arr = []
        for i in range(len(curr_line) - 1):
            line_arr.append(float(curr_line[i]))
            test_set.append(line_arr)
            test_labels.append(float(curr_line[-1]))
    classifier = LogisticRegression(solver='sag', max_iter=5000).fit(training_set, training_labels)
    test_accuracy = classifier.score(test_set, test_labels) * 100
    print('正确率：%f%%' % test_accuracy)


if __name__ == '__main__':
    dataArr, labelMat = loadDataSet()
    # print(dataArr)
    # print(labelMat)
    # plot_data_set()

    # print(gradAscent(dataArr, labelMat))

    # weights = gradAscent(dataArr, labelMat)
    # poltBestFit(weights)

    # weights = stocGradAscent0(array(dataArr), labelMat)
    # poltBestFit(weights)

    # weights = stocGradAscent1(array(dataArr), labelMat, 500)
    # poltBestFit(weights)

    # weights1, weights_array1 = stocGradAscent1(array(dataArr), labelMat)
    # weights2, weights_array2 = gradAscent(dataArr, labelMat)
    # plotWeights(weights_array1, weights_array2)

    # multiTest()

    colic_sklearn()

  1、加载数据集
  这个数据集有两列特征，最后一列数据即为分类标签。

  我们可以将第一列数据(X1)看作x轴上的值，第二列数据(X2)看作y轴上的值，最后一列数据即为分类标签，画出他们的分布。

  2、训练算法
  梯度上升迭代公式：

  将上述公式矢量化：

  截图即为求解出的回归系数[w0,w1,w2]。

  3、画出决策边界

  这个分类结果相当不错，从图上看只错分了两到四个点。但是，尽管例子简单且数据集很小，这个方法却需要大量的计算（300次乘法）。因此下一节将对该算法稍作改进，从而使它可以用在真实数据集上。
  4、随机梯度上升算法
  梯度上升算法在每次更新回归系数时都需要遍历整个数据集，该方法在处理100个左右的数据集时尚可，但如果有数十亿样本和成千上万的特征，那么该方法的计算复杂度就太高了。一种改进方法是一次仅用一个样本点来更新回归系数，该方法称为随机梯度上升算法。由于可以在新样本到来时对分类器进行增量式更新，因而随机梯度上升算法是一个在线学习算法。

  随机梯度上升算法与梯度上升算法在代码上很相似，但也有一些区别：第一，后者的变量 h 和误差 error 都是向量，而前者则全是数值；第二，前者没有矩阵的转换过程，所有变量的数据类型都是NumPy数组。从画出的划分结果来看，这里的分类器错分了三分之一的样本，结果并不完美。
  5、改进的梯度上升算法
  该算法第一个改进之处在于，alpha在每次迭代的时候都会调整，虽然alpha会随着迭代次数不断减小，但永远不会减小到0，因为这里还存在一个常数项。必须这样做的原因是为了保证在多次迭代之后新数据仍然具有一定的影响。如果需要处理的问题是动态变化的，那么可以适当加大上述常数项，来确保新的值获得更大的回归系数。另一点值得注意的是，在降低alpha的函数中，alpha每次减少1/(j+i)，其中j是迭代数，i是样本点的下标。
  第二个改进的地方在于更新回归系数(最优参数)时，只使用一个样本点，并且选择的样本点是随机的，每次迭代不使用已经用过的样本点。这样的方法，就有效地减少了计算量，并保证了回归效果。

  从上面的结果分布图。该分隔线达到了与GradientAscent()差不多的效果，但是所使用的计算量更少。默认迭代次数是150，可以通过 stocGradAscent() 的第3个参数来对此进行修改，例如修改为500次：

下面看看，回归系数与迭代次数的关系。

由于改进的随机梯度上升算法随机选取样本点，所以每次的运行结果是不同的，但是大体趋势是一样的。图中看出改进的随机梯度上升算法收敛效果更好。

一共有100个样本点，改进的随机梯度上升算法迭代次数为150。而上图显示15000次迭代次数的原因是，使用一次样本就更新一下回归系数。因此，迭代150次，相当于更新回归系数150*100=15000次。从上图左侧的改进随机梯度上升算法回归效果中可以看出，其实在更新2000次回归系数的时候，已经收敛了。相当于遍历整个数据集20次的时候，回归系数已收敛。训练已完成。

再看看上图右侧的梯度上升算法回归效果，梯度上升算法每次更新回归系数都要遍历整个数据集。从图中可以看出，当迭代次数为300多次的时候，回归系数才收敛。

三、示例：从疝气病症预测病马死亡率

数据集-训练集部分截图

分析过程

  这里的数据包含了368个样本和28个特征。这种病不一定源自马的肠胃问题，其他问题也可能引发马疝病。该数据集中包含了医院检测马疝病的一些指标，有的指标比较主观，有的指标难以测量，例如马的疼痛级别。另外需要说明的是，除了部分指标主观和难以测量外，该数据还存在一个问题，数据集中有30%的值是缺失的。下面将首先介绍如何处理数据集中的数据缺失问题，然后再利用Logistic回归和随机梯度上升算法来预测病马的生死。
  数据中的缺失值是一个非常棘手的问题，很多文献都致力于解决这个问题。那么，数据缺失究竟带来了什么问题？假设有100个样本和20个特征，这些数据都是机器收集回来的。若机器上的某个传感器损坏导致一个特征无效时该怎么办？它们是否还可用？答案是肯定的。因为有时候数据相当昂贵，扔掉和重新获取都是不可取的，所以必须采用一些方法来解决这个问题。下面给出了一些可选的做法：

  预处理数据做两件事：
   第一，如果测试集中一条数据的特征值已经确实，那么我们选择实数0来替换所有缺失值，因为本文使用Logistic回归。因此这样做不会影响回归系数的值。sigmoid(0)=0.5，即它对结果的预测不具有任何倾向性。
  第二，如果测试集中一条数据的类别标签已经缺失，那么我们将该类别数据丢弃，因为类别标签与特征不同，很难确定采用某个合适的值来替换。
  原始的数据集经过处理，保存为两个文件：horseColicTest.txt和horseColicTraining.txt。
  使用Logistic回归方法进行分类并不需要做很多工作，所需做的只是把测试集上每个特征向量乘以最优化方法得来的回归系数，再将该乘积结果求和，最后输入到Sigmoid函数中即可。如果对应的Sigmoid值大于0.5就预测类别标签为1，否则为0。

使用改进的随机梯度上升算法

使用梯度上升算法

  这个数据集使用梯度上升算法执行速度比改进的梯度上升法快很多，而且从运行结果对比来看，使用梯度上升算法错误率稳定且较低。
  所以我们可以得到如下结论：
  当数据集较小时，我们使用梯度上升算法。
  当数据集较大时，我们使用改进的随机梯度上升算法。

四、应用scikit-learn构建Logistic回归分类器

from sklearn.linear_model import LogisticRegression


def colic_sklearn():
    fr_train = open('./machinelearninginaction/Ch05/horseColicTraining.txt')
    fr_test = open('./machinelearninginaction/Ch05/horseColicTest.txt')
    training_set = []
    training_labels = []
    test_set = []
    test_labels = []
    for line in fr_train.readlines():
        curr_line = line.strip().split('\t')
        line_arr = []
        for i in range(len(curr_line) - 1):
            line_arr.append(float(curr_line[i]))
        training_set.append(line_arr)
        training_labels.append(float(curr_line[-1]))
    for line in fr_test.readlines():
        curr_line = line.strip().split('\t')
        line_arr = []
        for i in range(len(curr_line) - 1):
            line_arr.append(float(curr_line[i]))
            test_set.append(line_arr)
            test_labels.append(float(curr_line[-1]))
    classifier = LogisticRegression(solver='liblinear', max_iter=10).fit(training_set, training_labels)
    test_accuracy = classifier.score(test_set, test_labels) * 100
    print('正确率：%f%%' % test_accuracy)


if __name__ == '__main__':
    colic_sklearn()

solver='liblinear', max_iter=10

  可以看到，正确率又高一些了。更改solver参数，设置为sag，使用随机平均梯度下降算法，看一看效果。你会发现，有警告了。

solver='sag', max_iter=10'

  显而易见，警告是因为算法还没有收敛。更改max_iter=5000，再运行代码：

solver='sag', max_iter=5000

  迭代收敛后已经没有警告了。

五、小结

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

随机梯度上升算法与梯度上升算法的效果相当，但占用更少的计算资源。此外，随机梯度上升是一个在线算法，它可以在新数据到来时就完成参数更新，而不需要重新读取整个数据集来进行批处理运算。

优缺点

总结