《机器学习》及实战五、Logistic回归理论及实战

• Python版本： Python3.x
• 运行平台： Windows
• IDE： PyCharm
• 参考资料：《机器学习》（西瓜书）《机器学习实战》（王斌）
• 转载请标明出处：https://blog.csdn.net/tian121381/category_9748511.html
• 资料下载，提取码：719m
  

一 、前言

这⼀章，我们将首次接触到最优化算法。仔细想想就会发现，其实我们日常生活中遇到过很多最优化问题，比如如何在最短时间内从A点到达B点？如何投入最少工作量却获得最⼤的效益？ 如何设计发动机使得油耗最少而功率最大？可见，最优化的作用十分强大。接下来，我们介绍几个最优化算法，并利用它们训练出⼀个非线性函数用于分类。

二、 Logistic回归和Sigmoid函数的分类

1. Logistic回归

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

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

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^Tx。它表示将这两个数值向量的对应元素相乘然后全部加起来即得到z值。其中的向量x是分类器的输入数据，向量w也就是我们要找到的最佳参数（系数），从而使得分类尽可能地精确。为了寻找该最佳参数，需要用到最优化理论的⼀些知识。 下面首先介绍梯度上升的最优化方法（和梯度下降一样的），我们将学习到如何使用该⽅法求得数据集的最佳参数。

4. 梯度上升

梯度上升法基于的思想是：要找到某函数的最⼤值，最好的⽅法是沿着该函数的梯度方向探寻（因为梯度的方向是变化最大的方向）。如果梯度记为 ∇，则函数f(x,y)的梯度由下式表示：

举个例子哈，如下图。

梯度上升算法到达每个点后都会重新估计移动的方向。从P0开始，计算完该点的梯度，函数就根据梯度移动到下⼀点P1。在P1点，梯度再次被重新计算，并沿新的梯度⽅向移动到P2。如此循环迭代，直到满足停止条件。迭代的过程中，梯度算子总是保证我们能选取到最佳的移动方向。
图中的梯度上升算法沿梯度⽅向移动了⼀步。可以看到，梯度算子总是指向函数值增长最快的⽅向。这里所说的是移动⽅向，而未提到移动量的大小。该量值称为步长，记做 。用向量来表示的话，梯度上升算法的迭代公式如下（下面会给出代码例子的）：

该公式将⼀直被迭代执行，直至达到某个停止条件为止，比如迭代次数达到某个指定值或算法达到某个可以允许的误差范围。
梯度下降算法你最经常听到的应该是梯度下降算法，它与这里的梯度上升算法是⼀样的，只是公式中的加法需要变成减法。因此，对应的公式可以写成

梯度上升算法用来求函数的最⼤值，而梯度下降算法用来求函数的最小值。
现在，我们理论的东西讲的差不多了，开始进入代码中深入理解吧！

三、 Logistic回归分类器的简单应用

现在，我们有一个数据集（置顶资料下载处自行下载），我们要将采用梯度上升法找到 Logistic回归分类器在此数据集上的最佳回归系数。
先看看这个数据集（部分）什么样的吧！

可以看出，他有两个特征（前两列），1个标签，分了0，1两种情况。这时我们可以直接想到用二维坐标系能清楚的表现出他们的关系呀。
用代码实现一下！

1. 数据处理

引入相关包

"""
@Author:Yuuuuu、Tian
Thu Feb 20 16:28:58 2020
背景：使用梯度上升法找的最佳参数
    有100个样本，每个样本点包含两个数值型特征：X1,X2。在此数据集上，我们通过使用梯度上升法找到最佳的回归系数，也就是拟合出Logistic回归模型的最佳参数
"""
from matplotlib.font_manager import FontProperties
import numpy as np
import matplotlib.pyplot as plt
import math
import random

读取文件，提取数据

def loadDataSet():  #加载数据
    dataMat = []   #特征列表
    labelMat = []  #标签列表
    fr = open('0701 DataSet/testSet.txt')
    for line in fr.readlines():
        lineArr = line.strip().split()   #去回车，去空格
        dataMat.append([1.0,float(lineArr[0]), float(lineArr[1])])    #第一个1为偏移值，相当于b
        #print('a',dataMat)
        labelMat.append(int(lineArr[2]))
        #print('b',labelMat)
    fr.close()
    return dataMat,labelMat
