逻辑回归从理论到Python实现.整理

1. 逻辑回归的模型函数

前面我们讲了线性回归模型，将线性模型用于回归问题中。这篇我们讲一下线性模型用于分类任务。在二分类问题中，对于线性回归所产生的预测值

我们需将这个预测z转化为0/1值，最理想的是“单位阶跃函数”，即若预测z大于零就判为1，若预测z小于零则判为反例，预测值为临界值则可以任意判别。但是由于单位阶跃函数不连续，我们希望可以找到一个连续的，同时在一定程度上近似“单位阶跃函数“的函数，那么逻辑斯特函数(logistic function)则是一个很不错的替代函数：

如下图所示，这是一种"Sigmoid"函数，它z值转化为一个接近0或1的y值。

我们将

代入上式，最终得到逻辑回归的模型函数，如下

作者个人理解：“单位阶跃函数”是将预测值直接转换为分类，而“逻辑斯特函数”是将预测值转换为0-1之间的一个数y，这个值正好可以当成某一类别的概率值，从而实现分类。

2. 逻辑回归的目标函数

逻辑回归的模型函数已经知道，求逻辑回归的模型函数，其实就是确定该模型函数theta这个参数的值。在讲如何求这个参数之前，需要求得逻辑回归的目标函数，而在求解目标函数的过程中，需先讲解一下以下知识点：

一个事件的几率（odds）：指该事件发生与不发生的概率比值，若事件发生概率p，那么事件发生的几率就是

那么该事件的对数几率（log odds,亦称logit）就是：

将式子

转化为如下

也就是说，输出y的对数几率是由输入x的线性函数表示的模型，这就是逻辑回归模型。当

的值越接近正无穷，y的值是越接近1的，那么我们将y视为样本x为正例，即y=1的概率。那么上述公式可转化为如下

从而得

进一步可转化为

在统计学中，常常使用极大似然估计法来求解，即找到一组参数，使得在这组参数下，我们的数据的似然度（概率）最大，似然函数如下：

                   
取对数似然函数得

我们要最大化上式，也就等价于最小化

上式即可当做逻辑回归的目标函数

在逻辑回归模型中，我们最大化似然函数和最小化对数似然损失函数实际上是等价的。

3. 求解目标函数

上一块知识讲解中已经得到了目标函数，这一块知识讲解中，主要是讲解如何得到目标函数的最小值。自然而来的还是想到了梯度下降算法。我们求上述目标函数的梯度得

其中

那么对于

的一次迭代如下:

从而推出对于theta的一次迭代如下（并进行了向量化）

其中

上面将求和公式转化为了向量的形式，具体推导的过程已省略，可自行推导，也可以参考文末的梯度下降算法的向量化推导的参考链接。

4. 代码的实现

上面讲了逻辑回归的理论，下面我们通过代码来实现逻辑回归，例子主要是机器学习实战书上的例子。首先来看一下部分数据集（数据集的地址见文末给的链接地址），如下

-0.017612    14.053064   0
-1.395634    4.662541    1
-0.752157    6.538620    0
-1.322371    7.152853    0
...
-2.168791    0.143632    1
1.388610    9.341997    0
0.317029    14.739025   0

我们需要从文件中读取数据，并返回读取之后的列表结果，实现的代码如下

'''
从文件中读取数据
'''
def loadDataSet(filePath):
    dataList = []
    labelList = []
    f = open(filePath)
    for line in f.readlines():
        lineArr = line.strip().split()
        dataList.append([1.0, float(lineArr[0]), float(lineArr[1])])
        labelList.append(int(lineArr[2]))
    return dataList, labelList

之后我们采用梯度下降算法进行求解，在迭代过程中使用向量化，实现的代码如下：

'''
sigmoid函数
'''
def sigmoid(inX):
    return 1.0 / (1 + np.exp(-inX))


'''
梯度下降函数
'''
def gradAscent(dataList, labelList):
    dataMat = np.mat(dataList)  # x数据转化为矩阵
    labelMat = np.mat(labelList).T  # y数据转换为矩阵并转置为列向量
    m, n = np.shape(dataMat)  # 返回矩阵的大小
    alpha = 0.001  # 步长
    maxCycles = 500  # 迭代次数
    weights = np.ones((n, 1))  # 权重列向量
    # 梯度下降算法
    for k in range(maxCycles):
        h = sigmoid(dataMat * weights)
        error = h - labelMat
        weights = weights - alpha * dataMat.T * error
    return weights

