【机器学习】logistic回归模型

一、logistic回归概述

           我们都知道使用线性模型可以进行回归学习，但若要做的是分类任务改如何处理？只需要找一个单调可微函数将分类任务的真实标记y与线性回归模型的预测值联系起来。

        考虑二分类任务，其输出标记y的取值为0和1，而线性回归模型产生的预测值是实值，于是需将实值z转换为0/1值。通过Sigmoid函数引入非线性因素，可以实现实值z转换为0/1值，处理二分类问题。

（一）sigmoid函数

        在Logistic回归中，我们使用一个称为Logistic函数（也称为Sigmoid函数）的函数来建模概率。Logistic函数将输入值映射到[0, 1]之间的输出范围，表示数据属于正类的概率。该函数的数学表达式如下：

p(x) = 1 / (1 + exp(-z))

其中，p(x)表示数据属于正类的概率，x表示输入特征，z表示线性组合的参数和特征的乘积。

其函数曲线如下：

        Logistic回归通过最大似然估计来确定参数。其基本思想是，给定输入特征和对应的真实标签，我们寻找最佳参数值，使得模型对观测数据的分类概率与真实标签的概率分布最匹配。

        在训练阶段，我们使用训练数据来估计Logistic回归模型的参数。常用的方法是最大似然估计或梯度下降算法。在预测阶段，我们使用学习到的参数来计算新样本属于正类的概率，并根据设定的阈值进行分类决策。

（二）基于最优化方法的最佳回归系数确定

        1、极大似然估计

        极大似然估计是一种参数估计方法，用于确定给定模型的参数值，使得模型对观测数据的概率最大。在统计学中，它是一种常见的推断方法，应用广泛。

        假设我们有一些观测数据集合X，其概率分布为P(X | θ)，其中θ是待估计的参数。极大似然估计的目标是找到一个最优的θ，使得观测数据集合的概率P(X | θ)最大。

        具体来说，我们需要先定义模型的似然函数L(θ | X)，表示参数θ在给定观测数据X下的概率密度函数。似然函数与概率密度函数非常类似，但是概率密度函数是在已知参数情况下计算数据点的概率，而似然函数是在已知数据点情况下计算参数的概率。因此，似然函数可以表示为：

L(θ | X) = P(X | θ)

        针对似然函数L(θ | X)，我们需要找到一个最优的参数θ，使得L(θ | X)最大化。这可以通过求解似然函数的导数（也称为score函数）来实现。具体来说，我们需要找到一个θ，使得score函数等于0，即：

∂lnL(θ | X) / ∂θ = 0

其中，ln表示自然对数。这个方程的解就是最优的参数θ，它最大化了似然函数L(θ | X)。

极大似然估计的方法步骤：

• 确定待求解的未知参数，如均值、方差或特定分布函数的参数等；
• 计算每个样本的概率密度为；
• 根据样本的概率密度累乘构造似然函数：
• 通过似然函数最大化（求导为0），求解未知参数θ，为了降低计算难度，可采用对数加法替换概率乘法，通过导数为0/极大值来求解未知参数。

        极大似然估计是一种常见的参数估计方法。它具有许多优点，例如简单易用、渐进无偏、有效性高等。同时，它也存在一些问题，例如对异常值敏感、需要满足一些假设前提条件等。因此，在实际应用中需要结合具体问题和数据情况来选择合适的参数估计方法。

        2、梯度上升算法

        梯度上升法的目标是最大化似然函数。它首先计算似然函数对参数的梯度，然后按照该梯度的方向更新参数值。

具体步骤如下：

a. 初始化参数θ的值。

b. 计算似然函数关于参数θ的梯度。

c. 按照梯度的方向更新参数θ：θ = θ + α * gradient，其中α表示学习率，控制参数更新的步长。

d. 重复步骤b和c，直到满足终止条件（如达到最大迭代次数或梯度变化很小）。

        梯度上升法基本的思想是：要找到某函数的 最大值，最好的方法是沿着该函数的梯度方向探寻。如果梯度记为∇，则函数f(x,y)的梯度由下式表示：

        这个梯度意味着要沿x的方向移动 ，沿y的方向移动 。其中，函数f(x,y) 必须要在待计算的点上有定义并且可微。如下图： 

 梯度上升算法到达每个点后都会重新估计移动的方向。从P0开始，计算完该点的梯度，函数就根据梯度移动到下一点P1。在P1点，梯度再次被重新计算，并沿新的梯度方向移动到P2。如此循环迭代，直到满足停止条件。迭代的过程中，梯度算子总是保证我们能选取到最佳的移动方向。我们知道了移动的反向，那移动量的大小是多少。该量值称为步长，记做α。用向量来表示的话，梯度上升算法的迭代公式如下：

   

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

         3、梯度下降算法