#-------测试--------------------
dataMat,labelMat = loadDataSet()

结果：

[[1.0, -0.017612, 14.053064], [1.0, -1.395634, 4.662541], [1.0, -0.752157, 6.53862], [1.0, -1.322371, 7.152853], [1.0, 0.423363, 11.054677], [1.0, 0.406704, 7.067335], [1.0, 0.667394, 12.741452], [1.0, -2.46015, 6.866805], [1.0, 0.569411, 9.548755], [1.0, -0.026632, 10.427743], [1.0, 0.850433, 6.920334], [1.0, 1.347183, 13.1755], [1.0, 1.176813, 3.16702], [1.0, -1.781871, 9.097953], [1.0, -0.566606, 5.749003], [1.0, 0.931635, 1.589505], [1.0, -0.024205, 6.151823], [1.0, -0.036453, 2.690988], [1.0, -0.196949, 0.444165], [1.0, 1.014459, 5.754399], [1.0, 1.985298, 3.230619], [1.0, -1.693453, -0.55754], [1.0, -0.576525, 11.778922], [1.0, -0.346811, -1.67873], [1.0, -2.124484, 2.672471], [1.0, 1.217916, 9.597015], [1.0, -0.733928, 9.098687], [1.0, -3.642001, -1.618087], [1.0, 0.315985, 3.523953], [1.0, 1.416614, 9.619232], [1.0, -0.386323, 3.989286], [1.0, 0.556921, 8.294984], [1.0, 1.224863, 11.58736], [1.0, -1.347803, -2.406051], [1.0, 1.196604, 4.951851], [1.0, 0.275221, 9.543647], [1.0, 0.470575, 9.332488], [1.0, -1.889567, 9.542662], [1.0, -1.527893, 12.150579], [1.0, -1.185247, 11.309318], [1.0, -0.445678, 3.297303], [1.0, 1.042222, 6.105155], [1.0, -0.618787, 10.320986], [1.0, 1.152083, 0.548467], [1.0, 0.828534, 2.676045], [1.0, -1.237728, 10.549033], [1.0, -0.683565, -2.166125], [1.0, 0.229456, 5.921938], [1.0, -0.959885, 11.555336], [1.0, 0.492911, 10.993324], [1.0, 0.184992, 8.721488], [1.0, -0.355715, 10.325976], [1.0, -0.397822, 8.058397], [1.0, 0.824839, 13.730343], [1.0, 1.507278, 5.027866], [1.0, 0.099671, 6.835839], [1.0, -0.344008, 10.717485], [1.0, 1.785928, 7.718645], [1.0, -0.918801, 11.560217], [1.0, -0.364009, 4.7473], [1.0, -0.841722, 4.119083], [1.0, 0.490426, 1.960539], [1.0, -0.007194, 9.075792], [1.0, 0.356107, 12.447863], [1.0, 0.342578, 12.281162], [1.0, -0.810823, -1.466018], [1.0, 2.530777, 6.476801], [1.0, 1.296683, 11.607559], [1.0, 0.475487, 12.040035], [1.0, -0.783277, 11.009725], [1.0, 0.074798, 11.02365], [1.0, -1.337472, 0.468339], [1.0, -0.102781, 13.763651], [1.0, -0.147324, 2.874846], [1.0, 0.518389, 9.887035], [1.0, 1.015399, 7.571882], [1.0, -1.658086, -0.027255], [1.0, 1.319944, 2.171228], [1.0, 2.056216, 5.019981], [1.0, -0.851633, 4.375691], [1.0, -1.510047, 6.061992], [1.0, -1.076637, -3.181888], [1.0, 1.821096, 10.28399], [1.0, 3.01015, 8.401766], [1.0, -1.099458, 1.688274], [1.0, -0.834872, -1.733869], [1.0, -0.846637, 3.849075], [1.0, 1.400102, 12.628781], [1.0, 1.752842, 5.468166], [1.0, 0.078557, 0.059736], [1.0, 0.089392, -0.7153], [1.0, 1.825662, 12.693808], [1.0, 0.197445, 9.744638], [1.0, 0.126117, 0.922311], [1.0, -0.679797, 1.22053], [1.0, 0.677983, 2.556666], [1.0, 0.761349, 10.693862], [1.0, -2.168791, 0.143632], [1.0, 1.38861, 9.341997], [1.0, 0.317029, 14.739025]]
[0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0]

