《机器学习实战》学习笔记（四）之Logistic（上）基础理论及算法推导、线性回归，梯度下降算法

转载请注明作者和出处：http://blog.csdn.net/john_bh/
运行平台： Windows
Python版本： Python3.6
IDE： Sublime text3

一、概述

Logistic回归是统计学习中的经典分类方法，也是众多回归算法中的一员。回归算法有：线性回归、Logistic回归、多项式回归、逐步回归、令回归、Lasso回归等。我们常用Logistic回归模型做预测。通常，Logistic回归用于二分类问题，例如预测明天是否会下雨。它也可以用于多分类问题。利用Logistic回归进行分类的主要思想是：根据现有数据对分类边界线建立回归公式，以此进行分类。

Logistic回归的一般过程：

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

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

在学习Logistic回归之前我们应该了解线性模型和线性回归，那么什么是线性模型和线性回归呢？

2.1 线性模型（linear model）

举个例子：如下图图1中左图所示，x轴是房屋的面积，Y轴表示的是房屋的价格，我们希望建立一个y=ax+b的模型，画一条现出来如图1中右图所示：

图1

但是，往往房屋的价格不只是由面积大小决定的，比如说：地段位置，房屋楼层，朝向，户型等等属性特征，那么我们可以通过向量的形式来表示这个模型，xx={x1,x2,...,xn},其中xi是xx在第i个属性上的取值，线性模型试图学的一个通过属性的线性组合来进行预测的函数，即：

f(xx)=w1x1+w2x2+...+wnxn+b

一般用向量形式写成

f(xx)=wwTxx+b

其中ww={w1,w2,...,wn}，ww和b学的之后模型就可以确定。
线性模型形式简单、易于建模，另外ww直观地表达了个属性在预测中的重要性，因此先行模型有很好的可解释性。例如西瓜书中提到的西瓜问题：f西瓜（x）=0.2⋅x色泽+0.5⋅x根蒂+0.3⋅x敲声+1，则意味着可通过综合考虑色泽、根蒂和敲声来判断西瓜好不好，其中根蒂最要紧，而敲声比色泽更重要。

2.2 线性回归

给定数据集D={(xx1,y1),(xx2,y2),...,(xxm,ym)}，其中xxi={xxi1,xxi2,...,xxin,},yi∈R，线性回归试图学得一个线性模型以尽可能准确地预测实值输出标记。

f(xi)=wxi+b，使得f(xi)≃yi

如何确定w和b呢？关键在于如何衡量f(x)和y之间的差别，均方误差是回归任务中最常用的性能度量，因此可以视图让均方误差最小化，即：

基于均方误差最小化进行模型求解的方法称为最小二乘法（least square method），在线性回归中，最小二乘法就是找到一条直线，是所有的样本到直线上的欧氏距离之和最小。

2.3 Logistic回归

我们想要的的函数应该是，能接受所有的输入，然后预测出类别，即输出标记y∈{0,1}，然而线性回归模型产生的预测值z=θθTxx+b是实值，于是我们需要将实值z转换为0/1值。最理想的就是单位阶跃函数(unit-step function),

但是单位阶跃函数是不连续的，这个瞬间跳跃过程有时候很难处理，所以我们希望找到一个在一定程度上近似单位阶跃，并且它是单调可微的函数，对数几率函数（Logistic function）正是这样的一个函数。

g(z)=11+e−z,其中z=θTx

对数几率函数是一种Sigmoid函数，将z带入得到

g（z）=11+e−(θθTxx+b)

g′(z)=(11+e−z)′=e−z(1+e−z)2=(11+e−z)(1−11+e−z)=g(z)(1−g(z))

整合成一个公式如下：

hθ(x)=g(θTx)=11+e−θTx

这就是Logistic回归模型。

用概率表示：
正例(y=1)：hθ(x)=P(y=1|x;θ)
正例(y=0)：1−hθ(x)=P(y=0|x;θ)
将两个公式合成一个：P(y|x;θ)=(hθ(x))y(1−hθ(x))1−y
最大似然概率：

L(θ)=p(y|X;θ)=∏i=1mp(y(i)|X(i);θ)=∏i=1m(hθ(x(i)))y(i)(1−hθ(x(i)))1−y(i)

取对数，求对数似然函数：

l(θ)=logL(θ)=∑i=1m[y(i)loghθ(x(i))+(1−y(i))log(1−hθ(x(i)))]

其中，m为样本的总数，y(i)表示第i个样本的类别，x(i)表示第i个样本，需要注意的是θ是多维向量，x(i)也是多维向量。