梯度下降法： 梯度下降法的目标是最小化损失函数。与梯度上升法相反，它通过沿着梯度的反方向更新参数值来逐步减小损失函数。

具体步骤如下：

a. 初始化参数θ的值。

b. 计算损失函数关于参数θ的梯度。

c. 按照梯度的反方向更新参数θ：θ = θ - α * gradient，其中α表示学习率。

d. 重复步骤b和c，直到满足终止条件。

梯度下降算法，它与上述的梯度上升算法是一样的，只是公式中的加法需要变成减法。因此，对应的公式可以写成：

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

 

二、logistic实例分析

（一）准备数据

这次的数据同样与之前的一样下载，在testSet.txt文本中共有100个数据，每个数据包含两个数值型的特征X1和X2，与此同时还有其所属的类别，我们通过这个数据集来找出其最佳回归系数。所以我们可以编写代码先获取准备数据：

from numpy import *
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])])
        #将文本中的数据的前两个x1,x2和1.0作为x0一起放入列表中
        labelMat.append(int(lineArr[2]))
        #将标签放入标签列表
    return dataMat,labelMat
if __name__=='__main__':
	dataArr,labelMat = loadDataSet()
	print(dataArr)
    print(labelMat)

 
        这样我们就获取了文本中的数据，并且将其切分成我们所需要的形式。

        对于上述代码理解自然没有问题，但是其中数据列表中添加元素的时候，为什么除了数据本身具有的两个特征值，还加入了一个1.0作为X0，这个我们后面解释。

（二）使用梯度上升算法进行分类

我们前面提到Sigmoid函数的输入我们定义为Z，且Z为最佳系数与特征值的向量积，再结合我们之前推导出的梯度上升迭代公式，将迭代公式矢量化

其中Y是向量。我们就可以编写代码：

from numpy import *
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])])
        #将文本中的数据的前两个x1,x2和1.0作为x0一起放入列表中
        labelMat.append(int(lineArr[2]))
        #将标签放入标签列表
    return dataMat,labelMat
 
def sigmoid(inX):#sigmoid函数计算
    return 1.0/(1+exp(-inX))
 
def gradAscent(dataMatIn,classLabels):#梯度上升优化算法求最佳回归系数
    dataMatrix = mat(dataMatIn)
    #将训练特征数据转为numpy矩阵
    labelMat = mat(classLabels).transpose()
    #将标签列表转为numpy矩阵并将其转置
    m,n = shape(dataMatrix)
    #获得矩阵的行数和列数
    alpha = 0.001
    #向目标移动的步长，控制速率
    maxCycles = 500
    #迭代次数
    weights = ones((n,1))
    #初始化回归系数为1(以列数即特征数量为基准)
    for k in range(maxCycles):
        h = sigmoid(dataMatrix * weights)
        #梯度上升矢量化公式
        error = (labelMat - h)
        #公式计算
        weights = weights + alpha * dataMatrix.transpose()*error
    return weights.getA()
if __name__=='__main__':
    dataArr,labelMat = loadDataSet()
    weights = gradAscent(dataArr,labelMat)
    print(weights)

此处有所改动，为了方便后续的操作，所以gradAscent函数的返回值先将其从矩阵类型转为数组类型。

        由此我们便获得了其最佳回归系数。我们来解析代码，为什么要转置，我们清楚dataMatrix是mn的矩阵，而LabelMat是m1的矩阵，在参与Sigmoid函数运算时，dataMatrix与weights进行矩阵相乘，而weights是n*1的矩阵。我们知道矩阵相乘要满足前者列数要与后者的行数相等才能进行乘运算，所以要进行转置，所以在数据集列表前加入X0=1的一个原因就是方便计算。

（三）分析数据：画出决策边界 

        解出了回归系数，我们就确定了不同类别的数据之间的分隔线，那我们就可以将其画出来更深层次感受回归系数。用Python画函数图我们一样仍然用matplotlib。

