监督学习--逻辑回归

Yolanda Yan 9

已于 2023-02-27 11:03:14 修改

阅读量865

点赞数

分类专栏：机器学习文章标签：机器学习 python 逻辑回归

于 2020-08-04 22:23:05 首次发布

本文链接：https://blog.csdn.net/Amy9_Miss/article/details/107801881

版权

机器学习专栏收录该内容

21 篇文章 5 订阅

订阅专栏

逻辑回归（Logistic Regression，LR）是分类常用的算法。逻辑回归在西瓜书中又被称为对数几率回归。进行分类的主要思想是：根据现有数据对分类边界线建立回归公式，以此进行分类。训练分类器时的做法就是寻找最佳拟合参数，使用的是最优化算法。

适合数据类型：数值型

优点：计算代价不高，易于理解和实现

缺点：容易欠拟合，分类精度可能不高

原理

正负类区分

正负类没有明确区分，但是按经验来说

负类（0）：一般表示没有什么，比如没有癌症，非垃圾邮件
正类（1）：一般表示具有我们要寻找的东西，比如有癌症，是垃圾邮件

基于逻辑回归的分类

对于二分类的情况，其输出值 $\epsilon {0,1}$ ，线性回归模型产生的预测值 $z=\omega^Tx$ ，其中x为输入数据，w是回归系数。我们需要将实值z转换为0/1值。最理想的是“单位阶跃函数”，具体如下：
$g(z)=\begin{cases} 0,\quad z<0\\ 0.5， z=0\\1,\quad z>0\end{cases}$
即若预测值z大于0就为1类，小于0就为0类。预测值为临界值0则可任意判断，如下图中红色部分。

可以看到，单位阶跃函数不连续，所以需要找到一个函数有类似性质，且数学上更易处理，这个函数就是sigmoid函数，它单调可微，见上图蓝色曲线，可以看到，它将z值转化为0~1范围内的值。sigmoid函数的计算公式如下：
$g(z)=\frac{1}{1+e^{-z}}$
从下图可以看到，sigmoid函数在不同坐标尺度下的两条曲线。如果横坐标刻度足够大，sigmoid函数看起来很像一个阶跃函数。

为了实现逻辑回归分类器，可以在每个特征上都乘以一个回归系数，然后把所有的结果相加，将这个总和结果代入到g(z)函数中得到一个范围在0~1之间的值。如果值大于0.5则分为1类，如果值小于0.5则分为0类。所以，g(z)函数可以看成是一种概率估计。g(z)为预测为1类的概率，1-g(z)为预测为0类的概率。

函数的表达形式如下：
$P(y=1|x,w)=g(z)\\P(y=0|x,w)=1-g(z)\\P(正确)=g(w,x_i)^{y^{(i)}}*(1-g(w,x_i))^{(1-y^{(i)})}$

最大似然估计：想要找到一组w，使得预测出的结果全部正确的概率最大，所有样本预测正确的概率相乘得到的P(总体正确)最大，则数学表达如下：
$\prod_{i=1}^{m}{g(w,x_i)^{y^{(i)}}*(1-g(w,x_i))^{(1-y^{(i)})}}$
目标就是求上式中的w，使得 $L (w)$ 值最大。由于连乘的函数不好计算，可以通过两边同时取log的形式让其变成连加，即如下式所示：
$\sum_{i=1}^{m}{(y^{(i)}ln(g(w,x_i))+(1-y^{(i)})}ln(1-g(w,x_i)))$
则上面的就是逻辑回归的损失函数，为了寻找最佳系数w，需要用到一些最优化算法，如梯度上升法、随机梯度上升法以及改进的随机梯度上升法等。

1. 梯度上升法(Gradient Ascent)

梯度上升法的思想是：要找到某函数的最大值，最好的方法就是沿着该函数的梯度方向探索。如果梯度记为 $\bigtriangledown$ ,则函数f(x,y)的梯度由下式表示
$\bigtriangledown f(x,y)= \left( \begin{array}{c} \frac{\partial f(x,y)}{\partial x} \\ \frac{\partial f(x,y)}{\partial y} \end{array}\right)$
这个梯度意味着要沿x方向移动 $\frac{\partial f(x,y)}{\partial x}$ ，沿y方向移动 $\frac{\partial f(x,y)}{\partial y}$ 。其中，函数f(x,y)必须要在待计算的点上有定义并且可微。具体见下图

