开发: C++知识库 Java知识库 JavaScript Python PHP知识库人工智能区块链大数据移动开发嵌入式开发工具数据结构与算法开发测试游戏开发网络协议系统运维
教程: HTML教程 CSS教程 JavaScript教程 Go语言教程 JQuery教程 VUE教程 VUE3教程 Bootstrap教程 SQL数据库教程 C语言教程 C++教程 Java教程 Python教程 Python3教程 C#教程
数码: 电脑笔记本显卡显示器固态硬盘硬盘耳机手机 iphone vivo oppo 小米华为单反装机图拉丁

-> 人工智能 -> 《机器学习实战》—— Logistic 回归 -> 正文阅读

[人工智能]《机器学习实战》—— Logistic 回归

一、Logistic 回归

假设现在有一些数据点，我们用一条直线对这些点进行拟合（该线称为最佳拟合直线），这个拟合过程就称作回归。
在这里插入图片描述

利用Logistic回归进行分类的主要思想是：根据现有数据对分类边界线建立回归公式，以此进行分类。
这里的“回归”一词源于最佳拟合，表示要找到最佳拟合参数集，其背后的数学分析将在下一部分介绍。训练分类器时的做法就是寻找最佳拟合参数，使用的是最优化算法。

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

其实，Logistic本质上是一个基于条件概率的判别模型(Discriminative Model)。

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

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

我们想要的函数应该是，能接受所有的输入然后预测出类别。例如，在两个类的情况下，上述函数输出0或1。该函数称为海维塞德阶跃函数（Heaviside step function），或者直接称为单位阶跃函数。然而，海维塞德阶跃函数的问题在于：函数不连续，该函数在跳跃点上从0瞬间跳跃到1，这个瞬间跳跃过程有时很难处理。

$\begin{cases} \quad 0,\quad z < 0 \\ 0.5,\quad \ z = 0 \\ \quad 1,\quad z > 0 \\ \end{cases}$
在这里插入图片描述

幸好，另一个函数也有类似的性质，且数学上更易处理，这就是Sigmoid函数。Sigmoid函数具体的计算公式如下：
$\sigma(z) = \frac{1}{1 + e^{-z} }$

在这里插入图片描述
如图中上图给出了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$ ，由下面公式得出：
$w_0x_0 + w_1x_1 + w_2x_2 + \ldots + w_nx_n$

如果采用向量的写法，上述公式可以写成 $z = w^Tx$ ，它表示将这两个数值向量对应元素相乘然后全部加起来即得到z值。其中的向量x是分类器的输入数据，向量w也就是我们要找到的最佳参数（系数），从而使得分类器尽可能地精确。为了寻找该最佳参数，从而使用梯度上升算法进行求解。

梯度上升法基于的思想是：要找到某函数的最大值，最好的方法是沿着该函数的梯度方向探寻。如果梯度记为 $\nabla$ ，则函数 $f (x, y)$ 的梯度由下式表示：
$\nabla f(x,y) = \begin{Bmatrix} \LARGE \frac{ \partial f(x,y)}{\partial x } \\ \\ \LARGE \frac{ \partial f(x,y)}{\partial y } \end{Bmatrix}$

这个梯度意味着要沿 $x$ 的方向移动 $\large\frac{ \partial f(x,y)}{\partial x }$ ，沿 $y$ 的方向移动 $\large\frac{ \partial f(x,y)}{\partial x }$ 。其中，函数 $f (x, y)$ 必须要在待计算的点上有定义并且可微。

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

上图中的梯度上升算法沿梯度方向移动了一步。可以看到，梯度算子总是指向函数值增长最快的方向。这里所说的是移动方向，而未提到移动量的大小。该量值称为步长，记做 $\alpha$ 。用向量来表示的话，梯度上升算法的迭代公式如下：
$w+\alpha\nabla_wf(w)$
该公式将一直被迭代执行，直至达到某个停止条件为止，比如迭代次数达到某个指定值或算法达到某个可以允许的误差范围。

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

4.1 查看数据的分布情况