from numpy import *
import matplotlib.pyplot as plt
 
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])])
        #将文本中的数据的前两个x1,x2和1.0作为x0一起放入列表中
        labelMat.append(int(lineArr[2]))
        #将标签放入标签列表
    return dataMat,labelMat
 
def sigmoid(inX):#sigmoid函数计算
    return 1.0/(1+exp(-inX))
 
def gradAscent(dataMatIn,classLabels):#梯度上升优化算法求最佳回归系数
    dataMatrix = mat(dataMatIn)
    #将训练特征数据转为numpy矩阵
    labelMat = mat(classLabels).transpose()
    #将标签列表转为numpy矩阵并将其转置
    m,n = shape(dataMatrix)
    #获得矩阵的行数和列数
    alpha = 0.001
    #向目标移动的步长，控制速率
    maxCycles = 500
    #迭代次数
    weights = ones((n,1))
    #初始化回归系数为1(以列数即特征数量为基准)
    for k in range(maxCycles):
        h = sigmoid(dataMatrix * weights)
        #梯度上升矢量化公式
        error = (labelMat - h)
        #公式计算
        weights = weights + alpha * dataMatrix.transpose()*error
    return weights.getA()
 
def plotBestFit(weights):#绘制数据集和Logistic回归最佳拟合直线
    dataMat,labelMat = loadDataSet()
    #获取数据集和标签集
    dataArr = array(dataMat)
    #将数据集转为array数组
    n = shape(dataArr)[0]
    #获取数据集的行数，即数据的个数
    xcord1 = []
    ycord1 = []
    xcord2 = []
    ycord2 = []
    for i in range(n):#遍历数据
        if int(labelMat[i]) == 1:#若分类标签为1
            xcord1.append(dataArr[i,1])
            ycord1.append(dataArr[i,2])
            #添加正样本
        else:#类别为0
            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',alpha=.5)
    ax.scatter(xcord2,ycord2,s=30,c="green",alpha=.5)
    # 绘制散点图，前两个是数组，s是标量，c是颜色设置，marker是点的形状设置(
    #默认是圆，即o，s是正方形),alpha是透明度，在0-1之间
    x = arange(-3.0,3.0,0.1)
    #与range区别是可以为小数，指从-3到3以间隔为0.1递增
    y = (-weights[0]-weights[1]*x)/weights[2]
    ax.plot(x,y)
    #绘制函数图线
    plt.xlabel('X1')
    plt.ylabel('X2')
    #绘制label
    plt.show()
    #展示
if __name__='__main__':
	dataArr,labelMat = loadDataSet()
    weights = gradAscent(array(dataArr),labelMat)
	plotBestFit(weights)

代码中，将标签为1的数据和标签为0的数据区分画出，并且以X1的特征值作为X轴的坐标，X2的特征值作为Y轴的坐标，并且确定拟合直线的长度范围及增长趋势。对于y = (-weights[0]-weights[1]*x)/weights[2]其实也很好理解。

我们之前提到过，Sigmoid函数是以0.5作为分界线的，概率小于0.5的我们划分为0类别，大于0.5的我们划分为1类别。而我们所要画的拟合直线就是区分两者类别的直线，而当Sigmoid函数为0.5时,x的取值就是为0，也就代表着Z=0，所以我们可以写出公式0=w0x0+w1x1+w2x2，w是最佳系数，x是特征值，由于X0都是1.0，所以公式可以写成0=w0+w1x1+w2x2，上述说到X1作为x轴，X2作为y轴，所以转换之后可以写成0=w0+w1x+w2y,再带入相应的回归参数，就是代码中的公式，在-3-3之间算出y的值，即X2。那我们也可以发现，W0最后就是区分列别划分线的截距，即y=ax+b中的b，所以之前加入X0=1.0就是匹配这个方程式，W0*1.0不会改变其值。

三、实验总结

       Logistic回归的目的是寻找一个非线性函数Sigmoid的最佳拟合参数，求解过程可以由最优化算法来完成。在最优化算法中，最常用的就是梯度上升算法，而梯度上升算法又可以简化为随机梯度上升算法。

        随机梯度上升算法与梯度上升算法的效果相当，但占用更少的计算资源。此外，随机梯度上升是一个在线算法，它可以在新数据到来时就完成参数更新，而不需要重新读取整个数据集来进行批处理运算。