接下来，我们需要通过所有的数据来训练并且学习出来一组θ的值，使得l(θ)的值最小化。满足l(θ)的最小的θ值即是我们需要求解的模型。

怎么求解使l(θ)最小的θ值呢？这里使用使用梯度下降算法。如果面对的问题是求解使l(θ)最大的θ值，那么我们就需要使用梯度上升算法。

2.4 梯度下降算法

梯度下降算法基于的思想是：要找到某函数的最小值，最好的方法就是沿着该函数的梯度方向去寻找，如果梯度记为∇，则函数l(θ)的梯度可以表示为：

∇l(θ)=∂∂θjl(θ)

对l(θ)求θ的偏导，得出：

因此，梯度下降算法为：

θj=θj−α∂∂θjl(θ)，其中α是学习率、步长

梯度下降参数迭代(随机梯度下降)的更新法则：

θj:=θj−α(y(i)−hθ(x(i)))x(i)j

该公式将被一直迭代执行，知道达到某个停止条件为止，比如说迭代次数达到某个指定值或者达到某个可以允许的误差范围。

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

我们知道梯度下降算法，其实梯度上升算法和梯度下降算法是类似的，实施工事中的减法变成了加法，因此公式可以写成：

θj=θj+α∂∂θjl(θ)，其中α是学习率、步长

梯度下降算法是求的函数的最小值，而梯度上升算法是求的函数的最大值。

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

梯度上升的伪代码：

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

新建文件Logistic.py，具体代码实现如下：

# -*- coding:utf-8 -*-import matplotlib.pyplot as pltimport 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     #返回"""函数说明:sigmoid函数Parameters:    inX - 数据Returns:    返回sigmoid函数"""def sigmoid(inX):    return 1.0/(1+exp(-inX))"""函数说明:梯度上升算法Parameters:    dataMatIn - 数据集    classLabels - 数据标签Returns:    weights.getA() - 求得的权重数组(最优参数)"""def gradAscent(dataMatIn,classLables):    dataMatrix=np.mat(dataMatIn)    #转换成numpy的mat    lableMat=np.mat(classLables).transpose()   #转换成numpy的mat,并进行转置    m,n=np.shape(dataMatrix)    #返回dataMatrix的行数和列数：m为行数,n为列数。    alpha=0.001    #移动步长,或者叫做学习速率,控制更新的幅度。    maxCycles=500    #最大迭代次数    weights=np.ones((n,1))    #参数θ默认为1    for k in range(maxCycles):        h=sigmoid(dataMatrix * weights)     #梯度上升矢量化公式,通过sigmoid函数得到h_θ(x)        error = (lableMat - h)   #计算真实值和与测试值的差值，按照差值方向调整回归系数        weights = weights + alpha *dataMatrix.transpose() * error    #参数迭代生成新的参数θ    return weights    #返回参数向量"""函数说明:绘制数据集Parameters:    无Returns:    无"""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()      #显示if __name__ == '__main__':    #plotDataSet()    dataMat, labelMat = loadDataSet()    print(gradAscent(dataMat, labelMat))

程序运行结果如下图所示：

图3 最佳回归参数输出结果

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

我们已经解出了一组回归系数，它确定了不同类别数据之间的分隔线。现在开始绘制这个分隔线，在文件Logistic.py文件中添加函数plotBestFit(weights)，编写代码如下：

"""函数说明:绘制决策边界Parameters:    weights - 参数数组    dataMat - 数据集    labelMat - 数据标签Returns:    无"""def plotBestFit(weights,dataMat,labelMat):    dataArr=np.array(dataMat)   #转换成numpy的array数组    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])    #1为正样本        else:            xcord2.append(dataArr[i,1]);ycord2.append(dataArr[i,2])    #0为负样本    fig=plt.figure()    ax=fig.add_subplot(111)    ax.scatter(xcord1,ycord1,s=30,c='red',marker='s',alpha=.5)    #绘制正样本    ax.scatter(xcord2,ycord2,s=30,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()

得到的决策边界线如下图所示：

图4 Logistic回归最佳拟合直线

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

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

梯度上升算法在每次更新回归系数(最优参数)时，都需要遍历整个数据集。假设，我们使用的数据集一共有100个样本。那么，dataMatrix就是一个100*3的矩阵。每次计算h的时候，都要计算dataMatrix*weights这个矩阵乘法运算，要进行100*3次乘法运算和100*2次加法运算。同理，更新回归系数(最优参数)weights时，也需要用到整个数据集，要进行矩阵乘法运算。总而言之，该方法处理100个左右的数据集时尚可，但如果有数十亿样本和成千上万的特征，那么该方法的计算复杂度就太高了。因此，需要对算法进行改进，我们每次更新回归系数(最优参数)的时候，能不能不用所有样本呢？一次只用一个样本点去更新回归系数(最优参数)？这样就可以有效减少计算量了，这种方法就叫做随机梯度上升算法。