2. 数据可视化

数据已经提取出来，并且处理好了，下面就是可视化的过程了。解释在注释

def plotDataSet(weights):
    dataMat,labelMat = loadDataSet()  #加载数据
    dataArr = np.array(dataMat)   #转化为numpy格式，为画图做准备
    n = np.shape(dataMat)[0]      #样本个数
    #两种类型的x,y轴
    x1cord = []; y1cord = []
    x2cord = []; y2cord = []
    for i in range(n):
        if labelMat[i] == 1:     #假如为第一种情况
            x1cord.append(dataArr[i,1])    #因为0位置是1，偏置项，上文提到了
            y1cord.append(dataArr[i,2])
        else:
            x2cord.append(dataArr[i, 1])
            y2cord.append(dataArr[i, 2])
    #绘制图像
    fig = plt.figure()
    ax = fig.add_subplot(111)  #“111”表示“1×1网格，第一子图”，“234”表示“2×3网格，第四子图”。
    ax.scatter(x1cord,y1cord,s = 20,c = 'r', marker='s',alpha=0.5)   #绘制正样本
    ax.scatter(x2cord,y2cord,s = 20,c = 'b', alpha=0.5)   #绘制负样本
    plt.title('DataSet')  # 绘制title
    plt.xlabel('x');
    plt.ylabel('y')  # 绘制label
    plt.show()
#----------测试--------------------------------------------
if __name__ == '__main__':
    plotDataSet()

结果：

3. 梯度上升简单举例

在做梯度算法之前，我们先写一个简单的小例子,了解一下梯度上升（下降）的思路。
例如：去我们现在需要求下面函数的极值。

一般的做法是先求导数。

令导数为0，便求出了极值。这里的极值是极大值。
但是真实环境中的函数不会像上面这么简单，就算求出了函数的导数，也很难精确计算出函数的极值。此时我们就可以用迭代的方法来做。就像爬坡一样，一点一点逼近极值。这种寻找最佳拟合参数的方法，就是最优化算法。爬坡这个动作用数学公式表达即为：

其中α是学习率，后面是梯度。按梯度的方向逐步移动，为什么要按梯度呢？因为梯度方向是变化最大的方向呀。
下面我们来看一下代码具体实现的步骤。

def gradient():
    old = 0
    new = 0.1
    learn_rate = 0.01
    for i in range(1000):   #迭代次数
        new = old + learn_rate * (- 2 * old + 6)    #学习率乘梯度，向前移动的步长
        old = new
        print(new)
gradient()

结果：

……

自己运行一下，可以看出x的值依次迭代到靠近极大值3。
下面也是同理的。

4. 梯度上升的实际应用

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

下面的代码是梯度上升算法的具体实现。
这个例子中，我们是计算的真实类别与预测类别的差值，接下来就是按照该差值的方方向向调整回归系数，这里有个公式推导，先了解下这个推导过程再往下看。
代码如下。
先将sigmoid函数通过公式写出来方便调用。

def sigmoid(inX):
    result = 1.0 / (1 + np.exp(-inX))
    return result

然后开始写主体函数。更新参数的矢量化公式（上面公式推导最终结果，此例中只不过给矢量化了）如下。
θ = θ + αX^T(Y - g(Xθ)) 其中 g(Xθ) = sigmoid(WX) ，Y - g(Xθ)是真实类别与预测类别的差值

