赞
踩
导读:逻辑回归(Logistic regression)即逻辑模型,属于常见的一种分类算法。本文将从理论介绍开始,搞清楚什么是逻辑回归、回归系数、算法思想、工作原理及其优缺点等。进一步通过两个实际案例深化理解逻辑回归,以及在工程应用进行实现。(本文原创,转载必须注明出处: 决策树模型算法研究与案例分析)
逻辑回归
回归:假设现在有一些数据点,我们用一条直线对这些点进行拟合(这条直线称为最佳拟合直线),这个拟合的过程就叫做回归。
逻辑回归(Logistic Regression)是一种用于解决二分类(0 or 1)问题的机器学习方法,用于估计某种事物的可能性。比如某用户购买某商品的可能性,某病人患有某种疾病的可能性,以及某广告被用户点击的可能性等。 注意,这里用的是“可能性”,而非数学上的“概率”,logisitc回归的结果并非数学定义中的概率值,不可以直接当做概率值来用。该结果往往用于和其他特征值加权求和,而非直接相乘。
Sigmoid 函数
Sigmoid函数是一个常见的S型数学函数,在信息科学中,由于其单增以及反函数单增等性质,Sigmoid函数常被用作神经网络的阈值函数,将变量映射到0,1之间。在逻辑回归、人工神经网络中有着广泛的应用。Sigmoid函数的数学形式是:
对x求导可以推出如下结论:
下图给出了 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 函数得出的值为0.5,可以理解为给定数据和参数,数据被分入 1 类的概率为0.5。(注意:针对二分类问题,0.5不是唯一确定分类的值,你可以根据需求调整这个概率值。)
逻辑回归与线性回归的关系
逻辑回归(Logistic Regression)与线性回归(Linear Regression)都是一种广义线性模型(generalized linear model)。逻辑回归假设因变量 y 服从伯努利分布,而线性回归假设因变量 y 服从高斯分布。 因此与线性回归有很多相同之处,去除Sigmoid映射函数的话,逻辑回归算法就是一个线性回归。可以说,逻辑回归是以线性回归为理论支持的,但是逻辑回归通过Sigmoid函数引入了非线性因素,因此可以轻松处理0/1分类问题。
Sigmoid 函数的输入记为 z ,由下面公式得到:
如果采用向量的写法,上述公式可以写成 Sigmoid 函数计算公式向量形式 ,它表示将这两个数值向量对应元素相乘然后全部加起来即得到 z 值。其中的向量 x 是分类器的输入数据,向量 w 也就是我们要找到的最佳参数(系数),从而使得分类器尽可能地精确。为了寻找该最佳参数,需要用到最优化理论的一些知识。我们这里使用的是——梯度上升法(Gradient Ascent)。
梯度
对于梯度这个词一些人比较陌生,我们先看看维基百科的解释:在向量微积分中,标量场(向量场)中某一点的梯度指向在这点标量场增长最快的方向(当然要比较的话必须固定方向的长度),梯度的绝对值是长度为1的方向中函数最大的增加率,也就是说,其中 代表方向导数。
梯度形式化描述
考虑一座高度在 (x, y)点是 H(x, y)的山。H这一点的梯度是在该点坡度(或者说斜度)最陡的方向。梯度的大小告诉我们坡度到底有多陡。这个现象可以如下数学的表示。山的高度函数 H的梯度点积一个单位向量给出了表面在该向量的方向上的斜率。这称为方向导数。
理解梯度
为了大家更容易理解什么是梯度,我们介意向量的概念,向量是一个矢量具有大小和方向的。同样,梯度也可以类比为具备大小和方向的这么一个概念。其两者比较如下:(这里严格意义上讲是不成立的,便于大家理解。)
- 向量 = 值 + 方向
- 梯度 = 向量
- 梯度 = 梯度值 + 梯度方向
梯度上升
要找到某函数的最大值,最好的方法是沿着该函数的梯度方向探寻。如果梯度记为 ▽ ,则函数 f(x, y) 的梯度由下式表示:
这个梯度意味着要沿 x 的方向移动 ,沿 y 的方向移动。其中,函数f(x, y) 必须要在待计算的点上有定义并且可微。下图是一个具体的例子。
上图展示的,梯度上升算法到达每个点后都会重新估计移动的方向。从 P0 开始,计算完该点的梯度,函数就根据梯度移动到下一点 P1。在 P1 点,梯度再次被重新计算,并沿着新的梯度方向移动到 P2 。如此循环迭代,直到满足停止条件。迭代过程中,梯度算子总是保证我们能选取到最佳的移动方向。
上图中的梯度上升算法沿梯度方向移动了一步。可以看到,梯度算子总是指向函数值增长最快的方向。这里所说的是移动方向,而未提到移动量的大小。该量值称为步长,记作 α 。用向量来表示的话,梯度上升算法的迭代公式如下:
- 例如:y = w0 + w1x1 + w2x2 + ... + wnxn
- 梯度:参考上图的例子,二维图像,x方向代表第一个系数,也就是 w1,y方向代表第二个系数也就是 w2,这样的向量就是梯度。
- α:上面的梯度算法的迭代公式中的阿尔法,这个代表的是移动步长(step length)。移动步长会影响最终结果的拟合程度,最好的方法就是随着迭代次数更改移动步长。步长通俗的理解,100米,如果我一步走10米,我需要走10步;如果一步走20米,我只需要走5步。这里的一步走多少米就是步长的意思。
- ▽f(w):代表沿着梯度变化的方向,也可以理解该方向求导。
该公式将一直被迭代执行,直至达到某个停止条件为止,比如迭代次数达到某个指定值或者算法达到某个可以允许的误差范围。
梯度上升与梯度下降的区别
梯度下降是大家听的最多的,本质上梯度下降与梯度上升算法是一样的,只是公司中加法变减法,梯度下降的公式如下:
在求极值的问题中,有梯度上升和梯度下降两个最优化方法。梯度上升用于求最大值,梯度下降用于求最小值。如logistic回归的目标函数:代表的是概率,我们想求概率最大值,即对数似然函数的最大值,所以使用梯度上升算法。而线性回归的代价函数:代表的则是误差,我们想求误差最小值,所以用梯度下降算法。
根据现有数据对分类边界建立回归公司,以此进行分类。回归即最佳拟合。
- 每个回归系数初始化为 1
- 重复 R 次:
- 计算整个数据集的梯度
- 使用 步长 x 梯度 更新回归系数的向量
- 返回回归系数
- 收集数据: 采用任意方法收集数据
- 准备数据: 由于需要进行距离计算,因此要求数据类型为数值型。另外,结构化数据格式则最佳。
- 分析数据: 采用任意方法对数据进行分析。
- 训练算法: 大部分时间将用于训练,训练的目的是为了找到最佳的分类回归系数。
- 测试算法: 一旦训练步骤完成,分类将会很快。
- 使用算法: 首先,我们需要输入一些数据,并将其转换成对应的结构化数值;接着,基于训练好的回归系数就可以对这些数值进行简单的回归计算,判定它们属于哪个类别;
- 优点: 计算代价不高,易于理解和实现。
- 缺点: 容易欠拟合,分类精度可能不高。
- 适用数据类型: 数值型和标称型数据。
在一个简单的数据集上,采用梯度上升法找到 Logistic 回归分类器在此数据集上的最佳回归系数
开发流程
- 收集数据: 可以使用任何方法
- 准备数据: 由于需要进行距离计算,因此要求数据类型为数值型。另外,结构化数据格式则最佳
- 分析数据: 画出决策边界
- 训练算法: 使用梯度上升找到最佳参数
- 测试算法: 使用 Logistic 回归进行分类
- 使用算法: 对简单数据集中数据进行分类
本文采用100行的测试集文本。其中前两列是特征1,和特征2,第三类是对应的标签。(这里特征1,特征2作为测试使用没有实际意义,你可以理解为特征1 是水里游的,特征2是有鱼鳞。类别判断是否为鱼类。)
读取文本文件,加载数据集和类标签,这里将特征集第一列加1,便于后续回归系数的计算:
- '''加载数据集和类标签'''
- def loadDataSet(file_name):
- # dataMat为原始数据, labelMat为原始数据的标签
- dataMat,labelMat = [],[]
- fr = open(file_name)
- for line in fr.readlines():
- lineArr = line.strip().split(',')
- if len(lineArr) == 1:
- continue # 这里如果就一个空的元素,则跳过本次循环
- # 为了方便计算,我们将每一行的开头添加一个 1.0 作为 X0
- dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
- labelMat.append(int(lineArr[2]))
- return dataMat, labelMat
梯度上升算法
使用梯度上升训练算法模型,其代码实现如下:
- ''' 正常的梯度上升法,得到的最佳回归系数 '''
- def gradAscent(dataMatIn, classLabels):
- dataMatrix = mat(dataMatIn) # 转换为 NumPy 矩阵
- # 转化为矩阵[[0,1,0,1,0,1.....]],并转制[[0],[1],[0].....]
- # transpose() 行列转置函数
- # 将行向量转化为列向量 => 矩阵的转置
- labelMat = mat(classLabels).transpose() # 首先将数组转换为 NumPy 矩阵,然后再将行向量转置为列向量
- # m->数据量,样本数 n->特征数
- m, n = shape(dataMatrix) # 矩阵的行数和列数
- # print(m,n)
- alpha = 0.001 # alpha代表向目标移动的步长
- maxCycles = 500 # 迭代次数
- weights = ones((n, 1)) # 代表回归系数,ones((n,1)) 长度和特征数相同矩阵全是1
- for k in range(maxCycles):
- h = sigmoid(dataMatrix * weights) # 矩阵乘法
- # labelMat是实际值
- error = (labelMat - h) # 向量相减
- # 0.001* (3*m)*(m*1) 表示在每一个列上的一个误差情况,最后得出 x1,x2,xn的系数的偏移量
- weights = weights + alpha * dataMatrix.transpose() * error # 矩阵乘法,最后得到回归系数
- return array(weights)
其中sigmoid函数实现如下:
- ''' sigmoid跳跃函数 '''
- def sigmoid(ZVar):
- return 1.0 / (1 + exp(-ZVar))
代码分析:函数的两个参数是数据加载返回的特征集和标签类集合。对数据集进行mat矩阵话转化,而类标签集进行矩阵之后转置,便于行列式的计算。然后设定步长,和迭代次数。整个特征矩阵与回归系数乘积求sigmoid值,最后返回回归系数的值。运行结果如下:
- [[ 4.12414349]
- [ 0.48007329]
- [-0.6168482 ]]
思考?步长和迭代次数的初始值如何设定?
随机梯度上升算法
梯度上升算法在每次更新回归系数时都需要遍历整个数据集,该方法在处理 100 个左右的数据集时尚可,但如果有数十亿样本和成千上万的特征,那么该方法的计算复杂度就太高了。一种改进方法是一次仅用一个样本点来更新回归系数,该方法称为 随机梯度上升算法。由于可以在新样本到来时对分类器进行增量式更新,因而随机梯度上升算法是一个在线学习(online learning)算法。与 “在线学习” 相对应,一次处理所有数据被称作是 “批处理” (batch) 。其伪代码是:
- 所有回归系数初始化为 1
- 对数据集中每个样本
- 计算该样本的梯度
- 使用 alpha x gradient 更新回归系数值
- 返回回归系数值
随机梯度上升算法的代码实现如下:
- ''' 随机梯度上升'''
- # 梯度上升与随机梯度上升的区别?梯度下降在每次更新数据集时都需要遍历整个数据集,计算复杂都较高;随机梯度下降一次只用一个样本点来更新回归系数
- def stocGradAscent0(dataMatrix, classLabels):
- m, n = shape(dataMatrix)
- alpha = 0.01
- weights = ones(n) # 初始化长度为n的数组,元素全部为 1
- for i in range(m):
- # sum(dataMatrix[i]*weights)为了求 f(x)的值,f(x)=a1*x1+b2*x2+..+nn*xn
- h = sigmoid(sum(dataMatrix[i] * weights))
- # 计算真实类别与预测类别之间的差值,然后按照该差值调整回归系数
- error = classLabels[i] - h
- # 0.01*(1*1)*(1*n)
- weights = array(weights) + alpha * error * array(mat(dataMatrix[i]))
- return array(weights.transpose())
可以看到,随机梯度上升算法与梯度上升算法在代码上很相似,但也有一些区别: 第一,后者的变量 h 和误差 error 都是向量,而前者则全是数值;第二,前者没有矩阵的转换过程,所有变量的数据类型都是 NumPy 数组。
判断优化算法优劣的可靠方法是看它是否收敛,也就是说参数是否达到了稳定值,是否还会不断地变化?下图展示了随机梯度上升算法在 200 次迭代过程中回归系数的变化情况。其中的系数2,也就是 X2 只经过了 50 次迭代就达到了稳定值,但系数 1 和 0 则需要更多次的迭代。如下图所示:
针对波动问题,我们改进了之前的随机梯度上升算法,具体代码实现如下:
- ''' 改进版的随机梯度上升,使用随机的一个样本来更新回归系数'''
- def stocGradAscent1(dataMatrix, classLabels, numIter=150):
- m, n = shape(dataMatrix)
- weights = ones(n) # 创建与列数相同的矩阵的系数矩阵
- # 随机梯度, 循环150,观察是否收敛
- for j in range(numIter):
- dataIndex = list(range(m)) # [0, 1, 2 .. m-1]
- for i in range(m):
- # i和j的不断增大,导致alpha的值不断减少,但是不为0
- alpha = 4 / (1.0 + j + i) + 0.0001 # alpha随着迭代不断减小非0
- # random.uniform(x, y) 随机生成下一个实数,它在[x,y]范围内
- Index = int(random.uniform(0, len(dataIndex)))
- # sum(dataMatrix[i]*weights)为了求 f(x)的值, f(x)=a1*x1+b2*x2+..+nn*xn
- h = sigmoid(sum(dataMatrix[dataIndex[Index]] * weights))
- error = classLabels[dataIndex[Index]] - h
- weights = weights + alpha * error *array(mat(dataMatrix[dataIndex[Index]]))
- del (dataIndex[Index])
- # print(weights.transpose())
- return weights.transpose()
上面的改进版随机梯度上升算法改了两处代码。
边界可视化的代码实现如下:
- ''' 数据可视化展示 '''
- def plotBestFit(dataArr, labelMat, weights):
- n = shape(dataArr)[0]
- xcord1,xcord2,ycord1,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)
- """
- dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
- w0*x0+w1*x1+w2*x2=f(x)
- x0最开始就设置为1, x2就是我们画图的y值,而f(x)被我们磨合误差给算到w0,w1,w2身上去了
- 所以: w0+w1*x+w2*y=0 => y = (-w0-w1*x)/w2
- """
- y = (-weights[0] - weights[1] * x) / weights[2]
- ax.plot(x, y)
- plt.xlabel('X')
- plt.ylabel('Y')
- plt.show()
运行结果分别是:
梯度上升算法可视化结果图1-1:
随机梯度上升算法可视化结果:
优化随机梯度上升算法可视化结果:
结果分析:
图1-1的梯度上升算法在每次更新回归系数时都需要遍历整个数据集,虽然分类结果还不错该方法的计算复杂度就太高了。图1-2的随机梯度上升算法虽然分类效果不是很好(分类1/3左右),但是其迭代次数远远小于图1-1迭代次数(500次)。整体性能有所改进,但是其存在局部波动现象。基于此改进后的图1-3效果显示好很多。
代码实现如下:
- '''数据集决策可视化'''
- def simpleTest(file_name):
- # 1.收集并准备数据
- dataMat, labelMat = loadDataSet(file_name)
- # 2.训练模型, f(x)=a1*x1+b2*x2+..+nn*xn中 (a1,b2, .., nn).T的矩阵值
- dataArr = array(dataMat)
- weights = stocGradAscent1(dataArr, labelMat)
- # 数据可视化
- plotBestFit(dataArr, labelMat, weights)
使用 Logistic 回归来预测病毒性流感预测病人的死亡问题。这个数据集中包含了医院检测病毒性流感的一些指标,有的指标比较主观,有的指标难以测量,例如人的疼痛级别。
开发流程
收集数据: 给定数据文件
准备数据: 用 Python 解析文本文件并填充缺失值
分析数据: 可视化并观察数据
训练算法: 使用优化算法,找到最佳的系数
测试算法: 为了量化回归的效果,需要观察错误率。根据错误率决定是否回退到训练阶段,
通过改变迭代的次数和步长的参数来得到更好的回归系数
使用算法: 实现一个简单的命令行程序来收集马的症状并输出预测结果并非难事,
这可以作为留给大家的一道习题
训练数据已经给出,这里对文件处理即可,代码如下:
- '''加载数据集和类标签2'''
- def loadDataSet2(file_name):
- frTrain = open(file_name)
- trainingSet,trainingLabels = [],[]
- for line in frTrain.readlines():
- currLine = line.strip().split(',')
- # print(len(currLine))
- lineArr = []
- for i in range(len(currLine)-1):
- lineArr.append(float(currLine[i]))
- trainingSet.append(lineArr)
- trainingLabels.append(float(currLine[len(currLine)-1]))
- return trainingSet,trainingLabels
处理数据中的缺失值
假设有100个样本和20个特征,这些数据都是机器收集回来的。若机器上的某个传感器损坏导致一个特征无效时该怎么办?此时是否要扔掉整个数据?这种情况下,另外19个特征怎么办? 它们是否还可以用?答案是肯定的。因为有时候数据相当昂贵,扔掉和重新获取都是不可取的,所以必须采用一些方法来解决这个问题。下面给出了一些可选的做法:
现在,我们对下一节要用的数据集进行预处理,使其可以顺利地使用分类算法。在预处理需要做两件事:所有的缺失值必须用一个实数值来替换,因为我们使用的 NumPy 数据类型不允许包含缺失值。我们这里选择实数 0 来替换所有缺失值,恰好能适用于 Logistic 回归。这样做的直觉在于,我们需要的是一个在更新时不会影响系数的值。回归系数的更新公式如下:
weights = weights + alpha * error * dataMatrix[dataIndex[randIndex]]
如果 dataMatrix 的某个特征对应值为 0,那么该特征的系数将不做更新,即:weights = weights
另外,由于 Sigmoid(0) = 0.5 ,即它对结果的预测不具有任何倾向性,因此我们上述做法也不会对误差造成任何影响。基于上述原因,将缺失值用 0 代替既可以保留现有数据,也不需要对优化算法进行修改。此外,该数据集中的特征取值一般不为 0,因此在某种意义上说它也满足 “特殊值” 这个要求。
如果在测试数据集中发现了一条数据的类别标签已经缺失,那么我们的简单做法是将该条数据丢弃。这是因为类别标签与特征不同,很难确定采用某个合适的值来替换。采用 Logistic 回归进行分类时这种做法是合理的,而如果采用类似 kNN 的方法,则保留该条数据显得更加合理。
训练算法模型代码如下:
- '''测试Logistic算法分类'''
- def testClassier():
- # 使用改进后的随机梯度上升算法 求得在此数据集上的最佳回归系数 trainWeights
- file_name = './HorseColicTraining.txt'
- trainingSet,trainingLabels = loadDataSet2(file_name)
- trainWeights = stocGradAscent1(array(trainingSet), trainingLabels, 500)
- # 根据特征向量预测结果
- teststr = '2.000000,1.000000,38.300000,40.000000,24.000000,1.000000,1.000000,3.000000,1.000000,3.000000,3.000000,1.000000,0.000000,0.000000,0.000000,1.000000,1.000000,33.000000,6.700000,0.000000,0.000000'
- currLine = teststr.strip().split(',')
- lineArr = []
- for i in range(len(currLine)):
- lineArr.append(float(currLine[i]))
- res = classifyVector(array(lineArr), trainWeights)
- # 打印预测结果
- reslut = ['死亡','存活']
- print('预测结果是:',int(res))
分类函数代码如下:
- '''分类函数,根据回归系数和特征向量来计算 Sigmoid的值,大于0.5函数返回1,否则返回0'''
- def classifyVector(featuresV, weights):
- prob = sigmoid(sum(featuresV * weights))
- print(prob)
- if prob > 0.9: return 1.0
- else: return 0.0
为了量化回归的效果,需要观察错误率。根据错误率决定是否回退到训练阶段,通过改变迭代的次数和步长的参数来得到更好的回归系数
- '''打开测试集和训练集,并对数据进行格式化处理'''
- def colicTest():
- file_name = './HorseColicTraining.txt'
- trainingSet,trainingLabels = loadDataSet2(file_name)
- # 使用改进后的随机梯度上升算法 求得在此数据集上的最佳回归系数 trainWeights
- trainWeights = stocGradAscent1(array(trainingSet), trainingLabels, 500)
- frTest = open('./HorseColicTest.txt')
- errorCount = 0 ; numTestVec = 0.0
- # 读取 测试数据集 进行测试,计算分类错误的样本条数和最终的错误率
- for line in frTest.readlines():
- numTestVec += 1.0
- currLine = line.strip().split(',')
- lineArr = []
- for i in range(21):
- lineArr.append(float(currLine[i]))
- if int(classifyVector(array(lineArr), trainWeights)) != int(
- currLine[21]):
- errorCount += 1
- errorRate = (float(errorCount) / numTestVec)
- print("逻辑回归算法测试集的错误率为: %f" % errorRate)
- return errorRate
-
-
- # 调用 colicTest() 10次并求结果的平均值
- def multiTest():
- numTests = 10;errorSum = 0.0
- for k in range(numTests):
- errorSum += colicTest()
- print("迭代 %d 次后的平均错误率是: %f" % (numTests, errorSum / float(numTests)))
其运行结果如下:
逻辑回归算法测试集的错误率为: 0.298507
源码请进【机器学习和自然语言QQ群:436303759】文件下载:
本文版权归作者白宁超所有,本文原创,旨在学术和科研使用。文章同步如下:
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。