由于可以在新样本到来时对分类器进行增量更新，因而随机梯度上升算法是一个在线学习算法，与“在线学习算法”相对应，一次处理所有数据被称为“批处理”，也叫批量梯度下降/上升算法。

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

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

具体代码实现，在Logistic中添加函数stocGradAscent0()，

"""函数说明:随机梯度上升算法Parameters:    dataMatrix - 数据集    classLabels - 数据标签Returns:    weights - 求得的权重数组(最优参数)"""def stocGradAscent0(dataMatrix,classLables):    m,n=np.shape(dataMatrix)    #返回dataMatrix的行数和列数：m为行数,n为列数。    alpha=0.01    #移动步长,或者叫做学习速率,控制更新的幅度。    weights=np.ones(n)    #参数θ默认为1    for i in range(m):        h=sigmoid(np.sum(dataMatrix[i] * weights))     #梯度上升矢量化公式,通过sigmoid函数得到h_θ(x)        error = (classLables[i] - h)    #计算真实值和与测试值的差值，按照差值方向调整回归系数        weights = weights + alpha *dataMatrix[i] * error    #参数迭代生成新的参数θ    return weights    #返回参数向量

可以看出，随机梯度上升算法和梯度上升算法在代码上很相似，但也有不同：第一，后者的变量h和误差error都是向量，而前者则全是数值；第二，前者没有矩阵的转换过程，所有变量的数据类型都是NumPy数组。
运行结果如下图所示：

图5 随机梯度上升算法执行结果

由图可得知，拟合出来的直线效果还不错，但是并不像图4那样完美，这里的分类器错分了三分之一的样本。
判断一个优化算法的优劣可靠地方法是看它是否收敛，也就是说参数是否达到了稳定值，是否还在不断地变化。对此我们将函数stocGradAscent0()做一些修改，具体代码如下：

"""函数说明:改进的随机梯度上升算法Parameters:    dataMatrix - 数据集    classLabels - 数据标签    numIter - 迭代次数Returns:    weights - 求得的权重数组(最优参数)"""def stocGradAscent1(dataMatrix,classLables,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(np.sum(dataMatrix[randIndex] * weights))     # #选择随机选取的一个样本，计算h            error = (classLables[randIndex] - h)    #计算误差            weights = weights + alpha * error * dataMatrix[randIndex]      #更新回归系数            del(dataIndex[randIndex])    #删除已经使用的样本    return weights    #返回参数向量

该算法第一个改进之处在于，alpha在每次迭代的时候都会调整，并且，虽然alpha会随着迭代次数不断减小，但永远不会减小到0，因为这里还存在一个常数项。必须这样做的原因是为了保证在多次迭代之后新数据仍然具有一定的影响。如果需要处理的问题是动态变化的，那么可以适当加大上述常数项，来确保新的值获得更大的回归系数。另一点值得注意的是，在降低alpha的函数中，alpha每次减少1/(j+i)，其中j是迭代次数，i是样本点的下标。

第二个改进的地方，这里通过随机选取样本来更新回归系数，这种方法减少周期性的波动，每次迭代不使用已经用过的样本点。这样的方法，就有效地减少了计算量，并保证了回归效果。

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

图6 改进的随机梯度上升算法执行结果

3.4 讨论回归系数与迭代次数的关系

可以看到分类效果也是不错的。不过，从这个分类结果中，我们不容易看出迭代次数和回归系数的关系，也就不能直观的看到每个回归方法的收敛情况。因此，我们编写程序，绘制出回归系数和迭代次数的关系曲线如下图：

图6 改进的随机梯度上升算法执行结果

我们看到上图左侧的是梯度上升算法回归效果，梯度上升算法每次更新回归系数都要遍历整个数据集。从图中可以看出，总共迭代500次，当迭代次数为300多次的时候，回归系数才收敛。但是在大的波动停止后，还有一些小的周期性波动，产生这种现象的原因：存在一些不能正确分类的样本点（数据集并非线性可分），在每次迭代时，会引起系数的剧烈改变。

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

四、 小结

本文主要学习了线性回归，Logistic回归的基本原理和算法推导，还学习了梯度下降算法并且对梯度上升算法做了改进，那么下一篇文章将这次学的理论进行实践开发，运用Logistic回归实预测病马死亡率。