#梯度上升算法
def gradAscent(dataMatIn,classLables):
    dataMarix = np.mat(dataMatIn)    #创建矩阵，而不是数组了
    LabelMat = np.mat(classLables).transpose()
    m,n = np.shape(dataMarix)
    alpha = 0.001    #学习率      学习率需要多试，太小迭代次数不够，就到不了极值点。太大则步长太长，可能迈过了极值点。
    maxCycles = 500  #迭代次数
    weights = np.ones((n,1))   #设置n*1的权重矩阵
    weights_array = np.array([])
    for i in range(maxCycles):
        h = sigmoid(dataMarix * weights)   #公式
        error = LabelMat - h
        weights = weights + alpha * dataMarix.transpose() * error  ##梯度上升矢量化公式 θ = θ + αX^T(Y - g(Xθ)) 
        weights_array = np.append(weights_array, weights)
    # 因为plotBestFit()函数中有计算散点x,y坐标的部分，其中计算y的时候用到了weights，如果weights是矩阵的话，
    # weights[1]就是[[0.48007329]]（注意这里有中括号！），就不是一个数了，最终你会发现y的计算结果的len()只有1，
    # 而x的len()则是100
    weights_array = weights_array.reshape(maxCycles, n)
    return weights.getA(), weights_array #将矩阵转换为数组，返回权重数组，将自身返回成一个n维数组

该函数有两个参数。第⼀个参数是dataMatIn， 它是⼀个2维NumPy数组，每列分别代表每个不同的特征，每行则代表每个训练样本。我们现在采⽤的是 100个样本的简单数据集，它包含了两个特征X1和 X2，再加上第0维特征X0，所以dataMath⾥存放的将 是100×3的矩阵。我们获得输入数据并将它们转换成NumPy矩阵。NumPy对2维数组和矩阵都提供⼀些操作⽀持，如果混淆了数据类型和对应的操作，执⾏结果将与预期截然不同。第⼆个参数是类别标签，它是⼀个1×100的行向量。为了便于矩阵运算， 需要将该行向量转换为列向量，做法是将原向量转置，再将它赋值给labelMat。接下来的代码是得到矩阵大小，再设置⼀些梯度上升算法所需的参数。 变量alpha是向⽬标移动的步长，maxCycles是迭代次数。在for循环迭代完成后，将返回训练好的回归系数。需要强调的是，变量h不是⼀个数⽽是⼀个列向量，列向量的元素个数等于样本个数，这里是100。对应地，运算dataMatrix * weights代表的不是⼀次乘积计算，事实上该运算包含 了300次的乘积。
好，我们来看一下，经过梯度上升，参数W的值到底是多少了。
（注：θ是参数列向量，表示本例中的三个参数 [ W0(偏置项b) , W1 , W2] ^T

#测试
if __name__ == '__main__':
    dataMat, labelMat = loadDataSet()
    weights, weights_array = gradAscent(dataMat, labelMat)
    print(weights)

结果：

	[[ 4.12414349]
     [ 0.48007329]
     [-0.6168482 ]]

5. 通过图观察结果

直接引入上面图代码，添加些参数即可，如下。

def plotDataSet(weights):
    dataMat,labelMat = loadDataSet()  #加载数据
    dataArr = np.array(dataMat)   #转化为numpy格式
    n = np.shape(dataMat)[0]      #样本个数
    x1cord = []; y1cord = []
    x2cord = []; y2cord = []
    for i in range(n):
        if labelMat[i] == 1:     #假如为第一种情况
            x1cord.append(dataArr[i,1])
            y1cord.append(dataArr[i,2])
        else:
            x2cord.append(dataArr[i, 1])
            y2cord.append(dataArr[i, 2])
    #绘制图像
    fig = plt.figure()
    ax = fig.add_subplot(111)  #“111”表示“1×1网格，第一子图”，“234”表示“2×3网格，第四子图”。
    ax.scatter(x1cord,y1cord,s = 20,c = 'r', marker='s',alpha=0.5)   #绘制正样本
    ax.scatter(x2cord,y2cord,s = 20,c = 'b', alpha=0.5)   #绘制负样本
    #-----添加内容------绘制边界--------------------------
    #Z = w0 + W1X1 + W2X2    Z = 0是g(Z)是sigmoid的分解线，求的就是分界线，故令Z = 0    X2--->y
    x = np.arange(-3.0, 3.0, 0.1)
    y = (-weights[0] - weights[1] * x) / weights[2]
    #ax.plot(x,y)
    ax.plot(x, y)
    #----------------------------------------------------
    plt.title('DataSet')  # 绘制title
    plt.xlabel('x');
    plt.ylabel('y')  # 绘制label
    plt.show()