在这里插入图片描述

这个数据有两维特征，因此可以在一个二维平面上展示。首先将第一列数据 $X_1$ 看作x轴上的值，然后第二列数据 $X_2$ 看作y轴上的值，最后把最后一列数据即为分类标签。根据标签的不同，对这些点进行分类。

import matplotlib.pyplot as plt
import numpy as np

"""
Parameters:
    无
Returns:
    dataMat - 数据列表
    labelMat - 标签列表
"""
# 函数说明:加载数据
def loadDataSet():
    dataMat = []                                                        #创建数据列表
    labelMat = []                                                       #创建标签列表
    fr = open('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 plotDataSet():
    dataMat, labelMat = loadDataSet()                                   #加载数据集
    dataArr = np.array(dataMat)                                         #转换成numpy的array数组
    n = np.shape(dataMat)[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])    #1为正样本
        else:
            xcord2.append(dataArr[i,1]); ycord2.append(dataArr[i,2])    #0为负样本
    fig = plt.figure()
    ax = fig.add_subplot(111)                                           #添加subplot
    ax.scatter(xcord1, ycord1, s = 20, c = 'red', marker = 's',alpha=.5)#绘制正样本
    ax.scatter(xcord2, ycord2, s = 20, c = 'green',alpha=.5)            #绘制负样本
    plt.title('DataSet')                                                #绘制title
    plt.xlabel('x'); plt.ylabel('y')                                    #绘制label
    plt.show()                                                          #显示

在这里插入图片描述
从上图可以看出我们采用的数据的分布情况。

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

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

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

每个回归系数初始化为1
重复R次：
	计算整个数据集的梯度
	使用alpha × gradient更新回归系数的向量
	返回回归系数

'''
Parameters:
    inX - 数据
Returns:
    sigmoid函数
'''
# 函数说明:sigmoid函数
def sigmoid(inX):
    return 1.0 / (1 + np.exp(-inX))

'''
Parameters:
    dataMatIn - 数据集
    classLabels - 数据标签
Returns:
'''
# 函数说明:梯度上升算法
def gradAscent(dataMatIn, classLabels):
    dataMatrix = np.mat(dataMatIn)                                  #转换成numpy的mat
    labelMat = np.mat(classLabels).transpose()                      #转换成numpy的mat,并进行转置
    m, n = np.shape(dataMatrix)                                     #返回dataMatrix的大小。m为行数,n为列数。
    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.getA()                                           #将矩阵转换为数组，返回权重数组

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

该结果即为使用梯度上升算法找到的一组回归系数。

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

已经解出了一组回归系数。现在开始绘制这个分隔线:

'''
Parameters:
    weights - 权重参数数组
Returns:
    无
'''
# 函数说明:绘制数据集
def plotBestFit(weights):
    dataMat, labelMat = loadDataSet()                                   #加载数据集
    dataArr = np.array(dataMat)                                         #转换成numpy的array数组
    n = np.shape(dataMat)[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])    #1为正样本
        else:
            xcord2.append(dataArr[i,1]); ycord2.append(dataArr[i,2])    #0为负样本
    fig = plt.figure()
    ax = fig.add_subplot(111)                                           #添加subplot
    ax.scatter(xcord1, ycord1, s = 20, c = 'red', marker = 's',alpha=.5)#绘制正样本
    ax.scatter(xcord2, ycord2, s = 20, c = 'green',alpha=.5)            #绘制负样本
    x = np.arange(-3.0, 3.0, 0.1)
    y = (-weights[0] - weights[1] * x) / weights[2]
    ax.plot(x, y)
    plt.title('BestFit')                                                #绘制title
    plt.xlabel('X1'); plt.ylabel('X2')                                  #绘制label
    plt.show()

在这里插入图片描述
这个分类结果相当不错，从上图可以看出，只分错了几个点而已。但是，尽管例子简单切数据集很小，但是这个方法却需要大量的计算(300次乘法)。因此还是需要对算法稍作改进，从而减少计算量，使其可以应用于大数据集上。

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

