Logistic回归-CSDN博客

本文链接：https://blog.csdn.net/huanghuali_/article/details/121463380

文章目录

一、什么是Logistic？
二.基于Logistic回归和Sigmoid函数的分类
三、基于最优化方法的最佳回归系数确定
四、从疝气病预测病马的死亡率
- 1.测试算法：用Logistic回归进行分类

提示：参考书籍图灵程序设计丛书《机器学习实战》

一、什么是Logistic？

logistic回归又称logistic回归分析，主要在流行病学中应用较多，比较常用的情形是探索某疾病的危险因素，根据危险因素预测某疾病发生的概率等等。例如，想探讨胃癌发生的危险因素，可以选择两组人群，一组是胃癌组，一组是非胃癌组，两组人群肯定有不同的体征和生活方式等。这里的因变量就是是否胃癌，即“是”或“否”，为两分类变量，自变量就可以包括很多了，例如年龄、性别、饮食习惯、幽门螺杆菌感染等。自变量既可以是连续的，也可以是分类的。通过logistic回归分析，就可以大致了解到底哪些因素是胃癌的危险因素。

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

我们想要的函数应该是: 能接受所有的输入然后预测出类别。例如，在两个类的情况下，上述函数输出 0 或 1.或许你之前接触过具有这种性质的函数，该函数称为海维塞得阶跃函数(Heaviside step function)，或者直接称为单位阶跃函数。然而，海维塞得阶跃函数的问题在于: 该函数在跳跃点上从 0 瞬间跳跃到 1，这个瞬间跳跃过程有时很难处理。幸好，另一个函数也有类似的性质（可以输出 0 或者 1 的性质），且数学上更易处理，这就是 Sigmoid 函数。

sigmoid函数表达式：
在这里插入图片描述
下图给出了 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，由下面公式得出：
在这里插入图片描述
如果采用向量的写法，上述公式可以写成z=wTx, ，它表示将这两个数值向量对应元素相乘然后全部加起来即得到 z 值。其中的向量 x 是分类器的输入数据，向量 w 也就是我们要找到的最佳参数（系数），从而使得分类器尽可能地精确。为了寻找该最佳参数，需要用到最优化理论的一些知识。我们这里使用的是——梯度上升法。

1.梯度上升法

梯度上升法基于的思想是：要找到某函数的最大值，最好的方法是沿着该函数的梯度方向探寻。如果梯度记为 ▽ ，则函数 f(x, y) 的梯度由下式表示：
在这里插入图片描述
这个梯度意味着要沿 x 的方向移动

沿 y 的方向移动

其中，函数f(x, y) 必须要在待计算的点上有定义并且可微。下图是一个具体的例子如下图所示。在这里插入图片描述
上图展示的，梯度上升算法到达每个点后都会重新估计移动的方向。从 P0 开始，计算完该点的梯度，函数就根据梯度移动到下一点 P1。在 P1 点，梯度再次被重新计算，并沿着新的梯度方向移动到 P2 。如此循环迭代，直到满足停止条件。迭代过程中，梯度算子总是保证我们能选取到最佳的移动方向。

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

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

testSet.txt网上能找到文本（下图数据只有一部分截图）
在这里插入图片描述

# logistic回归
import numpy as np
from numpy import *
import logRegres

#训练算法：使用梯度上升找到最佳参数
def loadDataSet():      #遍历函数：打开文件并逐行读取
    dataMat = [];labelMat = []
    fr = open('testSet.txt')        #前两个值为X1,X2，第三个为类别标签
    for line in fr.readlines():
        lineArr = line.strip().split()
        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])]) #为方便计算，将x0值设为1.0
        labelMat.append(int(lineArr[2]))
    return dataMat, labelMat

#定义sigmoid函数
def sigmoid(inX):
    return 1.0 / (1 + exp(-inX))

#梯度上升：dataMatIn:2维NumPy数组,100*3矩阵;classLabels:类别标签，1*100行向量
def gradAscent(dataMatIn, classLabels):
    dataMatrix = mat(dataMatIn)  # 原始数据转换为NumPy矩阵
    labelMat = mat(classLabels).transpose()  #类标签转换为矩阵：100*1列向量
    m, n = shape(dataMatrix)
    alpha = 0.001       #向目标移动的步长
    maxCycles = 500     #迭代次数
    weights = ones((n, 1))
    for k in range(maxCycles):  # heavy on matrix operations
        h = sigmoid(dataMatrix * weights)  # 矩阵相乘 h为列向量
        error = (labelMat - h)  # 真实类别与预测类别的差值
        weights = weights + alpha * dataMatrix.transpose() * error  #
    return weights

3.分析数据：画出决策边界

#分析数据：画出决策边界
def plotBestFit(wei):
    import matplotlib.pyplot as plt
    # weights=wei.getA()
    dataMat, labelMat = loadDataSet()
    dataArr = array(dataMat)
    n = 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])
        else:
            xcord2.append(dataArr[i, 1]);ycord2.append(dataArr[i, 2])
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(xcord1, ycord1, s=30, c='red', marker='s')
    ax.scatter(xcord2, ycord2, s=30, c='green')
    x = arange(-3.0, 3.0, 0.1)
    # y = (-weights[0] - weights[1] * x) /weights[2]  #最佳拟合直线
    # 设置了sigmoid函数为0,0是两个类别（类别1和类别0）的分界处。我们设定0=w0x0+w1x1+w2x2然后解出X1和X2的关系式
    y = (-wei[0] - wei[1] * x) / wei[2]  #最佳拟合直线
    ax.plot(x, y)
    plt.xlabel('X1');
    plt.ylabel('X2');
    plt.show()

if __name__ =='__main__':
    dataArr,labelMat=loadDataSet()
    weights=gradAscent(dataArr,labelMat)
    plotBestFit(weights.getA())