结果：

可以看出，结果还是可以接受的。
但是，尽管例⼦简单且数据集很小，这个方法却需要⼤量的计算（一次迭代就运行300次乘法)。所以我们仍需要改进它。

6. 随机梯度上升算法

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

下面是具体代码，相关解释可以见注释。

#随机梯度上升算法
def stocGradAscent1(dataMatrix,classLabels,number = 150):
    m, n = np.shape(dataMatrix)
    weights = np.ones(n)
    weights_array = np.array([])
    for j in range(number):   #迭代次数
        dataIndex = list(range(m))    #每次迭代在0----（m-1）找随机样本
        for i in range(m):
            alpha = 4/(1.0+j+i)+0.01     #学习率  alpha在每次迭代的时候都会调整，并且，虽然alpha会随着迭代次数不断减小，但永远不会减小到0，也好理解，越快靠近的时候步幅要放小。
            randIndex = int(random.uniform(0, len(dataIndex)))  #每次迭代在0----（m-1）找唯一的随机样本索引
            h = sigmoid(sum(dataMatrix[randIndex] * weights))  # 选择随机选取的一个样本，计算h
            error = classLabels[randIndex] - h  # 计算误差
            weights = weights + alpha * error * dataMatrix[randIndex]  # 更新回归系数
            weights_array = np.append(weights_array, weights, axis=0)
            #del (dataIndex[randIndex])  # 删除已经使用的样本
    weights_array = weights_array.reshape(number * m, n)   #改变维度
    return weights, weights_array

可以看到，随机梯度上升算法与梯度上升算法在代码上很相似，但也有⼀些区别：第⼀，后者的变量h和误差error都是向量，而前者则全是数值；第⼆，前者没有矩阵的转换过程，所有变量的数据类型都是 NumPy数组。 那回归的效果怎么样呢？来看一下吧。

#测试随机梯度下降
if __name__ == '__main__':
    dataMat, labelMat = loadDataSet()
    weights, weights_array= stocGradAscent1(np.array(dataMat), labelMat)
    plotDataSet(weights)

结果：

回归的也是不错的，但我们不能直观的看到每个回归方法的具体性能。因此，我们编写程序，绘制出回归系数和迭代次数的关系曲线

7. 两种算法比较

相关解释我会写在代码注释中。

def plotWeights(weights_array1,weights_array2):
    """
    Parameters:
        weights_array1 - 回归系数数组1
        weights_array2 - 回归系数数组2
    Returns:
    """
    #设置汉字格式
    font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=14)

    #将fig画布分隔成1行1列,不共享x轴和y轴,fig画布的大小为(13,8)
    #当nrow=3,nclos=2时,代表fig画布被分为六个区域,axs[0][0]表示第一行第一列
    fig, axs = plt.subplots(nrows=3, ncols=2,sharex=False, sharey=False, figsize=(20,10))
    #----梯度下降--------------------------------------------------
    x1 = np.arange(0, len(weights_array1), 1)     #x轴取值
    #绘制w0与迭代次数的关系
    axs[0][0].plot(x1,weights_array1[:,0])   #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 = np.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'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')
    plt.show()
#-----查看结果----------------------------
if __name__ == '__main__':
    dataMat, labelMat = loadDataSet()
    weights1,weights_array1 = gradAscent(dataMat, labelMat)
    weights2,weights_array2 = stocGradAscent1(np.array(dataMat), labelMat)
    plotWeights(weights_array1, weights_array2)