梯度上升算法在每次更新回归系数时都需要遍历整个数据集，如果数据集很大，具有成千上万的特征，那么该方法的计算复杂度就太高了。

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

随机梯度上升算法可以写成如下的伪代码：

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

改进算法还增加了一个迭代次数作为第3个参数。如果该参数没有给定的话，算法将默认迭代150次。如果给定，那么算法将按照新的参数值进行迭代。代码如下：

'''
Parameters:
    dataMatrix - 数据数组
    classLabels - 数据标签
Returns:
    weights - 求得的回归系数数组(最优参数)
'''
def stocGradAscent0(dataMatrix, classLabels):
    m,n = np.shape(dataMatrix)                                          #返回dataMatrix的大小。m为行数,n为列数。
    alpha = 0.01
    weights = np.ones(n)                                                #参数初始化
    for i in range(m):
        h = sigmoid(sum(dataMatrix[i]*weights))                         #选择随机选取的一个样本，计算h
        error = classLabels[i] - h                                      #计算误差
        weights = weights + alpha*error*dataMatrix[i]                   #更新回归系数
    return weights

在这里插入图片描述
运行结果显示，红色将近一半分类错误，即使用随机梯度上升算法进行分类的错误率明显高于梯度上升算法，这是由于一些样本点不能被正确的分类（数据集并非线性可分），在每次迭代时会引发系数的剧烈改变。因此我们希望该算法能避免来回波动，从而收敛到某个值，且收敛速度也需要加快。

4.5 训练算法：改进随机梯度上升

将系数alpha在每次迭代时都进行调整，alpha会随着迭代次数不断减小，但永远不会减小到0。并通过随机选取样本来更新回归系数，该方法每次随机从列表中选出一个值，然后从列表中删除该值（再进行下一次迭代）

'''
Parameters:
    dataMatrix - 数据数组
    classLabels - 数据标签
    numIter - 迭代次数
Returns:
    weights - 求得的回归系数数组(最优参数)
'''
# 函数说明:改进的随机梯度上升算法
def stocGradAscent1(dataMatrix, classLabels, numIter=150):
    m,n = np.shape(dataMatrix)                                       #返回dataMatrix的大小。m为行数,n为列数。
    weights = np.ones(n)                                             #参数初始化
    for j in range(numIter):
        dataIndex = list(range(m))
        for i in range(m):
            alpha = 4/(1.0+j+i)+0.01                                 #降低alpha的大小，每次减小1/(j+i)。
            randIndex = int(random.uniform(0,len(dataIndex)))        #随机选取样本
            h = sigmoid(sum(dataMatrix[randIndex]*weights))          #选择随机选取的一个样本，计算h
            error = classLabels[randIndex] - h                       #计算误差
            weights = weights + alpha * error * dataMatrix[randIndex]#更新回归系数
            del(dataIndex[randIndex])                                #删除已经使用的样本
    return weights                                                   #返回

在这里插入图片描述
运行后的分类效果是很好的，只有两个点分类错误。

不过，从这个分类结果中，不好看出迭代次数和回归系数的关系，也就不能直观的看到每个回归方法的收敛情况。

4.6 回归系数与迭代次数的关系

编写程序，绘制出回归系数和迭代次数的关系曲线：

'''
Parameters:
    dataMatIn - 数据集
    classLabels - 数据标签
Returns:
'''
# 函数说明:梯度上升算法
def gradAscent(dataMatIn, classLabels):
    dataMatrix = np.mat(dataMatIn)                                  #转换成numpy的mat
    labelMat = np.mat(classLabels).transpose()                      #转换成numpy的mat,并进行转置
    m, n = np.shape(dataMatrix)                                     #返回dataMatrix的大小。m为行数,n为列数。
    alpha = 0.001                                                   #移动步长,也就是学习速率,控制更新的幅度。
    maxCycles = 500                                                 #最大迭代次数
    weights = np.ones((n,1))
    weights_array = np.array([])
    for k in range(maxCycles):
        h = sigmoid(dataMatrix * weights)                           #梯度上升矢量化公式
        error = labelMat - h
        weights = weights + alpha * dataMatrix.transpose() * error
        weights_array = np.append(weights_array, weights)
    weights_array = weights_array.reshape(maxCycles, n)
    return weights.getA(),weights_array                             #将矩阵转换为数组，返回权重数组

