机器学习 第五章 Logistic回归

Logistic回归是众多分类算法中的一员。
通常,Logistic回归用于二分类问题,例如预测明天是否会下雨。当然它也可以用于多分类问题,不过为了简单起见,本文暂先讨论二分类问题。首先,让我们来了解一下,什么是Logistic回归。

Logistic回归

假设现在有一些数据点,我们利用一条直线对这些点进行拟合(该线称为最佳拟合直线),这个拟合过程就称作为回归,如下图所示:

Logistic回归是分类方法,它利用的是Sigmoid函数阈值在[0,1]这个特性。

Logistic回归进行分类的主要思想是:根据现有数据对分类边界线建立回归公式,以此进行分类。其实,Logistic本质上是一个基于条件概率的判别模型(Discriminative Model)。

所以要想了解Logistic回归,我们必须先看一看Sigmoid函数 ,我们也可以称它为Logistic函数。
它的公式如下:
h θ ( x ) = g ( θ T x ) h_{\theta }(x)=g(\theta ^{T}x) hθ(x)=g(θTx) ,
z = [ θ 0 θ 1 . . . θ 0 ] [ x 0 x 1 . . . x n ] z=\begin{bmatrix} \theta _{0} &\theta _{1} & ... & \theta _{0} \end{bmatrix}\begin{bmatrix} x_{0}\\ x_{1} \\ ... \\ x_{n} \end{bmatrix} z=[θ0θ1...θ0] x0x1...xn ,
g ( z ) = 1 1 + e − z g(z)=\frac{1}{1+e^{-z}} g(z)=1+ez1 ,

整合成一个公式,就变成了如下公式:
h θ ( x ) = g ( θ T x ) = 1 1 + e − θ T x h_{\theta }(x)=g(\theta ^{T}x)=\frac{1}{1+e^{-\theta ^{T}x}} hθ(x)=g(θTx)=1+eθTx1 ,
下面这张图片,为我们展示了Sigmoid函数的样子。
在这里插入图片描述
z是一个矩阵,θ是参数列向量(要求解的),x是样本列向量(给定的数据集)。θ^T表示θ的转置。g(z)函数实现了任意实数到[0,1]的映射,这样我们的数据集([x0,x1,…,xn]),不管是大于1或者小于0,都可以映射到[0,1]区间进行分类。hθ(x)给出了输出为1的概率。比如当hθ(x)=0.7,那么说明有70%的概率输出为1。输出为0的概率是输出为1的补集,也就是30%。

如果我们有合适的参数列向量θ([θ0,θ1,…θn]^T),以及样本列向量x([x0,x1,…,xn]),那么我们对样本x分类就可以通过上述公式计算出一个概率,如果这个概率大于0.5,我们就可以说样本是正样本(分入1类),否则样本是负样本(分入0类)。
所以,Logistic回归也可以被看成是一种概率估计。

确定了分类器的函数形式之后,现在的问题变成了:最佳回归系数是多少?如何确定它们的大小?

梯度上升法

求得数据集的最佳参数的第一个最优化算法叫做梯度上升法
梯度上升法基于的思想是:要找到某函数的最大值,最好的方法是沿着该函数的梯度方向探寻。

从梯度向量几何意义上讲,就是函数变化增加最快的地方。
具体来说,对于函数f(x,y),在点(x0,y0),沿着梯度向量的方向就是 ( ∂ f ∂ x 0 , ∂ f ∂ y 0 ) T (\frac{ \partial f }{ \partial x_0 },\frac{ \partial f }{ \partial y_0 })^T (x0f,y0f)T的方向是f(x,y)增加最快的地方。或者说,沿着梯度向量的方向,更加容易找到函数的最大值。
反过来说,沿着梯度向量相反的方向,也就是 − ( ∂ f ∂ x 0 , ∂ f ∂ y 0 ) T -(\frac{ \partial f }{ \partial x_0 },\frac{ \partial f }{ \partial y_0 })^T (x0f,y0f)T的方向,梯度减少最快,也就是更加容易找到函数的最小值。