由于改进的随机梯度上升算法，随机选取样本点，所以每次的运行结果是不同的。但是大体趋势是一样的。我们改进的随机梯度上升算法收敛效果更好。为什么这么说呢？让我们分析一下。我们一共有100个样本点，改进的随机梯度上升算法迭代次数为150。而上图显示15000次迭代次数的原因是，使用一次样本就更新一下回归系数。因此，迭代150次，相当于更新回归系数150*100=15000次。简而言之，迭代150次，更新1.5万次回归参数。从上图右侧的改进随机梯度上升算法回归效果中可以看出，其实在更新2000次回归系数的时候，已经收敛了。相当于遍历整个数据集20次的时候，回归系数已收敛。训练已完成。
再让我们看看上图右侧的梯度上升算法回归效果，梯度上升算法每次更新回归系数都要遍历整个数据集。从图中可以看出，当迭代次数为300多次的时候，回归系数才收敛。
通过上面分析，随机梯度下降要好的多。并不是所有情况下都是随机的好，下个例子就可以说明问题。
到现在为止我们分析了回归系数的变化情况，但还没有达到本章的最终目标，即完成具体的分类任务。下⼀节将使⽤随机梯度上升算法来解决病马的生死预测问题。

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

本节将使用Logistic回归来预测患有疝病的马的存活问题。这里的数据包含368个样本和28个特征。疝病是描述马胃肠痛的术语。然而，这种病不⼀定源自马的胃肠问题， 其他问题也可能引发马疝病。该数据集中包含了医院检测马疝病的⼀指标，有的指标比较主观，有的指标难以测量，例如马的疼痛级别。
使用Logistic回归估计马疝病的死亡率流程

1. 收集数据：置顶资料下载处自行下载
2. 准备数据：用Python解析文本⽂件并填充缺失值。
3. 分析数据：可视化并观察数据。
4. 训练算法：使⽤优化算法，找到最佳的系数。
5. 测试算法：为了量化回归的效果，需要观察错误率。根据错误率决定是否回退到训练阶段，通过改变迭代的次数和步长等参数来得到更好的回归数。
6. 使用算法：实现⼀个简单的命令行程序来收集马的症状并输出预测结果。

1. 准备数据

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

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

预处理数据做两件事：
1. 如果测试集中一条数据的特征值已经确实，那么我们选择实数0来替换所有缺失值，因为本文使用Logistic回归。因此这样做不会影响回归系数的值。sigmoid(0)=0.5，即它对结果的预测不具有任何倾向性。
2. 如果测试集中一条数据的类别标签已经缺失，那么我们将该类别数据丢弃，因为类别标签与特征不同，很难确定采用某个合适的值来替换。
数据处理好了，开始代码吧！

2. 训练算法—随机梯度上升

首先，先导入相关包

"""
@Author:Yuuuuu、Tian
Sat Feb 22 09:12:09 2020
"""
import numpy as np
import random
import time

上面的例子是随机的算法好，那这个例子就从随机梯度算法开始吧！

def sigmoid(inX):
    return 1.0 / (1 + np.exp(-inX))
def stocGradAscent1(dataMatrix,classLabels,numIter = 150):
    m, n = np.shape(dataMatrix)
    weights = np.ones(n)
    for j in range(numIter):
        dataIndex = list(range(m))
        for i in range(m):
            alpha = 4 / (1.0 + i + j) + 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]
    return weights
    #分类器
def classfifyVector(inX,weights):
    prob = sigmoid((sum(inX * weights)))   #数组对应位置相乘，相加
    if prob > 0.5:
        return 1
    else:
        return 0

函数classifyVector()，它以回归系数和特征向量作为输入来计算对应的Sigmoid值。如果Sigmoid值⼤于0.5函数返回1，否则返回0。
相关算法写好了，开始我们的训练，测试吧！

def colicTest():
    #---数据处理---------------------------------
    frTrain = open('0702 DataSet/horseColicTraining.txt')
    frTest = open('0702 DataSet/horseColicTest.txt')   #打开测试集，训练集
    traingSet = []
    trainLabels = []
    for line in frTrain.readlines():
        currLine = line.strip().split('\t')
        lineArr = []
        for i in range(len(currLine) - 1):
            lineArr.append(float(currLine[i]))
        traingSet.append(lineArr)
        trainLabels.append(float(currLine[-1]))
    #--执行算法---------------------------------
    trainweights = stocGradAscent1(np.array(traingSet),trainLabels,numIter=500)
    #--测试集测试--------------------------------
    errorCount = 0
    testnum = 0.0
    #处理测试集数据
    for line in frTest.readlines():
        testnum += 1.0
        currLine = line.strip().split('\t')
        lineArr = []
        for i in range(len(currLine) - 1):
            lineArr.append(float(currLine[i]))
        if int(classfifyVector(np.array(lineArr),trainweights)) != int(currLine[-1]):
            errorCount += 1
    errorRate = (float(errorCount / testnum)) * 100
    print("测试集错误率为: %.2f%%" % errorRate)