求解梯度部分同样是对损失函数求偏导，过程如下:
$\begin{aligned}\bigtriangledown_{\omega}f(\omega)&=\frac{\partial l(w)}{\partial w_j}\\ &=\sum_{i=1}^{m}{(y^{(i)}\frac{1}{g(w,x_i)}\frac{\partial}{\partial w_j}g(w,x_i)-(1-y^{(i)})}\frac{1}{(1-g(w,x_i))}\frac{\partial}{\partial w_j}g(w,x_i))\\ &=\sum_{i=1}^{m}{(y^{(i)}\frac{1}{g(w^Tx_i)}-(1-y^{(i)})}\frac{1}{(1-g(w^Tx_i))})\frac{\partial }{\partial w_j}g(w^Tx_i)\\ &=\sum_{i=1}^{m}{(y^{(i)}\frac{1}{g(w^Tx_i)}-(1-y^{(i)})}\frac{1}{(1-g(w^Tx_i))})g(w^Tx_i)(1-g(w^Tx_i))\frac{\partial }{\partial w_j}w^Tx_i\\ &=\sum_{i=1}^{m}{(y^{(i)}(1-g(w^Tx_i))-(1-y^{(i)})}g(w^Tx_i))x_i\\ &=\sum_{i=1}^{m}{(y^{(i)}-g(w^Tx_i))x_i}\\ \end{aligned}$

梯度上升法达到每个点后都会重新估计移动的方向。梯度算子总是指向函数值增长最快的方向。这里说的是移动方向，未提到移动量的大小，该量值称为步长，记为 $\alpha$ ，梯度上升的迭代公式如下
$\omega := \omega + \alpha \bigtriangledown_{\omega}f(\omega)$
该公式将一直被迭代下去，直到达到某个停止条件为止，比如迭代次数达到某个指定值或算法达到某个可以允许的误差范围。

梯度下降法(Gradient Descent)

经常会听到的应该是梯度下降法，它和这里的梯度上升法一样的，只是公式中的加法需要改为减法。因此，梯度下降法的迭代公式如下
$\omega := \omega - \alpha \bigtriangledown_{\omega}f(\omega)$
梯度上升法用来求函数的最大值，而梯度下降法用来求函数的最小值。

Python实现

def sigmoid(inX):
    """返回sigmoid后的值"""
    return 1.0/(1+np.exp(-inX))

def gradAscent(dataSet,labels):
    """用最优化算法（梯度上升）获取最优回归系数"""
    dataMat = np.mat(dataSet)
    labelsMat = np.mat(labels).transpose()  # 将列表转换为列向量，即 n*1
    m,n = dataMat.shape
    numIters = 500  # 迭代次数
    alpha = 0.01   # 移动步长
    weights = np.ones((n,1))
    for i in range(numIters):
        h = sigmoid(dataMat*weights)  # 计算预测类别
        error = labelsMat-h  # 计算真实类别和预测类别的差值
        # 按照差值的方向调整回归系数，w= w + alpha*(y-predict)*data
        weights = weights+alpha*dataMat.transpose()*error   
    return weights

给出数据集和绘制最佳拟合直线的python实现，具体如下

def loadDataSet():
    """导入数据"""
    data = []
    labels = []
    fr = open('testSet.txt')
    for line in fr.readlines():
        currLine = line.strip().split()
        data.append([1.0, float(currLine[0]), float(currLine[1])])
        labels.append(int(currLine[2]))
    return data, labels

def plotBestFit(weights):
    """画出数据集和最佳拟合直线"""
    import matplotlib.pyplot as plt
    dataSet, labels = loadDataSet()
    dataMat = np.mat(dataSet)
    m = dataMat.shape[0]
    xcoord1 = []
    ycoord1 = []
    xcoord2 = []
    ycoord2 = []
    for i in range(m):
        if int(labels[i]) == 1:
            xcoord1.append(dataMat[i, 1])
            ycoord1.append(dataMat[i, 2])
        else:
            xcoord2.append(dataMat[i, 1])
            ycoord2.append(dataMat[i, 2])
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(xcoord1, ycoord1, s=30, c='red', marker='s')
    ax.scatter(xcoord2, ycoord2, s=30, c='green', marker='o')
    x = np.arange(-3.0, 3.0, 0.1)
    # 最佳拟合直线：对于h(z)=1/(1+e^(-z))函数，h(z)=0.5是两个分类的边界线，即z=0时
    # 则有z=w^T*x=0 -> w0*1+w1*x1+w2*x2=0 -> x2=(-w0-w1*x1)/w2
    y = (-weights[0] - weights[1] * x) / weights[2]
    ax.plot(x, y)
    plt.xlabel('X1')
    plt.ylabel('X2')
    plt.show()