'''
Parameters:
    dataMatrix - 数据数组
    classLabels - 数据标签
    numIter - 迭代次数
Returns:
    weights - 求得的回归系数数组(最优参数)
'''
# 函数说明:改进的随机梯度上升算法
def stocGradAscent1(dataMatrix, classLabels, numIter=150):
    m,n = np.shape(dataMatrix)                                       #返回dataMatrix的大小。m为行数,n为列数。
    weights = np.ones(n)                                             #参数初始化
    weights_array = np.array([])                                     #存储每次更新的回归系数
    for j in range(numIter):
        dataIndex = list(range(m))
        for i in range(m):
            alpha = 4/(1.0+j+i)+0.01                                 #降低alpha的大小，每次减小1/(j+i)。
            randIndex = int(random.uniform(0,len(dataIndex)))        #随机选取样本
            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(numIter * m, n)               # 改变维度
    return weights,weights_array                                        # 返回

'''
Parameters:
    weights_array1 - 回归系数数组1
    weights_array2 - 回归系数数组2
Returns:
    无
'''
# 函数说明:绘制回归系数与迭代次数的关系
def plotWeights(weights_array1,weights_array2):
    #设置汉字格式
    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)
    #绘制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'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')


    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'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()       

if __name__ == '__main__':
    dataMat, labelMat = loadDataSet()
    weights1,weights_array1 = stocGradAscent1(np.array(dataMat), labelMat)

    weights2,weights_array2 = gradAscent(dataMat, labelMat)
    plotWeights(weights_array1, weights_array2)

在这里插入图片描述

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

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

五、根据建立的气候模型来预测是否会崩溃

5.1 数据集

本数据集是在UCI Machine Learning Repository网站下载
在这里插入图片描述
一本有540条数据，每条数据有17个特征。现将数据集中80%数据作为训练集，剩下的20%数据作为测试集。

5.2 预测

'''
Parameters:
    inX - 特征向量
    weights - 回归系数
Returns:
    分类结果
'''
# 函数说明:分类函数
def classifyVector(inX, weights):
    prob = sigmoid(sum(inX*weights))
    if prob > 0.5: return 1.0
    else: return 0.0


"""
函数说明：使用Logistic分类器进行预测
Parameters:
    无
Returns:
    无
"""
# 函数说明:使用Python写的Logistic分类器做预测
def colicTest():
    frTrain = open('pop_failures_training.txt')                                          #打开训练集
    frTest = open('pop_failures_test.txt')                                               #打开测试集
    trainingSet = []; trainingLabels = []
    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]))
    trainingSet = np.array(trainingSet)
    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(len(currLine)-1):
            lineArr.append(float(currLine[i]))
        if int(classifyVector(np.array(lineArr), trainWeights))!= int(currLine[-1]):
            errorCount += 1
    errorRate = (float(errorCount)/numTestVec) * 100                                 #错误率计算
    print("测试集错误率为: %.2f%%" % errorRate)
    return errorRate
"""
函数说明：计算迭代numtests次后的错误率
"""
def multiTest():
    numtests = 10;errorsum = 0.0
    for k in range(numtests):
        errorsum += colicTest()
    print("%d次分类测试的平均错误率为：%.2f%%" % (numtests, errorsum/float(numtests)))

>>>
测试集错误率为: 6.48%
测试集错误率为: 8.33%
测试集错误率为: 7.41%
测试集错误率为: 8.33%
测试集错误率为: 8.33%
测试集错误率为: 6.48%
测试集错误率为: 6.48%
测试集错误率为: 6.48%
测试集错误率为: 9.26%
测试集错误率为: 8.33%
10次分类测试的平均错误率为：7.59%