如果梯度记为 ▽ \bigtriangledown ,则函数f(x,y)的梯度由下式表示: ▽ f ( x , y ) = ( ∂ f ( x , y ) ∂ x ∂ f ( x , y ) ∂ y ) \bigtriangledown f(x,y)=\binom{\frac{\partial f(x,y)}{\partial x}}{\frac{\partial f(x,y)}{\partial y}} f(x,y)=(yf(x,y)xf(x,y))
这个梯度意味着要沿x的方向移动 ∂ f ( x , y ) ∂ x \frac{\partial f(x,y)}{\partial x} xf(x,y),沿y的方向移动 ∂ f ( x , y ) ∂ y \frac{\partial f(x,y)}{\partial y} yf(x,y)。其中,函数f(x,y)必须要在待计算的点上有定义并且可微。

梯度上升
梯度上升算法到达每个点后都会重新估计移动的方向。从P0开始,计算完该点的梯度,函数就根据梯度移动到下一点P1。在P1点,梯度再次被计算,并沿着新的梯度方向移动到P2。如此循环迭代,直到满足停止条件。迭代的过程中,梯度算子总能保证我们能选取到最佳的移动方向。从上图可以看到,梯度算子总是指向函数增长最快的方向。这里说的是移动方向,而未提到移动量的大小。该量来表示的话,梯度算法的迭代公式如下:
w : = w + α ▽ w f ( w ) w:=w+\alpha \bigtriangledown _{w}f(w) w:=w+αwf(w)
该公式将一直被迭代执行,直至达到某个停止条件为止,比如迭代次数达到某个指定值或算法达到某个可以允许的误差范围。

梯度上升算法代码实现

1、数据准备

这就是一个简单的数据集,没什么实际意义。让我们先从这个简单的数据集开始学习。先看下数据集有哪些数据:

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

这个数据有两维特征,因此可以将数据在一个二维平面上展示出来。我们可以将第一列数据(X1)看作x轴上的值,第二列数据(X2)看作y轴上的值。而最后一列数据即为分类标签。根据标签的不同,对这些点进行分类。

那么,先让我们编写代码,看下数据集的分布情况:

# -*- coding:UTF-8 -*-
import matplotlib.pyplot as plt
import numpy as np

"""
函数说明:加载数据

Parameters:
    无
Returns:
    dataMat - 数据列表
    labelMat - 标签列表
"""
def loadDataSet():
    dataMat = []                                                        #创建数据列表
    labelMat = []                                                        #创建标签列表
    fr = open('D://Study//python//textSet.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                                            #返回

"""
函数说明:绘制数据集

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

运行结果:

从上图可以看出数据的分布情况。假设Sigmoid函数的输入记为z,那么z=w0x0 + w1x1 + w2x2,即可将数据分割开。其中,x0为全是1的向量,x1为数据集的第一列数据,x2为数据集的第二列数据。另z=0,则0=w0 + w1x1 + w2x2。横坐标为x1,纵坐标为x2。这个方程未知的参数为w0,w1,w2,也就是我们需要求的回归系数(最优参数)。

2、训练算法

在编写代码之前,让我们回顾下梯度上升迭代公式:
在这里插入图片描述
将上述公式矢量化:
在这里插入图片描述
根据矢量化的公式,编写代码如下:
运行结果如图所示:

可以看出,我们已经求解出回归系数[w0,w1,w2]。

通过求解出的参数,我们就可以确定不同类别数据之间的分隔线,画出决策边界。

3、绘制决策边界

我们已经解出了一组回归系数,它确定了不同类别数据之间的分隔线。现在开始绘制这个分隔线,编写代码如下:

# -*- coding:UTF-8 -*-
import matplotlib.pyplot as plt
import numpy as np

"""
函数说明:加载数据

Parameters:
    无
Returns:
    dataMat - 数据列表
    labelMat - 标签列表
"""
def loadDataSet():
    dataMat = []                                                        #创建数据列表
    labelMat = []                                                        #创建标签列表
    fr = open('D://Study//python//textSet.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 + np.exp(-inX))

"""
函数说明:梯度上升算法

Parameters:
    dataMatIn - 数据集
    classLabels - 数据标签
Returns:
    weights.getA() - 求得的权重数组(最优参数)
"""
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()                                                #将矩阵转换为数组,返回权重数组

"""
函数说明:绘制数据集

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

if __name__ == '__main__':
    dataMat, labelMat = loadDataSet()           
    weights = gradAscent(dataMat, labelMat)
    plotBestFit(weights)

运行结果:

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

改进的随机梯度上升算法

梯度上升算法在每次更新回归系数(最优参数)时,都需要遍历整个数据集。可以看一下我们之前写的梯度上升算法:

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.01                                                        #移动步长,也就是学习速率,控制更新的幅度。
    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(),weights_array                                    #将矩阵转换为数组,返回权重数组

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

1、随机梯度上升算法

让我们直接看代码:

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[dataIndex[randIndex]]*weights))  #选择随机选取的一个样本,计算h
            error = classLabels[dataIndex[randIndex]] - h                           #计算误差
            weights = weights + alpha * error * dataMatrix[dataIndex[randIndex]]   #更新回归系数
            del(dataIndex[randIndex])                                         #删除已经使用的样本
    return weights                                                      #返回

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