运行下面代码

dataSet, labels = loadDataSet()
weights = gradAscent(np.array(dataSet), labels)
plotBestFit(np.array(weights))

得到基于梯度上升法的最佳拟合直线图如下

可以看到，梯度上升法迭代了500次后的效果相当不错，只错分了几个点。但是这个方法计算量很大。

2. 随机梯度上升法 SGA (Stochastic Gradient Ascent)

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

Python实现

def stocGradAscent0(dataMat, labels):
    """用最优化算法（随机梯度上升法）获取最优回归系数"""
    m, n = dataMat.shape
    alpha = 0.01
    weights = np.ones(n)
    for i in range(m):
        h = sigmoid(sum(dataMat[i] * weights))
        error = labels[i] - h
        weights = weights + alpha * error * dataMat[i]  # alpha和error都是数值而非向量
    return weights

可以看到，随机梯度上升算法与梯度上升算法在代码上很相似，但也有一些区别：第一，后
者的变量h和误差error都是向量，而前者则全是数值；第二，前者没有矩阵的转换过程，所有
变量的数据类型都是N UmPy数组。

运行下面代码

dataSet, labels = loadDataSet()
weights = stocGradAscent0(np.array(dataSet), labels)
plotBestFit(np.array(weights))

得到基于随机梯度上升法的最佳拟合直线图如下

可以看到拟合效果没有基于梯度上升法的最佳拟合直线图那样完美，这里的分类器错分了三分之一的样本。但是直接这么比较不公平，因为基于梯度上升法的最佳拟合直线图是在整个数据集上迭代了500次才得到的效果。

一个判断优化算法优劣的可靠方法是看它是否收敛，也就是说参数是否达到稳定值，是否还会不断地变化？对此，对随机梯度上升法在整个数据集上迭代200次，绘制出三个回归系数的变化情况如下

上图为随机梯度上升算法在200次迭代过程中回归系数的变化情况。可以看到回归系数经过大量的迭代才能达到稳定值，而且仍然有局部的波动现象。另外需要注意的是，在大的波动停止后，还有一些小的周期性波动。产生这种现象的原因是存在一些不能正确分类的样本点（数据集并非线性可分），在每次迭代时会引发系数的剧烈改变，我们期望算法能避免来回波动，从而收敛到某个值。另外，收敛速度也需要加快。

3. 改进的随机梯度上升法 SGA

第一个改进点：alpha在每次迭代的时候都会调整，这会缓解数据的波动或者高频波动。并且，虽然alpha会随着迭代次数不断减小，但永远不会减小到0，因为这里还存在一个常数项。必须这样做的原因是为了保证在多次迭代之后新数据仍然具有一定的影响。如果需要处理的问题是动态变化的，那么可以适当加大上述常数项，来确保新的值获得更大的回归系数。

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

Python实现

def stocGradAscent1(dataMat, labels, numIters=150):
    """用最优化算法（改进的随机梯度上升）获取最优回归系数"""
    m, n = dataMat.shape
    weights = np.ones(n)
    for j in range(numIters):
        dataIndex = list(np.arange(m))
        for i in range(m):
            alpha = 4.0 / (j + i + 1.0) + 0.01  # 每次迭代时需要调整alpha,缓解数据波动
            # 随机选取样本来更新回归系数
            randIndex = int(np.random.uniform(0, len(dataIndex)))
            h = sigmoid(sum(dataMat[randIndex] * weights))
            error = labels[randIndex] - h
            weights = weights + alpha * error * dataMat[randIndex]
            del dataIndex[randIndex]

    return weights

运行下面代码

dataSet, labels = loadDataSet()
weights = stocGradAscent1(np.array(dataSet), labels)
plotBestFit(np.array(weights))