通过函数gradAscent可以得到一组回归系数，它确定了不同类别数据之间的分隔线，下面我们将这个分割线画出来，代码如下：

'''
画分类的分割线
'''
def plotBestFit(dataList, labelList, weights):
    dataArr = np.array(dataList)
    m = np.shape(dataArr)[0]  # 获取样本的数量
    xcord0 = []
    ycord0 = []
    xcord1 = []
    ycord1 = []
    # 对不同标记的数据进行分类
    for i in range(m):
        if int(labelList[i]) == 0:
            xcord0.append(dataArr[i, 1])
            ycord0.append(dataArr[i, 2])
        else:
            xcord1.append(dataArr[i, 1])
            ycord1.append(dataArr[i, 2])
    # 画散点图
    plt.scatter(xcord0, ycord0, s=30, c='red', marker='s')
    plt.scatter(xcord1, ycord1, s=30, c='green')
    # 画分割线
    x1 = np.arange(-3.0, 3.0, 0.1)
    x2 = (-weights[0] - weights[1] * x1) / weights[2]
    plt.plot(x1, x2)
    plt.xlabel('X1')
    plt.ylabel('X2')
    plt.show()

上述代码中，比较令人不解的是x2 = (-weights[0] - weights[1] * x1) / weights[2]这个是怎么来的？因为对于样本数据来说

是两个分类（类别1和类别0）的分界处，根据该式代入x_1即可求x_2从而这条线可以在图中表示出来。运行如下代码：

dataList, labelList = loadDataSet("testSet.txt") 
weights = gradAscent(dataList, labelList)  
plotBestFit(dataList, labelList, weights.getA())

输出结果如下图所示，从图中可以看到分类效果相当不错，从图上看值只分错了几个点。

然而跟线性回归中使用梯度下降算法求最小值类似，这种方法需要大量的计算，每次更新回归系数都需要遍历整个数据集，假如数据量一多，该方法的计算复杂度就很高了。所以我们对上述的批量梯度下降算法进行改变，改进方法是一次仅用一个样本点来更新回归系数，该方法称为随机梯度下降算法。下面是随机梯度下降算法的实现代码

'''
随机梯度下降算法
'''
def stocGradAscent0(dataList, labelList):
    dataArr = np.array(dataList)  # x数据转化为矩阵
    m, n = np.shape(dataArr)
    alpha = 0.01
    weights = np.ones(n)
    for i in range(m):
        h = sigmoid(np.sum(dataArr[i]*weights))
        error = (h-labelList[i])
        weights = weights-alpha*error*dataArr[i]
    return weights

运行如下代码，得到了该方法的分类效果图：

dataList, labelList = loadDataSet("testSet.txt")
weights = stocGradAscent0(dataList, labelList)
plotBestFit(dataList, labelList, weights)

相比批量梯度下降算法而言，效果较差。直接比较效果是有点不公平的，因为两者使用的数据集次数不同。而判断一个优化算法优劣的可靠方法是看它是否收敛，也就是说参数是否达到稳定值。通过下面这段代码，来绘制三个回归系数的变化情况：

'''
记录随机梯度下降过程中的历史权重值
'''
def stocGradAscent0WeightHistory(dataList, labelList):
    dataMat = np.array(dataList)  # x数据转化为矩阵
    m, n = np.shape(dataMat)
    alpha = 0.01
    times = 100
    weights = np.ones(n)
    weightsHistory = np.zeros((times*m, n) )
    for j  in range(times):
        for i in range(m):
            h = sigmoid(np.sum(dataMat[i] * weights))
            error = (h - labelList[i])
            weights = weights - alpha * error * dataMat[i]
            weightsHistory[j*m+i, :] = weights
    return weightsHistory
'''
绘制历史权值图
'''
def plotSDGError(weightsHistory):
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    fig = plt.figure()
    ax = fig.add_subplot(311)
    ax.plot(weightsHistory[:,0])
    plt.ylabel('X0')
    ax = fig.add_subplot(312)
    ax.plot(weightsHistory[:,1])
    plt.ylabel('X1')
    ax = fig.add_subplot(313)
    ax.plot(weightsHistory[:,2])
    plt.ylabel('X2')
    plt.xlabel(u'迭代次数')
    plt.show()


def main():
    dataList, labelList = loadDataSet("testSet.txt")
    weightsHistory = stocGradAscent0WeightHistory(dataList, labelList)
    plotSDGError(weightsHistory)

if __name__ == '__main__':
    main()

最终绘制的三个回归系数的变化情况如下：