colicTest()函数，是用于打开测试集和训练集，并对数据进行格式化处理的函数。该函数首先导入训练集，同前面⼀样，数据的最后⼀列仍然是类别标签。数据最初有三个类别标签，分别代表马的三 种情况：“仍存活”、“已经死亡”和“已经安乐死”。这里为了方便，将“已经死亡”和“已经安乐死”合并成“未能存活”这个标签。数据导⼊之后，便可以使⽤函数stocGradAscent1()来计算回归系数向量。这⾥可以自由设定迭代的次数，例如在训练集上使⽤500次迭代，实验结果表明这比默认迭代150次的效果更好。 在系数计算完成之后，导⼊测试集并计算分类错误率。整体看来，colicTest()具有完全独立的功能，多次运行得到的结果可能稍有不同，这是因为其中有随机的成分在里面。如果在stocGradAscent1()函数中回归系数已经完全收敛，那么结果才将是确定的。
来测试一下这个方法的性能。

# #进行测试
if __name__ == '__main__':
    s = time.clock()
    colicTest()
    e = time.clock()
    print('消耗时间：%f秒'%float(e - s))

结果：

测试集错误率为: 47.76%
消耗时间：1.458792秒

错误率还是蛮高的，而且耗时1.4s，并且每次运行的错误率也是不同的，错误率高的时候可能达到40%多。为啥这样？首先，因为数据集本身有30%的数据缺失，这个是不能避免的。另一个主要原因是，我们使用的是改进的随机梯度上升算法，因为数据集本身就很小，就几百的数据量。用改进的随机梯度上升算法显然不合适。让我们再试试梯度上升算法，看看它的效果如何？

3. 梯度上升算法

直接复制上个例子的即可，都一样的。

def gradAscent(dataMatIn,classLabels):
    dataMat = np.mat(dataMatIn)
    labelMat = np.mat(classLabels).transpose()
    m, n = np.shape(dataMat)
    alpha = 0.01
    numIter = 500
    weights = np.ones((n, 1))
    for i in range(numIter):
        h = sigmoid(dataMat * weights)
        error = labelMat - h
        weights = weights + alpha * dataMat.transpose() * error
    return weights.getA()  #将矩阵转换为数组，并返回

同上，再写一下训练，测试的函数

def colicTest():
    #---数据处理---------------------------------
    frTrain = open('0702 DataSet/horseColicTraining.txt')
    frTest = open('0702 DataSet/horseColicTest.txt')   #打开测试集，训练集
    traingSet = []
    trainLabels = []
    for line in frTrain.readlines():
        currLine = line.strip().split('\t')
        lineArr = []
        for i in range(len(currLine) - 1):
            lineArr.append(float(currLine[i]))
        traingSet.append(lineArr)
        trainLabels.append(float(currLine[-1]))
    #--执行算法---------------------------------
    trainweights = gradAscent(np.array(traingSet),trainLabels)
    #--测试集测试--------------------------------
    errorCount = 0
    testnum = 0.0
    #处理测试集数据
    for line in frTest.readlines():
        testnum += 1.0
        currLine = line.strip().split('\t')
        lineArr = []
        for i in range(len(currLine) - 1):
            lineArr.append(float(currLine[i]))
        if int(classfifyVector(np.array(lineArr),trainweights[:,0])) != int(currLine[-1]):
            errorCount += 1
    errorRate = (float(errorCount / testnum)) * 100
    print("测试集错误率为: %.2f%%" % errorRate)
#-----进行测试--------------------------------------------
if __name__ == '__main__':
    s = time.clock()
    colicTest()
    e = time.clock()
    print('消耗时间：%f秒'%float(e - s))