得到基于改进的随机梯度上升法的最佳拟合直线图如下

从上图可以发现，改进后的梯度上升法没有出现周期性的波动，这归功于样本随机选择机制。而且改进后的梯度上升法仅仅对数据集做了20次遍历，回归系数就收敛了，收敛速度更快。

示例：从疝气病症预测病马的死亡率

使用逻辑回归来预测患有疝病的马的存活问题。

数据情况

该数据集包含300个训练样本和68个测试样本，每个样本包含了医院检测马疝病的一些指标（共计28个），有的指标比较主观，有的指标难以测量，例如马的疼痛级别。
数据集来源点这里下载，下载下图这两个文件

需要把.data和.test文件里的问号替换为0，否则后面的代码会报错

将处理完的.data保存为horseColicTrain.txt，.test保存为horseColicTest.txt，具体如下

数据预处理

预处理1：需要将所有的缺失值用一个实数值来替换，这里选择实数0来替换所有缺失值，恰好能适用于逻辑回归。这是因为它在更新时不会影响权重系数的值，而且sigmoid (0) = 0.5，它对结果的预测不具有任何倾向性，因此上述做法也不会对误差项造成任何影响。基于上述原因，将缺失值用0代替既可以保留现有数据，也不需要对优化算法进行修改。此外，该数据集中的特征取值一般不为0 , 因此在某种意义上说它也满足“特殊值”这个要求。

预处理2：如果在测试数据集中发现了一条数据的类别标签已经缺失，那么我们的简单做法是将该条数据丢弃。这是因为类别标签与特征不同，很难确定采用某个合适的值来替换。

Python实现

数据最初有三个类别标签，分别代表马的三种情况：“仍存活”、“已经死亡”和“已经安乐死”。这里为了方便，将“已经死亡”和“已经安乐死”合并成“未能存活”这个标签。

def classifyVector(inX, weights):
    """
    逻辑回归的分类函数
    :param inX: 某条给定数据
    :param weights: 最优回归系数
    :return: 预测的分类值
    """
    prob = sigmoid(sum(inX * weights))
    if prob > 0.5:
        return 1.0  # 表示未能存活
    else:
        return 0.0  # 表示仍存活


def colicTest():
    """导入训练和测试集，进行逻辑回归，计算测试集分类错误率"""
    frTrain = open('horseColicTrain.txt')  # 训练集
    frTest = open('horseColicTest.txt')    # 测试集
    # 获取训练集输入数据和标签
    trainingSet = []
    trainingLabels = []
    for line in frTrain.readlines():
        currLine = line.strip().split(' ')  # 以空格为分隔符
        lineArr = []
        del currLine[2]  # 删除数据第三列，即医院登记号
        for i in range(21):  # 删除后数据第22列表示马活着、死掉或者安乐死
            lineArr.append(float(currLine[i]))
        trainingSet.append(lineArr)
        # 对已经安乐死和已经死亡归类于未能存活（这里是正例，用1表示；存活认为是反例用0表示）
        if int(float(currLine[21])) == 2 or int(float(currLine[21])) == 3:
            trainingLabels.append(1.0)
        else:
            trainingLabels.append(0.0)
    # 对训练集数据调用随机梯度上升算法，得到最优回归系数
    trainingweights = stocGradAscent1(np.array(trainingSet), trainingLabels,500)
    # 获取测试集，并计算分类错误数
    errorCount = 0
    numTestVec = 0.0
    for line in frTest.readlines():
        numTestVec += 1
        currLine = line.strip().split(' ')
        lineArr = []
        del currLine[2]
        for i in range(21):
            lineArr.append(float(currLine[i]))
        if int(float(currLine[21])) == 2 or int(float(currLine[21])) == 3:
            curlabel = 1.0
        else:
            curlabel = 0.0
        if int(classifyVector(np.array(lineArr), trainingweights)) != int(
                curlabel):
            errorCount += 1
    # 计算测试集的分类错误率
    errorRate = float(errorCount) / numTestVec
    print('the error rate of this test is: %f' % (errorRate))
    return errorRate


def multiTest():
    """计算迭代10次的平均错误率"""
    numTest = 10
    errorSum = 0.0
    for k in range(numTest):
        errorSum += colicTest()
    print('after %d iterations the average error rate is: %f' % (
    numTest, errorSum / float(numTest)))