编写代码如下,看下改进的随机梯度上升算法分类效果如何:

# -*- coding:UTF-8 -*-
from matplotlib.font_manager import FontProperties
import matplotlib.pyplot as plt
import numpy as np
import random

"""
函数说明:加载数据

Parameters:
    无
Returns:
    dataMat - 数据列表
    labelMat - 标签列表
"""
def loadDataSet():
    dataMat = []                                                        #创建数据列表
    labelMat = []                                                        #创建标签列表
    fr = open('D://Study//python//textSet.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 + np.exp(-inX))

"""
函数说明:绘制数据集

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

"""
函数说明:改进的随机梯度上升算法

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[dataIndex[randIndex]]*weights))         #选择随机选取的一个样本,计算h
            error = classLabels[dataIndex[randIndex]] - h                        #计算误差
            weights = weights + alpha * error * dataMatrix[dataIndex[randIndex]]   #更新回归系数
            del(dataIndex[randIndex])                                         #删除已经使用的样本
    return weights                                                            #返回

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

运行结果:

2、回归系数与迭代次数的关系

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

# -*- coding:UTF-8 -*-
from matplotlib.font_manager import FontProperties
import matplotlib.pyplot as plt
import numpy as np
import random


"""
函数说明:加载数据

Parameters:
    无
Returns:
    dataMat - 数据列表
    labelMat - 标签列表

"""
def loadDataSet():
    dataMat = []                                                        #创建数据列表
    labelMat = []                                                        #创建标签列表
    fr = open('D://Study//python//textSet.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 + np.exp(-inX))

"""
函数说明:梯度上升算法

Parameters:
    dataMatIn - 数据集
    classLabels - 数据标签
Returns:
    weights.getA() - 求得的权重数组(最优参数)
    weights_array - 每次更新的回归系数

"""
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.01                                                        #移动步长,也就是学习速率,控制更新的幅度。
    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 - 求得的回归系数数组(最优参数)
    weights_array - 每次更新的回归系数

"""
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[dataIndex[randIndex]]*weights))          #选择随机选取的一个样本,计算h
            error = classLabels[dataIndex[randIndex]] - h                           #计算误差
            weights = weights + alpha * error * dataMatrix[dataIndex[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'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'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)

运行结果如下:

由于改进的随机梯度上升算法,随机选取样本点,所以每次的运行结果是不同的。但是大体趋势是一样的。我们改进的随机梯度上升算法收敛效果更好。为什么这么说呢?让我们分析一下。我们一共有100个样本点,改进的随机梯度上升算法迭代次数为150。而上图显示15000次迭代次数的原因是,使用一次样本就更新一下回归系数。因此,迭代150次,相当于更新回归系数150*100=15000次。简而言之,迭代150次,更新1.5万次回归参数。从上图左侧的改进随机梯度上升算法回归效果中可以看出,其实在更新2000次回归系数的时候,已经收敛了。相当于遍历整个数据集20次的时候,回归系数已收敛。训练已完成。

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

改进的随机梯度上升算法,在遍历数据集的第20次开始收敛。而梯度上升算法,在遍历数据集的第300次才开始收敛。

总结

Logistic回归的一般过程:

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

Logistic回归的目的是寻找一个非线性函数Sigmoid的最佳拟合参数,求解过程可以由最优化算法完成。

  • 1
    点赞
  • 2
    收藏
    觉得还不错? 一键收藏
  • 1
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值