logistic回归是机器学习中最经常使用最经典的分类方法之一,有人称之为逻辑回归或者逻辑斯蒂回归。虽然他称为回归模型,可是却处理的是分类问题,这主要是由于它的本质是一个线性模型加上一个映射函数Sigmoid,将线性模型获得的连续结果映射到离散型上。它经常使用于二分类问题,在多分类问题的推广叫softmax。html
本文首先阐述Logistic回归的定义,而后介绍一些最优化算法,其中包括基本的梯度上升法和一个改进的随机梯度上升法,这些最优化算法将用于分类器的训练,最好本文将给出一个Logistic回归的实例,预测一匹病马是否能被治愈。git
在咱们的平常生活中遇到过不少最优化问题,好比如何在最短期内从A点到达B点?如何投入最少工做量却得到最大的效益?如何设计发动机使得油耗最少而功率最大?可见,最优化的做用十分强大,因此此处咱们介绍几个最优化算法,并利用它们训练出一个非线性函数用于分类。github
如今假设有一些数据点,咱们用一条直线对这些点进行拟合(该线称为最佳拟合直线),这个拟合过程就称做回归。利用logistic回归进行分类的主要思想是:根据现有数据对分类边界线创建回归公式,以此进行分类,这里的“回归”一词源于最佳拟合,表示要找到最佳拟合参数集。训练分类器时的作法就是寻找最佳拟合参数,使用的是最优化算法,下面咱们首先介绍一下这个二值型输出分类器的数学原理。算法
优势:计算代码很少,易于理解和实现,计算代价不高,速度快,存储资源低数组
缺点:容易欠拟合,分类精度可能不高网络
适用数据类型:数值型和标称型数据app
咱们想要的函数应该是:能接受全部的输入,而后预测出类型。例如,在两个类的状况下,上述函数输出0或1。该函数称为海维赛德阶跃函数(Heaviside step function),或者直接称为单位阶跃函数。然而,海维赛德阶跃函数的问题在于:该函数在跳跃点上从0瞬间跳跃到1,这个瞬间跳跃过程有时很难处理。幸亏,另外一个函数也有相似的性质(能够输出0或者1),且数学上更易处理,这就是Sigmoid函数。Sigmoid函数具体的计算公式以下:dom
自变量取值为任意实数,值域[0, 1]机器学习
图5-1给出了Sigmoid函数在不一样坐标尺度下的两条曲线图。当x为0时,Sigmoid函数值为0.5。随着x的增大,对应的Sigmoid值将逼近于1;而随着x的减小,Sigmoid值将逼近于0.若是横坐标刻度足够大,Sigmoid函数看起来很像一个阶跃函数。ide
解释Sigmoid函数:将任意的输入映射到了 [0, 1]区间,咱们在线性回归中能够获得一个预测值,再将该值映射到 Sigmoid函数中这样就完成了由值到几率的转换,也就是分类任务。
所以,为了实现Logistic回归分类器,咱们能够在每一个特征上都乘以一个回归系数,而后把全部的结果值相加,将这个总和带入Sigmoid函数中,进而获得一个范围在0~1之间的数值。任何大于0.5的数据被分入1类,小于0.5即被纳入0类,因此,Logistic回归也能够被当作是一种几率估计。
肯定了分类器的函数形式以后,如今的问题变成了:最佳回归系数是多少?如何肯定其大小。
Sigmoid函数的输入记为z,由下面公式获得:
若是采用向量的写法,上述公式能够写成 z = wTx ,它表示将这两个数值向量对应元素相乘,而后所有加起来即获得z值。
其中的向量x是分类器的输入数据,向量w也就是咱们要找到的最佳参数(系数),从而使得分类器尽量的准确,为了寻找该最佳参数,须要用到最优化理论的一些知识。
而后再看看咱们的Logistic回归模型的公式:
这里假设 W>0,Y与X各维度叠加的图形关系,以下图所示(x为了方便取1维):
下面首先学习梯度上升的最优化方法,咱们将学习到如何使用该方法求得数据集的最佳参数,接下来,展现如何绘制梯度上升法产生的决策边界图,该图将梯度上升法的分类效果可视化的呈现出来,最后咱们将学习随机梯度上升算法,以及如何对其进行修改以得到很好地结果。
可能咱们最常听到的是梯度降低算法,它与这里的梯度上升算法是同样的,只是公式中的 加法须要变成减法,梯度上升算法用来求函数的最大值,而梯度降低算法是用来求函数的最小值
梯度上升法基于的思想是:要找到某函数的最大值,最好的方法是沿着该函数的梯度方向探寻,若是梯度记为,则函数 f(x,y) 的梯度由下面式子表示:
这个梯度意味着要沿着x的方向移动,沿着y方向移动
,其中函数f(x,y)必需要在待计算的点上有定义而且可微,一个具体的函数例子见图5-2:
上图中的梯度上升算法沿梯度方向移动了一步,能够看出,梯度算子老是指向函数值增加最快的方向。这里所说的移动方向,而未提到移动量的大小。该量值称为步长,记为。用向量来表示的话,梯度算法的迭代公式以下:
该公式将一直被迭代执行,直至达到某个中止条件为止,好比迭代次数达到某个指定值或算法达到某个能够容许的偏差范围。
基于上面的内容,咱们来看一个Logistic回归分类器的应用例子,从图5-3能够看到咱们采用的数据集。
在LR中,应用极大似然估计法估计模型参数,因为Sigmoid函数的特性,咱们能够作以下的假设:
上式即为在已知样本X和参数θ的状况下。样本X属性正类(y=1)和负类(y=0)的条件几率,将两个公式合并成一个,以下:
假定样本与样本之间相互独立,那么整个样本集生成的几率即为全部样本生成几率的乘积(也就是n个独立样本出现的似然函数以下):
为了简化问题,咱们对整个表达式求对数(即为LR 损失函数):
知足似然函数(θ)的最大的θ值即时咱们须要求解的模型。
那么梯度上升法就像爬坡同样,一点一点逼近极值,而上升这个动做用数学公式表达即为:
其中,α 为步长。
回到Logistic回归问题,咱们一样对函数求偏导。
对这个公式进行分解,先看:
咱们能够看到,对函数求偏导,分解为三部分,而后咱们对这三部分分布求导。
其中:
再由:
可得:
接下来:
最后:
综合三部分即获得:
若是上面链式分解很差理解的话,能够看下面直接求导(结果是同样的):
注意上面是将梯度上升求最大值,转换为梯度降低了,本质没变。
所以梯度迭代公式为:
若是为梯度降低,咱们注意符号的变化,以下:
上图有100个样本点,每一个点包含两个数值型特征:X1和X2,在此数据集上,咱们将经过使用梯度上升法找到最佳回归系数,也就是拟合出Logistic回归模型的最佳参数。
因此咱们的目标:创建分类器,求解出theta参数
设定阈值,根据阈值判断结果
梯度上升法的伪代码以下:
每一个回归系数初始化为1 重复R次: 计算整个数据集的梯度 使用alpha * gradient 更新回归系数的向量 返回回归系数
testSet.txt的文件内容请去个人GitHub下载。地址:https://github.com/LeBron-Jian/MachineLearningNote
下面具体实现梯度上升算法的代码:
#_*_coding:utf-8_*_ from numpy import * # 读取数据 def loadDataSet(filename): ''' 对于testSet.txt,每行前两个值分别是X1和X2,第三个值数据对应的类别标签 并且为了设置方便,该函数还将X0的值设置为1.0 :return: ''' dataMat = [] labelMat = [] fr = open(filename) for line in fr.readlines(): lineArr = line.strip().split() dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])]) labelMat.append(int(lineArr[2])) return dataMat,labelMat def sigmoid(inX): return 1.0/(1+exp(-inX)) def gradAscent(dataMatIn,classLabels): ''' :param dataMatIn: 是一个2维Numpy数组,每列分别表明每一个不一样的特征 每行则表明每一个训练样本。 :param classLabels: 是类别标签,是一个1*100的行向量,为了便于矩阵运算,须要将行向量 转换为列向量,就是矩阵的转置,再将其赋值与labelMat。 :return: ''' dataMatrix = mat(dataMatIn) labelMat = mat(classLabels).transpose() # labelMat = mat(classLabels).T m,n = shape(dataMatrix) # alpha是向目标移动的步长 alpha = 0.001 # 迭代次数 maxCycles = 500 weights = ones((n,1)) for k in range(maxCycles): h = sigmoid(dataMatrix*weights) error = (labelMat-h) weights = weights + alpha*dataMatrix.transpose()*error return weights
注意上面代码(为何将x0的值设置为1呢,主要是咱们要使用矩阵进行运算,能够看下图比较明显)
测试结果以下:
if __name__ == '__main__': filename = 'testSet.txt' dataArr,labelMat = loadDataSet(filename) weights_res = gradAscent(dataArr,labelMat) print(weights_res) ''' [[ 4.12414349] [ 0.48007329] [-0.6168482 ]] '''
上面已经解出了一组回归系数,它肯定了不一样类别数据之间的分割线,那么怎样画出该分割线,从而使得优化的过程便于理解呢?下面代码来解决这个问题。
画出数据集和Logistic回归最佳拟合直线的函数代码:
def plotBestFit(wei): import matplotlib.pyplot as plt weights = wei.getA() dataMat,labelMat = loadDataSet(filename) 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] ax.plot(x,y) plt.xlabel('X1') plt.ylabel('X2') plt.show()
输出的结果和代码以下图所示:
if __name__ == '__main__': filename = 'testSet.txt' dataArr,labelMat = loadDataSet(filename) weights_res = gradAscent(dataArr,labelMat) print(weights_res) plotBestFit(weights_res)
梯度上升算法在500次迭代后获得的Logistic回归最佳拟合直线
这个分类结果至关不错,从图上看,只错分了四个点。可是,尽管例子简单且数据集很小,这个方法却须要大量的计算(300次乘法),所以下一节将对算法稍做改进,从而使它能够用在真实数据集上。
梯度上升算法在每次更新回归系数时都需遍历整个数据集,该方法在处理100个左右的数据集尚可,可是如有数十亿样本和成千上万的特征,那么该方法的计算复杂度就过高了。一种改进方法是一次仅用一个样本点来更新回归系数,该方法称为随机梯度上升算法。因为能够在新样本到来时对分类器进行增量式更新,于是随机梯度上升算法是一个在线学习算法,与“在线学习”相对应,一次处理全部数据被称做是“批处理”。
随机梯度上升算法能够写成以下的伪代码:
全部回归系数初始化为1 对数据集中每一个样本 计算该样本的梯度 使用alpha*gradient 更新回归系数值 返回回归系数值
如下是随机梯度上升算法的实现代码:
# 随机梯度上升算法 def stocGradAscent0(dataMatrix,classLabels): m,n = shape(dataMatrix) alpha = 0.01 weights = ones(n) 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__': filename = 'testSet.txt' dataArr,labelMat = loadDataSet(filename) weights_res = stocGradAscent0(array(dataArr),labelMat) print(weights_res) plotBestFit(weights_res)
图5-5 随机梯度上升算法在上述数据集上的执行结果,最佳拟合直线并不是最佳分类线
能够看出,拟合出来的直线效果还不错,但并不像,可是不像上面那个完美,这里的分类器错分了三分之一的样本。
直接比较结果两个结果是不公平的,后者的结果是在整个数据集上迭代了500次才获得的。一个判断优化算法优劣的可靠方法是看它是否收敛,也就是说参数是否达到了稳定值,是否还会不断地变化?对此,咱们在上面的的随机梯度算法上作了些修改,使其在整个数据集上运行200次,最终绘制的三个回归系数的变化状况以下图所示:
上图展现了随机梯度上升算法在200次迭代过程当中回归系数的变化状况,其中的系数2,也就是图5-5中的X2只通过了50次迭代就达到了稳定值,但系数1和0则须要更屡次的迭代。另外值得注意的是,在大的波动中止后,还有一些小的周期性波动。不难理解,产生这种现象的缘由是存在一些不能正确分类的样本点(数据集并不是现象可分),在每次迭代时会引起系数的剧烈改变。咱们指望算法能避免来回波动,从而收敛到某个值。另外,收敛速度也须要加快。
改进的随机梯度上升算法代码以下:
# 改进的随机梯度上升算法 def stocGradAscent1(dataMatrix,classLabels,numIter=150): m,n = shape(dataMatrix) weights = ones(n) for j in range(numIter): dataIndex = list(range(m)) for i in range(m): alpha = 4/(1.0+j+i)+0.01 randIndex = int(random.uniform(0,len(dataIndex))) h = sigmoid(sum(dataMatrix[randIndex]*weights)) error = classLabels[randIndex] - h weights = weights + alpha *error*dataMatrix[randIndex] del(dataIndex[randIndex]) return weights
上述代码大致上与以前的随机梯度上升算法一致,修改了两处,一处是alpha在每次迭代的时候都会调整,这会环节以前的数据波动或者高频波动。另外,虽然alpha会随着迭代次数不断减小,但永远不会减小到0,。必须这样作的缘由是为了保证在屡次迭代以后新数据仍然具备必定的影响。若是要处理的问题是动态变化的,那么能够适当增长常数项,来确保新的值得到更大的回归系数。另一点值得注意的是,在下降alpha的函数中,alpha每次减小1/(j+1),其中j是迭代次数,i是样本点的下标,这样当j<max(i)的时候,alpha就不是严格降低的,避免参数的严格降低也常见于模拟退火算法等其余优化算法中。程序的第二个改进地方就是经过随机选取样原本更新回归系数,这种方法将减小周期性的波动,而且改进的算法还增长一个迭代次数做为第三个参数,若是该参数没有给定的话,算法将默认迭代150次,若是给定,那么算法将按照新的参数值进行迭代。
与stocGradAscent1()相似,下图显示了每次迭代时各个回归系数的变化状况。
比较5-7和5-6能够看到两点不一样,第一点是,图5-7中的系数没有像5-6里那样出现周期性的波动,这归功于stocGradAscent1()里的样本随机选择机制,第二点是5-7的水平轴比5-6的短了不少,这是因为stocGradAscent1()能够收敛的更快,此次咱们仅仅对数据集作了20次遍历,以前是500次。
下面看看在同一个数据集上的分类效果,将程序运行能够看到:
if __name__ == '__main__': filename = 'testSet.txt' dataArr,labelMat = loadDataSet(filename) weights_res = stocGradAscent1(array(dataArr),labelMat) print(weights_res) plotBestFit(weights_res)
该分割线达到了与GradientAscent()差很少的效果,可是所使用的计算量更少。
默认的迭代次数是150次,可是咱们经过stocGradAscent()的第三个参数来对此进行修改,例如:
weights_res = stocGradAscent1(array(dataArr),labelMat,500)
迄今为止咱们分析了回归系数的变化状况,可是尚未达到目的,,即完成具体的分类任务,下面咱们将使用随机梯度上升算法来解决病马的生死预测问题。,
本次将使用logistic回归来预测患有疝气的马的存活问题,这里的数据包含了368个样本和28个特征。(疝气指的是马胃肠痛的术语,然而这种病不必定源于马 的肠胃问题,其余问题也能够引起疝气)该数据集包含了医院检测马疝气病的一些指标,有的指标比较主观,有的指标难以预测,例如马的疼痛级别
(1)收集数据:使用给定数据文件
(2)准备数据:用Python解析文本文件并填充缺失值
(3)分析数据:可视化并观察数据
(4)训练算法:使用优化算法,找到最佳的系数
(5)测试算法:为了量化回归的结果,须要观察错误率,根据错误率决定是否回退到训练阶段,经过改变迭代的次数和步长等参数来获得更好的回归系数
(6)使用算法:实现一个简单的命令行程序来收集马的症状并输出预测结果
数据中的缺失值是个很是棘手的问题,有不多文献都致力于解决这个问题。那么,数据缺乏究竟带来了什么问题?假设有100个样本和20个特征,这些数据都是机器收集回来的,若机器上的某个传感器损坏致使一个特征无效时该怎么办?此时是否要扔掉整个数据?这种状况下,另外19个特征怎么办?它们是否还可用?答案是确定的。由于有时候数据至关昂贵,扔掉和从新获取都是不可取的,因此必须采起一些方法来解决这个问题。
下面给出了一些可选的作法:
如今,咱们对下一节要用的数据集进行预处理,使其能够顺利地使用分类算法。在预处理阶段须要作两件事:第一,全部的缺失值必须用一个实数值来替换,由于咱们使用的Numpy数据类型不容许包含缺失值。这里选择实数0来替换,由于咱们使用的Numpy数据类型不容许包括缺失值,这里选择实数0来替换全部缺失值,刚好能适应于Logistic回归。这样作的直觉在于,咱们须要的是一个在更新时不会影响系数的值,回归系数的更新公式以下:
weights = weights + alpha *error*dataMatrix[randIndex]
若是dataMatirx的某特征对应值为0,那么该特征的系数将不作更新,即:
weights = weights
另外,因为sigmoid(0) = 0.5 ,即对它结果的预测不具备任何倾向性,所以上述作法也不会对偏差项形成任何影响,基于上述缘由,将缺失值用0代替既能够保留现有数据,也不须要对优化算法进行修改,此外,数据集中的特征值通常不取0,所以在某种意义上说它也知足“特殊值”这个要求。
预处理中作的第二件事,是若是在测试数据集中发现了一条数据的类别标签以及缺失,那么咱们的简单作法是将该条数据丢弃。这是由于类别标签与特征不一样,很难肯定采用某个合适的值来替换,采用Logistic回归进行预处理以后保存两个文件,horseColicTest.txt和horseColicTraining.txt。若是想对于原始数据和预处理以后的数据作个比较。
咱们有一个“干净”可用的数据集和一个不错的优化算法,下面将这些部分融合在一块儿训练出一个分类器,而后利用该分类器来预测病马的生死问题。
horseColicTest.txt的数据 horseColicTraining.txt的数据 请去个人GitHub下载:https://github.com/LeBron-Jian/MachineLearningNote
完整代码:
#-*_coding:utf-8_*_ import math from numpy import * # Sigmoid函数的计算 def sigmoid(inX): return 1.0/(1+exp(-inX)) # 改进的随机梯度上升算法 def stocGradAscent1(dataMatrix,classLabels,numIter=150): m,n = shape(dataMatrix) weights = ones(n) for j in range(numIter): dataIndex = list(range(m)) for i in range(m): alpha = 4/(1.0+j+i)+0.01 randIndex = int(random.uniform(0,len(dataIndex))) h = sigmoid(sum(dataMatrix[randIndex]*weights)) error = classLabels[randIndex] - h weights = weights + alpha *error*dataMatrix[randIndex] del(dataIndex[randIndex]) return weights # Logistic回归分类函数 def classifyVetor(inX,weights): prob = sigmoid(sum(inX*weights)) if prob > 0.5: return 1.0 else: return 0.0 def colicTest(filetrain,filetest): frTrain = open(filetrain) frTest = open(filetest) trainingSet = [] trainingLabeles = [] for line in frTrain.readlines(): currLine = line.strip().split('\t') lineArr = [] for i in range(21): lineArr.append(float(currLine[i])) trainingSet.append(lineArr) trainingLabeles.append(float(currLine[21])) trainWeights = stocGradAscent1(array(trainingSet),trainingLabeles,500) errorCount = 0 numTestVec = 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(classifyVetor(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 multTest(filetrain,filetest): numTests = 6 errorSum = 1.0 for k in range(numTests): errorSum += colicTest(filetrain,filetest) print('after %d iterations the average error rate is %f'%(numTests,errorSum/float(numTests))) if __name__ == '__main__': filetrain = 'horseColicTraining.txt' filetest = 'horseColicTest.txt' multTest(filetrain,filetest)
运行结果:
the error rate of this test is 0.373134 the error rate of this test is 0.358209 the error rate of this test is 0.402985 the error rate of this test is 0.432836 the error rate of this test is 0.462687 the error rate of this test is 0.343284 after 6 iterations the average error rate is 0.562189
从结果来看,6次迭代以后的平均错误率为0.56,事实上,这个结果还不错,由于数据有30%的数据缺失,固然若是调整colicTest()中的迭代次数和stochGradAscent1()中的步长,平均错误率就能够降到20%左右。
这里再补充一个示例,咱们将创建一个逻辑回归模型来预测一个学生是否被大学录取。假设你是一个大学习的管理员,你想根据两次考试的结果来决定每一个申请人的录取机会。你又之前的申请人的历史数据,你能够用它做为逻辑回归的训练集。对于每个培训例子,你有两个考试的申请人的分数和录取决定。为了作到这一点,咱们将创建一个分类模型,根据考试成绩估计入学几率。
数据和代码请参考个人GitHub:https://github.com/LeBron-Jian/MachineLearningNote
首先,拿到数据,并查看部分数据:
import numpy as np import pandas as pd import matplotlib.pyplot as plt path = 'LogiReg_data.txt' pdData = pd.read_csv(path, header=None, names=['Exam 1', 'Exam 2', 'Admitted']) print(pdData.head(10))
部分数据以下:
拿到数据后,咱们将其按照结果分为录取和不被录取的集合,并展现:
positive = pdData[pdData['Admitted'] == 1] negative = pdData[pdData['Admitted'] == 0] fig, ax = plt.subplots(figsize=(5, 5)) ax.scatter(positive['Exam 1'], positive['Exam 2'], s=30, c='b', marker='o', label='Admitted') ax.scatter(negative['Exam 1'], negative['Exam 2'], s=30, c='r', marker='x', label='Not Admitted') ax.legend() ax.set_xlabel('Exam 1 Score') ax.set_ylabel('Exam 2 Score')
图示以下:
下面使用逻辑回归来创建分类器,咱们须要求解出三个参数 theta0, theta1, theta2。 而后设定阈值,根据阈值判断录取结果。
要完成的模块函数:
首先,完成Sigmoid函数:
def sigmoid(z): return 1 / (1 + np.exp(-z)) # 若是想查看Sigmoid函数,能够运行下面代码 nums = np.arange(-10, 10, step=1) fig, ax = plt.subplots(figsize=(12, 4)) ax.plot(nums, sigmoid(nums), 'r')
而后对数据集进行划分,分为训练数据和标签,其中咱们添加了一列,而且初始化了theta值
# 增长一行1 在第0列 pdData.insert(0, 'Ones', 1) print(pdData.head()) # 设置X和Y 即训练数据和训练变量 orig_data = pdData.values cols = orig_data.shape[1] # (100, 4) X = orig_data[:, 0:cols-1] # 等价于 np.matrix(X.values) Y = orig_data[:, cols-1:cols] # 等价于 np.matrox(data.iloc[:, 3:4].value # 传递numpu矩阵 并初始化参数矩阵 theta theta = np.zeros([1, 3]) # [[0. 0. 0.]]
咱们能够看一下添加了一列后,整体的数据:
其余的数据都可以本身打印。
下面是损失函数:
代码整理以下:
def cost(X, y, theta): left = np.multiply(-y, np.log(model(X, theta))) right = np.multiply(1-y, np.log(1 - model(X, theta))) return np.sum(left - right) / (len(X))
下面是计算梯度,梯度函数以下(计算梯度就是咱们对theta求偏导):
代码可写成下面方式:
def gradient(X, y, theta): grad = np.zeros(theta.shape) error = (model(X, theta)-y).ravel() for j in range(len(theta.ravel())): # for each parmeter term = np.multiply(error, X[:, j]) grad[0, j] = np.sum(term) / len(X) return grad
下面比较三种不一样梯度降低方法,三种梯度降低分别是 批量梯度降低,随机梯度降低,小批量梯度降低。以下图:
咱们的方法为1,根据迭代次数中止 2,根据损失函数中止 3,根据梯度降低中止
(注意:咱们每次进行迭代更新的时候,都会对数据进行洗牌,一把状况数据都会有某种规律)
# 比较三种不一样梯度降低方法 STOP_ITER = 0 STOP_COST = 1 STOP_GRAD = 2 def stopCriterion(type, value, threshold): # 设定三种不一样的中止策略 if type == STOP_ITER: return value > threshold elif type == STOP_COST: return abs(value[-1]-value[-2]) < threshold elif type == STOP_GRAD: return np.linalg.norm(value) < threshold # 洗牌 def shuffleData(data): np.random.shuffle(data) cols = data.shape[1] X = data[:, 0:cols-1] y = data[:, cols-1:] return X, y def descent(data, theta, batchSize, stopType, thresh, alpha): # 梯度降低求解 init_time = time.time() i = 0 # 迭代次数 k = 0 # batch X, y = shuffleData(data) grad = np.zeros(theta.shape) # 计算的梯度 costs = [cost(X, y, theta)] # 损失值 while True: grad = gradient(X[k:k+batchSize], y[k:k+batchSize], theta) k += batchSize # 取batch数量个数据 if k>=n: k = 0 X, y = shuffleData(data) # 从新洗牌 theta = theta = alpha*grad # 参数更新 costs.append(cost(X, y, theta)) # 计算新的损失 i += 1 if stopType == STOP_ITER: value = i elif stopType == STOP_COST: value = costs elif stopType == STOP_GRAD: value = grad if stopCriterion(stopType, value, thresh): break return theta, i-1, costs, grad, time.time()-init_time
下面咱们定义一个功能性的函数,此函数的做用就是执行一次梯度降低,并完成更新,选择中止策略。其中主要代码就第一句。。。
def runExpe(data, theta, batchSize, stopType, thresh, alpha): theta, iter, costs, grad, dur = descent(data, theta, batchSize, stopType, thresh, alpha) name = 'Original' if (data[:, 1] > 2).sum() > 1 else 'Scaled' name += 'data - learning rate: {} -'.format(alpha) if batchSize == n: strDescType = 'Gradient' elif batchSize == 1: strDescType = 'Stochastic' else: strDescType = 'Min-batch({})'.format(batchSize) name += strDescType + 'descent - Stop: ' if stopType == STOP_ITER: strStop = "{} iterations".format(thresh) elif stopType == STOP_COST: strStop = "costs change < {}".format(thresh) else: strStop = 'gradient norm < {}'.format(thresh) name += strStop print("***{}\nTheta: {} - Iter:{}-Last cost: {:03.2f}--Duration: {:03.2f}s".format( name, theta, iter, costs[-1], dur )) fig, ax = plt.subplots(figsize=(12, 4)) ax.plot(np.arange(len(costs)), costs, 'r') ax.set_xlabel('Iterations') ax.set_ylabel('Cost') ax.set_title(name.upper() + '- Error vs. Iteration') return theta
下面按照不一样的中止策略画图展现:
1,设定迭代次数 5000
2,设定阈值1e-6,差很少须要 110000次迭代
3,设置阈值 0.05,差很少须要 40000次迭代
if __name__ == '__main__': # 设定迭代次数 选择的梯度降低方法是基于全部样本的 n = 100 # runExpe(orig_data, theta, n, STOP_ITER, thresh=5000, alpha=0.000001) # 根据损失值中止 设定阈值 1e-6 差很少须要110 000 次迭代 # runExpe(orig_data, theta, n, STOP_COST, thresh=0.000001, alpha=0.001) # 根据梯度变换中止 设定阈值 0.05,差很少须要 40000次迭代 # runExpe(orig_data, theta, n, STOP_GRAD, thresh=0.05, alpha=0.001) # 对比不一样的梯度降低方法 stochastic descent 这个不收敛 # runExpe(orig_data, theta, 1, STOP_ITER, thresh=5000, alpha=0.001) # 咱们将学习率调小,提升thresh再试试,发现收敛了 # runExpe(orig_data, theta, 1, STOP_ITER, thresh=15000, alpha=0.000001) # Mine-batch descent 这个也不收敛 ,咱们将alpha调小 0.000001 发现收敛了 runExpe(orig_data, theta, 16, STOP_ITER, thresh=15000, alpha=0.001)
浮动仍然比较大,咱们来尝试下对数据进行标准化,将数据按其属下(按列进行)减去其均值,而后除以其方差。最后获得的结果是,对每一个属性/每列来讲全部数据都汇集在 0 附加,方差值为1。
scaled_data = orig_data.copy() scaled_data[:, 1:3] = pp.scale(orig_data[:, 1:3]) runExpe(scaled_data, theta, n, STOP_ITER, thresh=5000, alpha=0.01)
下面两幅图,咱们分别将学习率设置为 0.001 0.01
· 咱们发现,这里已经降到了0.3如下,因此说作数据预处理很是重要。
scaled_data = orig_data.copy() scaled_data[:, 1:3] = pp.scale(orig_data[:, 1:3]) # runExpe(scaled_data, theta, n, STOP_ITER, thresh=5000, alpha=0.01) runExpe(scaled_data, theta, n, STOP_GRAD, thresh=0.02, alpha=0.01)
咱们修改一下,根据梯度变换中止,发现更多的迭代次数会使得损失降低的更多
runExpe(scaled_data, theta, 1, STOP_GRAD, thresh=0.002/5, alpha=0.01)
随机梯度降低更快,可是咱们须要迭代的次数也须要更多,因此仍是用 batch的比较合适。
精度
# 设置阈值 def predict(x, theta): return [1 if x >= 0.5 else 0 for x in model(X, theta)] if __name__ == '__main__': # 设定迭代次数 选择的梯度降低方法是基于全部样本的 n = 100 scaled_data = orig_data.copy() scaled_X = scaled_data[:, :3] y = scaled_data[:, 3] predictions = predict(scaled_X, theta) correct = [1 if ((a == 1 and b == 1) or (a == 0 and b == 0)) else 0 for a, b in zip(predictions, y)] accuracy = (sum(map(int, correct)) % len(correct)) print('accuracy = {0}%'.format(accuracy)) print(len(correct))
不知道为何,我求出来的精度不高。。。。百分之60。
Logistic回归的目的是寻找一个非线性函数Sigmoid的最佳拟合参数,求解过程能够由最优化算法来完成。在最优化算法中,最经常使用的就是梯度上升算法,而梯度上升算法又能够简化为随机梯度上升算法。
随机梯度上升算法与梯度上升算法的效果至关,可是占用更少的计算资源。此外,随机梯度上升是一个在线算法,它能够在新数据到来时就完成参数更新,不须要从新读取整个数据集来进行批处理运算。
机器学习的一个重要问题就是如何处理缺失数据,这个问题没有标准答案,取决于实际应用中的需求。
前面也说了Logistic回归模型主要用于二分类,那么下面说一下多分类问题中的推广——softmax回归。
softmax与Logistic回归的主要区别就是,Logistic处理二分类问题,只有一组权重参数θ,而softmax处理多分类问题,若是有k个类别,那么softmax就有k组权值参数。每组权值对应一种分类,经过k组权值求解出样本数据对应每一个类别的几率,最后取几率最大的类别做为该数据的分类结果,它的几率函数为:
softmax函数常常用于神经网络的最后一层,用于对神经网络已经处理好的特征进行分类。
总结来讲:逻辑回归真的真的很好用!
下面让咱们看一下Sklearn的Logistic回归分类器
英文的Sklearn文档地址:请点击我
sklearn.linear_model模块提供了不少模型供咱们使用,好比Logistic回归、Lasso回归、贝叶斯脊回归等,可见须要学习的东西还有不少不少。本次,咱们使用LogisticRegressioin。
让咱们先看一下LogisticRegression这个函数,一共有14个参数
参数说明以下:
penalty:惩罚项,str类型,可选参数为l1和l2,默认为l2。用于指定惩罚项中使用的规范。newton-cg、sag和lbfgs求解算法只支持L2规范。L1G规范假设的是模型的参数知足拉普拉斯分布,L2假设的模型参数知足高斯分布,所谓的范式就是加上对参数的约束,使得模型更不会过拟合(overfit),可是若是要说是否是加了约束就会好,这个没有人能回答,只能说,加约束的状况下,理论上应该能够得到泛化能力更强的结果。
dual:对偶或原始方法,bool类型,默认为False。对偶方法只用在求解线性多核(liblinear)的L2惩罚项上。当样本数量>样本特征的时候,dual一般设置为False。
tol:中止求解的标准,float类型,默认为1e-4。就是求解到多少的时候,中止,认为已经求出最优解。
c:正则化系数λ的倒数,float类型,默认为1.0。必须是正浮点型数。像SVM同样,越小的数值表示越强的正则化。
fit_intercept:是否存在截距或误差,bool类型,默认为True。若是使用中心化的数据,能够考虑设置为False,不考虑截距。注意这里是考虑,通常仍是要考虑截距
intercept_scaling:仅在正则化项为”liblinear”,且fit_intercept设置为True时有用。float类型,默认为1。
class_weight:用于标示分类模型中各类类型的权重,能够是一个字典或者balanced字符串,默认为不输入,也就是不考虑权重,即为None。若是选择输入的话,能够选择balanced让类库本身计算类型权重,或者本身输入各个类型的权重。举个例子,好比对于0,1的二元模型,咱们能够定义class_weight={0:0.9,1:0.1},这样类型0的权重为90%,而类型1的权重为10%。若是class_weight选择balanced,那么类库会根据训练样本量来计算权重。某种类型样本量越多,则权重越低,样本量越少,则权重越高。当class_weight为balanced时,类权重计算方法以下:n_samples / (n_classes * np.bincount(y))。n_samples为样本数,n_classes为类别数量,np.bincount(y)会输出每一个类的样本数,例如y=[1,0,0,1,1],则np.bincount(y)=[2,3]。
那么class_weight有什么做用呢? 在分类模型中,咱们常常会遇到两类问题:
random_state:随机数种子,int类型,可选参数,默认为无,仅在正则化优化算法为sag,liblinear时有用。
solver:优化算法选择参数,只有五个可选参数,即newton-cg,lbfgs,liblinear,sag,saga。默认为liblinear。solver参数决定了咱们对逻辑回归损失函数的优化方法,有四种算法能够选择,分别是:
总结:
max_iter:算法收敛最大迭代次数,int类型,默认为10。仅在正则化优化算法为newton-cg, sag和lbfgs才有用,算法收敛的最大迭代次数。
multi_class:分类方式选择参数,str类型,可选参数为ovr和multinomial,默认为ovr。ovr即前面提到的one-vs-rest(OvR),而multinomial即前面提到的many-vs-many(MvM)。若是是二元逻辑回归,ovr和multinomial并无任何区别,区别主要在多元逻辑回归上。
verbose:日志冗长度,int类型。默认为0。就是不输出训练过程,1的时候偶尔输出结果,大于1,对于每一个子模型都输出。
warm_start:热启动参数,bool类型。默认为False。若是为True,则下一次训练是以追加树的形式进行(从新使用上一次的调用做为初始化)。
n_jobs:并行数。int类型,默认为1。1的时候,用CPU的一个内核运行程序,2的时候,用CPU的2个内核运行程序。为-1的时候,用全部CPU的内核运行程序。
除此以外,LogisticRegression也有一些方法供咱们使用:
有一些方法和MultinomialNB的方法都是相似的,所以再也不累述。
了解了基本知识点,咱们就能够编写Sklearn分类器的代码了,代码以下:
#_*_ codingLutf-8_*_ from sklearn.linear_model import LogisticRegression def colicSklearn(filetrain,filetest): frTrain = open(filetrain) frTest = open(filetest) trainingSet = [] trainingLabels = [] testSet = [] testLabels = [] for line in frTrain.readlines(): currLine = line.strip().split('\t') lineArr = [] for i in range(len(currLine)-1): lineArr.append(float(currLine[i])) trainingSet.append(lineArr) trainingLabels.append(float(currLine[-1])) for line in frTest.readlines(): currLine = line.strip().split('\t') lineArr = [] for i in range(len(currLine)-1): lineArr.append(float(currLine[i])) testSet.append(lineArr) testLabels.append(float(currLine[-1])) classifier = LogisticRegression(solver='liblinear',max_iter=20).fit(trainingSet,trainingLabels) test_accurcy = classifier.score(testSet,testLabels)*100 print("正确率为%s%%"%test_accurcy) if __name__ == '__main__': filetrain = 'horseColicTraining.txt' filetest = 'horseColicTest.txt' colicSklearn(filetrain,filetest)
执行结果以下:
正确率为73.13432835820896%
能够看到,正确率又搞了,更改solver参数,好比设置为sag,使用随机平均梯度降低算法,看一看效果,你会发现:正确率高了,可是发出了警告,发出警告是由于算法尚未收敛,更改迭代次数便可。
当咱们迭代到5000的时候,就不会报错了,因此说,对于小数据集,sag算法须要迭代上千次才收敛,而liblinear只须要不到10次。
因此,咱们须要根据数据集状况,选择最优化算法。
首先,LR将线性模型利用Sigmoid函数进一步作了非线性映射。将分类超平面两侧的正负样本点经过压缩函数转化成了以 0.5 为分类的两类:类别0 和类别1。
在上面咱们讲述了线性边界与LR分布函数(即Sigmoid函数)的映射对应关系;一样,对于非线性断定边界Sigmoid函数也使用,以下:
而后,再去考虑如何去得到模型参数,就是断定边界的参数怎么得到,这里是利用MLE进行求解的,具体求解过程能够百度。
代码以下:
#_*_coding:utf-8_*_ import numpy as np import matplotlib.pyplot as plt def sigmoid(z): return 1.0/(1.0 + np.exp(-z)) z = np.arange(-7, 7, 0.1) print(z) phi_z = sigmoid(z) def cost_1(z): return -np.log(sigmoid(z)) def cost_0(z): return -np.log(1-sigmoid(z)) z = np.arange(-10, 10, 0.1) phi_z = sigmoid(z) c1 = [cost_1(x) for x in z] plt.plot(phi_z, c1, label='J(w) if y=1') c0 = [cost_0(x) for x in z] plt.plot(phi_z, c0, linestyle='--', label='J(w) if y=0') plt.ylim(0.0, 5.1) plt.xlim([0, 1]) plt.xlabel('$\phi$(z)') plt.ylabel('J(w)') plt.legend(loc='best') plt.tight_layout() # plt.savefig('./figures/log_cost.png', dpi=300) plt.show()
结果以下:
最后,考虑LR进行三分类(多分类)时,是特征的线性组合和Sigmoid函数符合的函数进行几率计算和分类的。
代码:
from IPython.display import Image # Added version check for recent scikit-learn 0.18 checks from distutils.version import LooseVersion as Version from sklearn import __version__ as sklearn_version from sklearn import datasets import numpy as np iris = datasets.load_iris() # http://scikit-learn.org/stable/auto_examples/datasets/plot_iris_dataset.html X = iris.data[:, [2, 3]] print(X.shape) y = iris.target # 取species列,类别 if Version(sklearn_version) < '0.18': from sklearn.cross_validation import train_test_split else: from sklearn.model_selection import train_test_split # train_test_split方法分割数据集 X_train, X_test, y_train, y_test = train_test_split( X, y, test_size=0.3, random_state=0) from sklearn.preprocessing import StandardScaler sc = StandardScaler() # 初始化一个对象sc去对数据集做变换 sc.fit(X_train) # 用对象去拟合数据集X_train,而且存下来拟合参数 X_train_std = sc.transform(X_train) X_test_std = sc.transform(X_test) from sklearn.linear_model import LogisticRegression def sigmoid(z): return 1.0 / (1.0 + np.exp(-z)) lr = LogisticRegression(C=1000.0, random_state=0) # http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html#sklearn.linear_model.LogisticRegression lr.fit(X_train_std, y_train) # 计算该预测实例点属于各种的几率 lr.predict_proba(X_test_std[0, :].reshape(1, -1)) # Output:array([[ 2.05743774e-11, 6.31620264e-02, 9.36837974e-01]]) # 验证predict_proba的做用 c = lr.predict_proba(X_test_std[0, :].reshape(1, -1)) # c[0, 0] + c[0, 1] + c[0, 2] # Output:0.99999999999999989 # 查看lr模型的特征系数 lr = LogisticRegression(C=1000.0, random_state=0) # http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html#sklearn.linear_model.LogisticRegression lr.fit(X_train_std, y_train) print(lr.coef_) # Output:[[-7.34015187 -6.64685581] # [ 2.54373335 -2.3421979 ] # [ 9.46617627 6.44380858]] # 验证predict_proba工做原理 Zz = np.dot(lr.coef_, X_test_std[0, :].T) + lr.intercept_ np.array(sigmoid(Zz)) / sum(np.array(sigmoid(Zz))) # Output:array([ 2.05743774e-11, 6.31620264e-02, 9.36837974e-01]) # 此结果就是预测实例点各种的几率
对鸢尾花进行分类代码:
from sklearn import datasets from numpy import * from sklearn.linear_model import LogisticRegression from sklearn.model_selection import train_test_split def colicSklearn(): iris = datasets.load_iris() X = iris.data Y = iris.target trainingSet,testSet,trainingLabels,testLabels = train_test_split(X,Y,test_size=0.25,random_state=40) classifier = LogisticRegression(solver='sag', max_iter=5000).fit(trainingSet, trainingLabels) test_accurcy = classifier.score(testSet, testLabels) * 100 print("正确率为%s%%" % test_accurcy) if __name__ == '__main__': colicSklearn()
结果:
正确率为100.0%
代码:
#_*_coding:utf-8_*_ ''' 下面这个例子,从数据产生,到数据提取,数据标准化 模型训练和评估来讲明各个API的调用过程 ''' print(__doc__) from sklearn.preprocessing import StandardScaler from sklearn.linear_model import LinearRegression from sklearn.model_selection import train_test_split import matplotlib.pyplot as plt import numpy as np import pandas as pd # import matplotlib as mpl import pylab as mpl # 设置字符集,防止中文乱码 mpl.rcParams['font.sans-serif'] = [u'simHei'] mpl.rcParams['axes.unicode_minus'] = False # 定义目标函数经过改函数产生对应的y # y = 1 *x[0] + 2 *x[1] + ... (n+1) *x[n] def l_model(x): params = np.arange(1, x.shape[-1] +1) y = np.sum(params *x) + np.random.randn(1) * 0.1 return y # 定义数据集 x = pd.DataFrame(np.random.rand(500,6)) # print(x) y = x.apply(lambda x_rows:pd.Series(l_model(x_rows)),axis=1) # print(y) # 划分数据集 x_train, x_test, y_train, y_test = train_test_split(x,y,test_size=0.3,random_state=2) # 数据标准化 ss = StandardScaler() x_train_s = ss.fit_transform(x_train) x_test_s = ss.fit_transform(x_test) # 输出下元数据的标准差和平均数 print(ss.scale_) print(ss.mean_) # 训练模型 lr = LinearRegression() lr.fit(x_train_s , y_train) # 训练后的输入端模型系数,若是label有两个,即y值有两列,那么是一个2D的array print(lr.coef_) # 截距 print(lr.intercept_) # 用模型预测,并计算得分 y_predict = lr.predict(x_test_s) test_accuracy = lr.score(x_test_s, y_test) print("正确率为%s%%" % test_accuracy) # 预测值和实际值画图比较 t = np.arange(len(x_test_s)) # 建一个画布,facecolor是背景色 plt.figure(facecolor='W') plt.plot(t, y_test, 'r-', linewidth= 2, label = '真实值') plt.plot(t, y_predict, 'b-', linewidth= 1, label = '预测值') # 显示图例,设置图例的位置 plt.legend(loc= 'upper left') plt.title("线性回归预测真实值之间的关系", fontsize = 20) # 加网格 plt.grid(b = True) plt.show()
结果展现:
下面这个例子,从数据产生,到数据提取,数据标准化 模型训练和评估来讲明各个API的调用过程 [0.28299035 0.3138177 0.29522932 0.30125841 0.28402977 0.28989605] [0.53282546 0.51045318 0.48925819 0.5202345 0.4989613 0.54335117] [[0.29202621 0.5787407 0.87280326 1.12786483 1.50043119 1.73271735]] [10.38441274] 正确率为0.9706613427158242%
传送门:请点击我
若是点击有误:https://github.com/LeBron-Jian/MachineLearningNote
参考:https://www.cnblogs.com/bonelee/p/7253508.html
https://blog.csdn.net/c406495762/article/details/77851973