结果;

测试集错误率为: 28.36%
消耗时间：0.027854秒

通过结果可知，运行时间有了巨大的提升啊。
由这咱们就可以做个小结，当数据集较小时，我们使用梯度上升算法，当数据集较大时，我们使用改进的随机梯度上升算法。

五、 病马存活率 —（SKlearn）

每次最后，我们都要用skleran模块中的函数来实现一下我们的例子，你可能会问，有这么简单的方式，为什么还要写上面的东西呢？我们的目的是学习，要明白原理，这样面对问题解决起来才能游刃有余。
背景就不介绍了，和上面一样的。说说sklearn的Logistic回归函数吧！官方介绍
看一下官方的介绍，谷歌浏览器遇到英文会自带翻译，也不会看不懂。
只说下我们用到的参数。

• solver：优化算法选择参数，只有五个可选参数，即newton-cg, lbfgs, liblinear, sag, saga。默认为liblinear。liblinear适用于小数据集，而sag和saga适用于大数据集因为速度更快。
• max_iter：算法收敛最大迭代次数，int类型，默认为10。仅在正则化优化算法为newton-cg, sag和lbfgs才有用，算法收敛的最大迭代次数。
  相关只是了解了开始写代码吧！很简单就一段，解释还是看注释。

from sklearn.linear_model import LogisticRegression
def colicSklearn():
    frTrain = open('0702 DataSet/horseColicTraining.txt')
    frTest = open('0702 DataSet/horseColicTest.txt')  # 打开测试集，训练集
    trainingSet = []
    trainingLabels = []
    testSet = []
    testLabels = []
    for line in frTrain.readlines():
        currLine = line.strip().split('\t')
        lineArr = []
        for i in range(len(currLine) - 1):
            lineArr.append(float(currLine[i]))
        trainingSet.append(lineArr)
        trainingLabels.append(float(currLine[-1]))     #训练集数据，标签处理
    for line in frTest.readlines():
        currLine = line.strip().split('\t')
        lineArr = []
        for i in range(len(currLine) - 1):
            lineArr.append(float(currLine[i]))
        testSet.append(lineArr)
        testLabels.append(float(currLine[-1]))  #测试集数据，标签处理

    #Logistic回归
    classifier = LogisticRegression(solver='liblinear',max_iter=15).fit(trainingSet,trainingLabels)  #使用liblinear优化，最大迭代15次，用训练集拟合
    test_accurcy = classifier.score(testSet,testLabels) * 100    #用测试集测试
    print('正确率:%f%%' % test_accurcy)
#查看结果
if __name__ == '__main__':
    colicSklearn()

结果：

正确率:73.134328%

我们试试其它的优化算法，solver = liblinear在max_iter=10时有警告

说明还没有收敛，改为15就没有警告了。
更改solver='sag’max_iter=100时有警告，说明没有收敛，依次改500 1000 1500 2000。2000才收敛，可以看到，对于我们这样的小数据集，sag算法需要迭代上千次才收敛，而liblinear只需要不到20次。

六、 总结

Logistic回归的目的是寻找⼀个非线性函数Sigmoid的最佳拟合参数，求解过程可以由最优化算法来完成。 在最优化算法中，最常用的就是梯度上升算法，而梯度上升算法又可以简化为随机梯度上升算法。 随机梯度上升算法与梯度上升算法的效果相当，但占用更少的计算资源。此外，随机梯度上升是⼀个在线算法，它可以在新数据到来时就完成参数更新，⽽不需要重新读取整个数据集来进⾏批处理运算。 机器学习的⼀个重要问题就是如何处理缺失数据。这 个问题没有标准答案，取决于实际应用中的需求。现 有⼀些解决方案，每种方案都各有优缺点。 下⼀章将介绍与Logistic回归类似的另⼀种分类算法： ⽀持向量机，它被认为是目前最好的现成的算法之⼀，也是我认为的最难的算法了。好难啊。。。。。(●ˇ∀ˇ●)

• 本文结合各位大牛所思所想，不胜感激！
• 如有错误，请不吝指正！
• 喜欢的请点赞呀(★ ω ★)