测试结果：
在这里插入图片描述

4.训练算法：随机梯度上升

# 随机梯度上升
def stocGradAscent0(dataMatrix, classLabels):
    m,n = shape(dataMatrix)
    alpha = 0.01
    weights = ones(n)   #initialize to all ones
    for i in range(m):
        h = sigmoid(sum(dataMatrix[i]*weights))
        error = classLabels[i] - h
        weights = weights + alpha * error * dataMatrix[i]
    return weights
    
if __name__ =='__main__':
    # 随机梯度上升
    dataArr, labelMat = loadDataSet()
    weights=stocGradAscent0(array(dataArr),labelMat)
    logRegres.plotBestFit(weights)

测试结果：
在这里插入图片描述

判断优化算法优劣的可靠方法是看它是否收敛，也就是说参数是否达到了稳定值，是否还会不断地变化？我们在程序清单中随机梯度上升算法上做了些修改，使其在整个数据集上运行 200 次，如下图所示:
在这里插入图片描述
运行随机梯度上升算法，在数据集的一次遍历中回归系数与迭代次数的关系图。回归系数经过大量迭代才能达到稳定值。并且仍然有局部的波动现象。

改进的随机梯度上升算法

#改进的随机梯度上升算法
def stoGradAscent1(dataMatrix, classLabels, numIter = 150):
    #m为行 n为列
    m,n = shape(dataMatrix)
    weights = ones(n)
    # 随机梯度, 循环150,观察是否收敛
    for j in range (numIter):
        # [0, 1, 2 .. m-1]
        dataIndex = range(m)
        for i in range(m):
            #alpha会随着迭代次数不断减小，但不会是0
            alpha = 4 / (1.0 + j + i) + 0.01            #alpha 每次迭代需要调整
            randIndex = int(random.uniform(0, len(dataIndex)))  #随机选取更新
            #随机选取更新
            h = sigmoid(sum(dataMatrix[randIndex] * weights))
            error = classLabels[randIndex] - h
            weights = weights + alpha * error * dataMatrix[randIndex]
            del (list(dataIndex)[randIndex])
    return weights
    
if __name__ =='__main__':
    # 改进的随机梯度上升算法
    dataArr, labelMat = loadDataSet()
    weights=stoGradAscent1(array(dataArr),labelMat)
    logRegres.plotBestFit(weights)

测试结果：
在这里插入图片描述

每次迭代时各个回归系数的变化情况
在这里插入图片描述
该方法比采用固定的alpha收敛速度更快。系数没有出现周期性的波动，主要归功于stocGradAscent1()的样本随机机制；stocGradAscent1()收敛更快。这次仅对数据集做了20次遍历，而之前的方法是500次。

四、从疝气病预测病马的死亡率

疝病是描述马胃肠痛的术语。然而，这种病不一定源自马的胃肠问题，其他问题也可能引发马疝病。这个数据集中包含了医院检测马疝病的一些指标，有的指标比较主观，有的指标难以测量，例如马的疼痛级别。

1.测试算法：用Logistic回归进行分类

horseColicTest.txt网上能找到文本（下图数据只有一部分截图）
在这里插入图片描述
horseColicTraining.txt网上能找到文本（下图数据只有一部分截图）

#从疝气病症预测病马的死亡率
#sigmoid函数分类
#以回归系数和特征向量作为输入来计算对应的Sigmoid值。
def classifyVector(inX, weights):
    prob = sigmoid(sum(inX*weights))            #sigmoid值大于0.5函数返回1，否则返回0
    if prob > 0.5: return 1.0
    else: return 0.0

#用于打开测试集和训练集
def colicTest():
    frTrain = open('horseColicTraining.txt')             # 训练数据
    frTest = open('horseColicTest.txt')                  # 测试数据
    trainingSet = []; trainingLabels = []
    for line in frTrain.readlines():
        currLine = line.strip().split('\t')
        lineArr =[]
        for i in range(21):
            lineArr.append(float(currLine[i]))
        trainingSet.append(lineArr)
        trainingLabels.append(float(currLine[21]))
    trainWeights = stoGradAscent1(np.array(trainingSet), trainingLabels, 1000) #通过随机梯度上升算法确定回归系数
    errorCount = 0; numTestVec = 0.0
    for line in frTest.readlines():
        numTestVec += 1.0
        currLine = line.strip().split('\t')
        lineArr =[]
        for i in range(21):
            lineArr.append(float(currLine[i]))
        if int(classifyVector(np.array(lineArr), trainWeights))!= int(currLine[21]):             #判断错误
            errorCount += 1
    errorRate = (float(errorCount)/numTestVec)
    print ("the error rate of this test is: %f" % errorRate)
    return errorRate

# 真实类别与预测类别的差值
def multiTest():
    numTests = 10; errorSum=0.0
    for k in range(numTests):
        errorSum += colicTest()
    print ("after %d iterations the average error rate is: %f" % (numTests, errorSum/float(numTests)))
   
   if __name__ =='__main__':
    # 从疝气病症预测病马的死亡率
    multiTest()

测试结果：平均错误率差不多30%
在这里插入图片描述

训练集中最后一列是类别标签。数据有两类类别标签：仍存活和未能存活（已经死亡and已经安乐死）。用stoGradAscent1()来计算回归系数向量。这里可以自由设定迭代的次数，例如在训练集上使用500次迭代，实验结果表明这比默认迭代150次的效果好。在系数计算完成之后，导入在测试集并计算分类错误率。整体看来，colicTest()具有完全独立的功能，多次运行得到的结果稍有不同，这是因为其中有随机成分在里面。如果stoGradAscent1()函数中回归系数已经完全收敛，那么结果才将是确定的。