从图中可以看到在大的波动之后,有些系数还会存在一些小的周期性波动。这种现象的原因是因为存在一些不能正确分类的样本点，在每次迭代时会引发系数的剧烈改变。我们希望算法能避免来回波动，从而收敛于某个值，同时是加快收敛速度，因此我们将上述的随机梯度下降算法进行更改。更改为如下：

'''
改进的随机梯度下降算法
'''
def stocGradAscent1(dataList, labelList, numIter=150):
    dataArr = np.array(dataList)
    m,n = np.shape(dataArr)
    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        # 步长为动态变化
            rand = int(np.random.uniform(0, len(dataIndex)))
            choseIndex = dataIndex[rand]
            h = sigmoid(np.sum(dataArr[choseIndex]*weights))
            error = h-labelList[choseIndex]
            weights = weights-alpha*error*dataArr[choseIndex]
            del(dataIndex[rand])
    return weights

相对第一个的随机梯度下降算法，上述算法一方面使alpha在每次迭代的时候都会进行调整，缓解了数据波动或者高频波动。虽然alpha会随着迭代次数不断减小，但是不会减少至0，而这样子做也保证了在多次迭代之后新数据仍然具有一定的影响。另一方面，通过随机选取样本来更新回归系数，从而减少周期性的波动。同时每次随机从列表中选出一个值之后，就从列表中删除该值。

上述算法对机器学习实战中的使用的算法有进行更改，作者个人认为书上在选出一个值并删除该值之后的实现方式不对。按照书上实现的算法，比如我们dataIndex序列是[0, 1, 2, 3, 4, 5]，那么假如第一次随机产生的数是0，那么使用的则是dataArr[0]这个样本，之后使用del删除了0下标这个数。那么dataIndex序列变成了[1, 2, 3, 4, 5]。之后再随机产生一个数，比如还是0，那么使用的还是dataArr[0]这个样本，之后del删除了下标为0的这个数，那么序列就变成了[2, 3, 4, 5]了。那么这种实现方式跟书的作者意思就有点不一致了。书中作者的意思应该是从[0, 1, 2, 3, 4, 5]这个列表中随机选取一个值，比如选取了0，然后使用第0个样本，之后删除0这个数，那么下次就变成从[1, 2, 3, 4, 5]这个列表中随机选取一个值，比如选取了1，然后使用第1个样本，之后删掉1这个数。

运行如下代码，得到的效果如下图所示，该效果达到了跟批量梯度下降算法差不多的效果，但是所使用的计算机更少

dataList, labelList = loadDataSet("testSet.txt")
weights = stocGradAscent1(dataList, labelList)
plotBestFit(dataList, labelList, weights)

使用更改之后的算法，来显示每次迭代时各个回归系数的变化情况，实现的代码如下

'''
记录改进的随机梯度下降过程中的历史权重值
'''
def stocGradAscent1WeightHistory(dataList, labelList, numIter=150):
    dataArr = np.array(dataList)
    m, n = np.shape(dataArr)
    weights = np.ones(n)
    weightsHistory = np.zeros((numIter*m, n))
    for j in range(numIter):
        dataIndex = list(range(m))
        for i in range(m):
            alpha = 4 / (1.0 + j + i) + 0.01  # 步长为动态变化
            rand = int(np.random.uniform(0, len(dataIndex)))
            choseIndex = dataIndex[rand]
            h = sigmoid(np.sum(dataArr[choseIndex] * weights))
            error = h - labelList[choseIndex]
            weights = weights - alpha * error * dataArr[choseIndex]
            weightsHistory[j*m+i, :] = weights
            del (dataIndex[rand])
    return weightsHistory

运行如下代码，得到的效果图如图所示。从图中可以看到没有出现周期性的波动，这归功于样本随机选择机制，同时收敛的速度也更快。

dataList, labelList = loadDataSet("testSet.txt")
weightsHistory = stocGradAscent1WeightHistory(dataList, labelList, 40)
plotSDGError(weightsHistory)

附：数据集及代码

本文所有数据集及代码，均可在https://github.com/DawnGuoDev/MachineLearningStudy/tree/master/LogisticRegression中获取。

本文参考与作者推荐

1. 机器学习，周志华

2. 机器学习实战，[美]Peter Harrington

3. CSDN_梯度下降法的向量化推导

4. CSDN_梯度上升算法与随机梯度上升算法的实现

不甘于「本该如此」，「多选参数 」值得关注

「多选参数 」的参数们正等着您加入，然后一起学习！