一,引言 python
前面讲到的基本都是分类问题,分类问题的目标变量是标称型数据,或者离散型数据。而回归的目标变量为连续型,也便是回归对连续型变量作出预测,最直接的办法是依据输入写出一个目标值的计算公式,这样,对于给定的输入,利用该公式能够计算出相应的预测输出。这个公式称为回归方程,而求回归方程显然就是求该方程的回归系数,而一旦有了这些回归系数,再给定输入,就能够将这些回归系数乘以输入值,就获得了预测值。算法
二,线性回归json
线性回归,简单而言,就是将输入项分别乘以一些常量,再将结果加起来获得输出。假设输入数据存放在矩阵x中,而回归系数存放在向量w中,那么对于给定的数据x1,预测结果将会经过y=xTw给出。那么,如何才可以找出最佳的回归系数向量w呢?api
很容易想到使用最小化偏差的w,可是这里的偏差为预测y值和真实y值的差值,使用该偏差的简单累加将会出现正差值和负差值的相互抵消,因此,咱们能够采用平方偏差来进行度量。即:数组
Σ(yi-xiTW)2,i=1,2,3......N(N为样本总数)app
这样,用矩阵表示能够写成(y-Xw)T(y-Xw).由于要求函数的极小值,再对w求导得,xT(y-Xw),则令其等于0,便可获得w的最优解:dom
w*=(XTX)-1XTy函数
须要注意的是,公式中,出现的求逆运算,而对于任意一个矩阵而言,不必定可逆,因此,咱们在实际写代码过程当中,须要事先肯定矩阵是否可逆,不然极可能会形成程序出现严重的错误。学习
有了回归系数的求解公式,咱们就能够写代码来根据样本数据拟合出最佳的回归系数向量w了:测试
form numpy import * #解析文件中的数据为适合机器处理的形式 def loadDataSet(filename): numFeat=len(open(filename).readline().split('\t'))-1 dataMat=[];labelMat=[] fr=open(filename) for line in fr.readlines(): lineArr=[] curLine=line.strip().split('\t') for i in range(numFeat): lineArr.extend(float(curLine[i])) dataMat.append(lineArr) labelMat.append(float(curLine[-1])) return dataMat,labelMat #标准线性回归算法 #ws=(X.T*X).I*(X.T*Y) def standRegres(xArr,yArr): #将列表形式的数据转为numpy矩阵形式 xMat=mat(xArr);yMat=mat(yArr).T #求矩阵的内积 xTx=xMat.T*xMat #numpy线性代数库linalg #调用linalg.det()计算矩阵行列式 #计算矩阵行列式是否为0 if linalg.det(xTx)==0.0: print('This matrix is singular,cannot do inverse') return #若是可逆,根据公式计算回归系数 ws=xTx.I*(xMat.T*yMat) #能够用yHat=xMat*ws计算实际值y的预测值 #返回归系数 return ws
这里,须要说明的是,Numpy提供了一个线性代数的库linalg,其中包含有不少有用的函数,其中就有代码中用到的计算矩阵行列式值得函数linalg.det(),若是行列式值为0,则矩阵不可逆,若是不为0,那么就能够顺利求出回归系数。
此外,咱们知道线性回归的方程的通常形式为:y=wx+b;即存在必定的偏移量b,因而,咱们能够将回归系数和特征向量均增长一个维度,将函数的偏移值b也算做回归系数的一部分,占据回归系数向量中的一个维度,好比w=(w1,w2,...,wN,b)。相应地,将每条数据特征向量的第一个维度设置为1.0,即全部特征向量的第一个维度值均为1.0,这样,最后获得的回归函数形式为y=w[0]+w[1]*x1+...+w[N]*xN.
经过上面的数据咱们能够获得一个拟合的线性回归模型,可是咱们还须要验证该模型的好坏。若是将全部样本的真实值y保存于一个数列,将全部的预测值yHat保存在另一个数列,那么咱们能够经过计算这两个序列的相关系数,来度量预测值与真实值得匹配程度。咱们能够经过NumPy库中提供的corrcoef()函数计算两个序列的相关性。好比以下代码能够获得真实值和预测值序列的相关性,相关性最好为1,最差为0,表示不相关。
#计算预测值序列
yHat=xMat*ws corrcoef(yHat.T,yMat)
三,局部加权线性回归
线性回归一个比较容易出现的问题是有可能出现欠拟合现象,欠拟合显然不能取得最好的预测效果。由于,咱们求的是均方偏差最小的模型,因此能够在估计中引入一些误差,从而下降预测的均方偏差,达到误差和方差的折中,从而找到最佳的模型参数,即回归系数。
解决上述问题的一个方法便是局部加权线性回归(LWLR)。即在算法中,为每个待预测的数据点附近的赋予必定的权重,越靠近预测点的数据点分配的权重越高。这里,咱们采用高斯核函数为预测点附近的数据点分配权重,即:
w(i,i)=exp(|x(i)-x|/-2*k2),其中参数k能够由用户本身定义。
显然,有上面高斯核函数可知,矩阵W是一个只含对角元素的矩阵。这样,回归系数的解的形式变为:w*=(xTWx)-1(xTWy)。对上述代码稍做修改以后获得以下局部加权线性回归函数
#局部加权线性回归 #每一个测试点赋予权重系数 #@testPoint:测试点 #@xArr:样本数据矩阵 #@yArr:样本对应的原始值 #@k:用户定义的参数,决定权重的大小,默认1.0 def lwlr(testPoint,xArr,yArr,k=1.0): #转为矩阵形式 xMat=mat(xArr);yMat=mat(yArr).T #样本数 m=shape(xMat)[0] #初始化权重矩阵为m*m的单位阵 weights=mat(eye((m))) #循环遍历各个样本 for j in range(m): #计算预测点与该样本的误差 diffMat=testPoint-xMat[j,:] #根据误差利用高斯核函数赋予该样本相应的权重 weights[j,j]=exp(diffMat*diffMat.T/(-2.0*k**2)) #将权重矩阵应用到公式中 xTx=xMat.T*(weights*xMat) #计算行列式值是否为0,即肯定是否可逆 if linalg.det(xTx)==0.0 print('This matrix is singular,cannot do inverse') return #根据公式计算回归系数 ws=xTx.I*(xMat.T*(weights*yMat)) #计算测试点的预测值 return testPoint*ws #测试集进行预测 def lwlrTest(testArr,xArr,yArr,k=1.0): #测试集样本数 m=shape(testArr)[0] #测试集预测结果保存在yHat列表中 yHat=zeros(m) #遍历每个测试样本 for i in range(m): #计算预测值 yHat[i]=lwlr(testArr[i],xArr,yArr,k) return yHat
须要说明的是,当为某一预测点附件的数据点分配权重时,由高斯核函数公式可知,与预测点约相近的点分配获得的权重越大,不然权重将以指数级衰减,与预测点足够远的数据点权重接近0,那么这些数据点将不会在该次预测中起做用。固然,用户能够经过本身设定参数k来控制衰减的速度,若是k值较大,衰减速度较慢,这样就会有更多的数据点共同参与决策;不然,若是参数k很是小,那么衰减速度极快,参与预测点决策的数据点就不多。因此,咱们在实验中,应该多选择几组不一样的k值,获得不一样的回归模型,从而找到最优的回归模型对应的k值。
下图从上到下,依次是k取1.0,0.01及0.003状况下,对应的拟合模型。
图中,当k=1.0(上图)时,意味着全部的数据等权重,这样,拟合结果与普通的线性回归拟合结果一致,显然极可能出现欠拟合状况。当k=0.01时(中图),显然获得了最好的拟合效果,由于此时掌握了数据的潜在模式,在预测时剔除了部分不重要数据的影响,增大了重要数据的权重。而k=0.003(下图),能够看出预测时,归入了太多的噪声点,拟合的直线与数据点过于贴近,出现了过拟合的现象。
此外,局部加权线性回归也存在必定的问题,相对于普通的线性回归,因为引入了权重,大大增长了计算量,虽然取得了不错的拟合效果,但也相应地付出了计算量的代价。咱们发现,在k=0.01时,大多的数据点的权重都接近0,因此,若是咱们能避免这些计算,将必定程度上减小程序运行的时间,从而缓解计算量增长带来的问题。
四,示例:预测鲍鱼的年龄
接下来,就利用上述算法对真实的数据进行测试,并计算预测的偏差。鲍鱼的年龄能够经过鲍鱼壳的层数推算获得。在运行代码前,咱们还须要添加计算预测偏差的函数代码:
#计算平方偏差的和 def rssError(yArr,yHatArr): #返回平方偏差和 return ((yArr-yHatArr)**2).sum()
固然,为了获得更好的效果,咱们有必要采起交叉验证的方法,选取多个样本集来进行测试,找出预测偏差最小的回归模型。
五,缩减系数"理解"数据的方法
试想,若是此时数据集样本的特征维度大于样本的数量,此时咱们还能采起上面的线性回归方法求出最佳拟合参数么?显然不可能,由于当样本特征维度大于样本数时,数据矩阵显然是非满秩矩阵,那么对非满秩矩阵求逆运算会出现错误。
为了解决这个问题,科学家提出了岭回归的概念,此外还有一种称为"前向逐步回归"的算法,该算法能够取得很好的效果且计算相对容易。
1 岭回归
简单而言,岭回归便是在矩阵xTx上加入一个λI从而使得矩阵非奇异,进而能对矩阵xTx+λI求逆。其中矩阵I是一个单位矩阵,即对角线上元素皆为1,其余均为0。这样,回归系数的计算公式变为:
w*=(xTx+λI)-1(xTy)
公式中经过引入该惩罚项,从而减小不重要的参数,更好的理解和利用数据。此外,增长了相关约束:Σwi2<=λ,即回归系数向量中每一个参数的平方和不能大于λ,这就避免了当两个或多个特征相关时,可能出现很大的正系数和很大的负系数。
上面的岭回归就是一种缩减方法,经过此方法能够去掉不重要的参数,更好的理解数据的重要性和非重要性,从而更好的预测数据。
在岭回归算法中,经过预测偏差最小化来获得最优的λ值。数据获取以后,将数据随机分红两部分,一部分用于训练模型,另外一部分则用于测试预测偏差。为了找到最优的λ,能够选择多个不一样λ值重复上述测试过程,最终获得一个使预测偏差最小的λ。
#岭回归 #@xMat:样本数据 #@yMat:样本对应的原始值 #@lam:惩罚项系数lamda,默认值为0.2 def ridgeRegres(xMat,yMat,lam=0.2): #计算矩阵内积 xTx=xMat.T*xMat #添加惩罚项,使矩阵xTx变换后可逆 denom=xTx+eye(shape(xMat)[1])*lam #判断行列式值是否为0,肯定是否可逆 if linalg.det(denom)==0.0: print('This matrix is singular,cannot do inverse') return #计算回归系数 ws=denom.I*(xMat.T*yMat) return ws #特征须要标准化处理,使全部特征具备相同重要性 def ridgeTest(xArr,yArr): xMat=mat(xArr);yMat=mat(yArr).T #计算均值 yMean=mean(yMat,0) yMat=yMat-yMean xMeans=mean(xMat,0) #计算各个特征的方差 xVar=var(xMat,0) #特征-均值/方差 xMat=(xMat-xMeans)/xVar #在30个不一样的lamda下进行测试 numTestpts=30 #30次的结果保存在wMat中 wMat=zeros((numTestpts,shape(xMat)[1])) for i in range(numTestpts): #计算对应lamda回归系数,lamda以指数形式变换 ws=ridgeRegres(xMat,yMat,exp(i-10)) wMat[i,:]=ws.T return wMat
须要说明的几点以下:
(1)代码中用到NumPy库中的eye()方法来生成单位矩阵。
(2)代码中仍保留了判断行列式是否为0的代码,缘由是当λ取值为0时,又回到了普通的线性回归,那么矩阵极可能出现不可逆的状况
(3)岭回归中数据须要进行标准化处理,即数据的每一维度特征减去相应的特征均值以后,再除以特征的方差。
(4)这里,选择了30个不一样的λ进行测试,而且这里的λ是按照指数级进行变化,从而能够看出当λ很是小和很是大的状况下对结果形成的影响
下图示出了回归系数与log(λ)之间的关系:
能够看出,当λ很是小时,系数与普通回归同样。而λ很是大时,全部回归系数缩减为0。这样,能够在中间的某处找到使得预测结果最好的λ。
2 逐步前向回归
逐步前向回归是一种贪心算法,即每一步都尽量的减少偏差。从一开始,全部的权重都设为1,而后每一步所作的决策是对某个权重增长或减小一个较小的数值。算法的伪代码为:
数据标准化,使其分布知足均值为0,和方差为1
在每轮的迭代中:
设置当前最小的偏差为正无穷
对每一个特征:
增大或减少:
改变一个系数获得一个新的w
计算新w下的偏差
若是偏差小于当前最小的偏差:设置最小偏差等于当前偏差
将当前的w设置为最优的w
将本次迭代获得的预测偏差最小的w存入矩阵中
返回屡次迭代下的回归系数组成的矩阵
#前向逐步回归 #@eps:每次迭代须要调整的步长 def stageWise(xArr,yArr,eps=0.01,numIt=100): xMat=mat(xArr);yMat=mat(yArr) yMean=mean(yMat,0) yMat=yMat-yMean #将特征标准化处理为均值为0,方差为1 xMat=regularize(xMat) m,n=shape(xMat) #将每次迭代中获得的回归系数存入矩阵 returnMat=zeros((numIt,m)) ws=zeros((n,1));wsTest=ws.copy();wsMat=ws.copy() for i in range(numIt): print ws.T #初始化最小偏差为正无穷 lowestError=inf; for j in range(n): #对每一个特征的系数执行增长和减小eps*sign操做 for sign in [-1,1]: wsTest=ws.copy() wsTest[j]+=eps*sign #变化后计算相应预测值 yTest=xMat*wsTest #保存最小的偏差以及对应的回归系数 rssE=rssError(yMat.A,yTest.A) if rssE<lowestError: lowestError=rssE wsMat=wsTest ws=wsMat.copy() returnMat[i,:]=ws return returnMat
下面看一次迭代的算法执行过程,此时学习步长eps=0.01:
能够看到,上述结果中w1和w6都是0,这代表了这两个维度的特征对预测结果不构成影响,即所谓的不重要特征。此外,在eps=0.01状况下,一段时间后系数就已经饱和并在特定值之间来回震荡,这显然是步长太大的缘故,咱们能够据此调整较小的步长。
前向逐步回归算法也属于缩减算法。它主要的优势是能够帮助咱们更好地理解现有的模型并做出改进。当构建出一个模型时,能够运行该算法找出重要的特征,这样就可能及时中止对不重要特征的收集。一样,在算法测试的过程当中,可使用s折交叉验证的方法,选择偏差最小的模型。
此外,当咱们无论是应用岭回归仍是前向逐步回归等缩减算法时,就相应的为模型增长了误差,与此同时也就减少了模型的方差。而最优的模型每每是在模型误差和方差的折中时得到。不然,当模型方差很大误差很小时模型复杂度很大而出现过拟合现象,而方差很小误差很大时而容易出现欠拟合现象。所以,权衡模型的误差和方差能够作出最好的预测。
六,岭回归应用示例:预测乐高玩具套装的价格
咱们知道乐高玩具是一种拼装类玩具,由不少大小不一样的塑料插件组成。一种乐高玩具套装基本上在几年后就会停产,但乐高收藏者之间仍会在停产后彼此交易。这样,咱们能够拟合一个回归模型,从而对乐高套装进行估价。显然这样作十分有意义。
算法流程:
1 收集数据:用google shopping的api收集数据
2 准备数据:从返回的json数据中抽取价格
3 分析数据:可视化并观察数据
4 训练算法:构建不一样的模型,采用岭回归和普通线性回归训练模型
5 测试算法:使用交叉验证来测试不一样的模型,选择效果最好的模型
1 收集数据:使用Google 购物的API来获取玩具套装的相关信息和价格,能够经过urllib2发送http请求,API将以JSON格式返回须要的产品信息,python的JSON解析模块能够帮助咱们从JSON格式中解析出所须要的数据。收集数据的代码以下:
#收集数据 #添加时间函数库 from time import sleep #添加json库 import json #添加urllib2库 import urllib2 #@retX:样本玩具特征矩阵 #@retY:样本玩具的真实价格 #@setNum:获取样本的数量 #@yr:样本玩具的年份 #@numPce:玩具套装的零件数 #@origPce:原始价格 def searchForSet(retX,retY,setNum,yr,numPce,origPrc): #睡眠十秒 sleep(10) #拼接查询的url字符串 myAPIstr='get from code.google.com' searchURL='https://www.googleapis.com/shopping/search/v1/public/products?\ key=%s&country=US&q=lego+%d&alt=json' %(myAPIstr,setNum) #利用urllib2访问url地址 pg=urllib2.urlopen(searchURL) #利用json打开和解析url得到的数据,数据信息存入字典中 retDict=json.load(pg.read()) #遍历数据的每个条目 for i in range(len(retDict['items'])): try: #得到当前条目 currItem=retDict['items'][i] #当前条目对应的产品为新产品 if currItem['product']['condition']=='new': #标记为新 newFlag=1 else:newFlag=0 #获得当前目录产品的库存列表 listOfInv=currItem['product']['inventories'] #遍历库存中的每个条目 for item in listOfInv: #获得该条目玩具商品的价格 sellingPrice=item['price'] #价格低于原价的50%视为不完整套装 if sellingPrice>origPrc*0.5: print('%d\t%d\t%d\t%f\t%f',%(yr,numPce,newFlag,\ origPrc,sellingPrice)) #将符合条件套装信息做为特征存入数据矩阵 retX.append([yr,numPce,newFlag,origPce]) #将对应套装的出售价格存入矩阵 retY.append(sellingPrice) except:print('problem with item %d',%i) #屡次调用收集数据函数,获取多组不一样年份,不一样价格的数据 def setDataCollect(retX,retY): searchForSet(retX,retY,8288,2006,800,49.99) searchForSet(retX,retY,10030,2002,3096,49.99) searchForSet(retX,retY,10179,2007,5195,499.99) searchForSet(retX,retY,10181,2007,3428,199.99) searchForSet(retX,retY,10189,2008,5922,299.99) searchForSet(retX,retY,10196,2009,3263,249.99)
须要指出的是,咱们在代码中发送http请求,须要从NumPy导入urllib2模块;使用json解析得到的数据时须要导入json模块;同时为避免多个函数同时访问网站,在程序开始时先睡眠必定的时间,用于缓冲,即须要调用time模块的sleep。
此外,因为套装是由多个小插件组成,因此存在插件损失的状况,因此须要过滤掉这样的玩具套装,咱们能够经过必定的过滤条件(通常套装有缺失可能有标识,可经过关键词筛选,或者经过贝叶斯估计估计套装完整性),这里用的是价格来判断,若是当前价格不到原始价格的一半,那么这样的套装必然有缺陷,损坏和缺失等。由于,一旦具备收藏价值的玩具套装停产,必然价格相对上涨,因此出现这种价格反常的状况代表该产品套装存在必定的缺陷,能够将其过滤掉。
2 训练算法:创建模型
上面有了数据,那么咱们就能够开始完成代码进行模型训练了,这里采用岭回归来训练模型,而且采用交叉验证的方法来求出每一个λ对应的测试偏差的均值,最后分析选出预测偏差最小的回归模型。
#训练算法:创建模型 #交叉验证测试岭回归 #@xArr:从网站中得到的玩具套装样本数据 #@yArr:样本对应的出售价格 #@numVal:交叉验证次数 def crossValidation(xArr,yArr,numVal=10): #m,n=shape(xArr) #xArr1=mat(ones((m,n+1))) #xArr1[:,1:n+1]=mat(xArr) #获取样本数 m=len(yArr) indexList=range(m) #将每一个回归系数对应的偏差存入矩阵 errorMat=zeros((numVal,30)) #进行10折交叉验证 for i in range(numVal): trainX=[];trainY=[] testX=[];testY=[] #混洗索引列表 random.shuffle(indexList) #遍历每一个样本 for j in range(m): #数据集90%做为训练集 if j<m*0.9: trainX.append(xArr1[indexList[j]]) trainY.append(yArr[indexList[j]]) #剩余10%做为测试集 else: testX.append(xArr1[indexList[j]]) testY.append(yArr[indexList[j]]) #利用训练集计算岭回归系数 wMat=ridgeRegres(trainX,trainY) #对于每个验证模型的30组回归系数 for k in range(30): #转为矩阵形式 matTestX=mat(testX);matTrainX=mat(trainX) #求训练集特征的均值 meanTrain=mean(matTrainX,0) #计算训练集特征的方差 varTrain=val(matTrainX,0) #岭回归须要对数据特征进行标准化处理 #测试集用与训练集相同的参数进行标准化 matTestX=(matTestX-meanTrain)/varTrain #对每组回归系数计算测试集的预测值 yEst=matTestX*mat(wMat[k,:]).T+mean(trainY) #将原始值和预测值的偏差保存 errorMat[i,k]=rssError(yEst.T.A,array(testY)) #对偏差矩阵中每一个lamda对应的10次交叉验证的偏差结果求均值 meanErrors=mean(errorMat,0) #找到最小的均值偏差 minMean=float(min(meanErrors)) #将均值偏差最小的lamda对应的回归系数做为最佳回归系数 bestWeigths=wMat[nonzero(meanErrors==minMean)] xMat=mat(xArr);yMat=mat(yArr).T meanX=mean(xMat,0);valX=val(xMat,0) #数据标准化还原操做 unReg=bestWeigths/valX print('the best model from Ridge Regression is :\n',unReg) print('with constant term :',-1*sum(multiply(meanX,unReg))+mean(yMat))
一样,这里须要说明的有如下几点:
(1) 这里对于数据集采用随机的方式(random.shffle())选取训练集和测试集,训练集占数据总数的90%,测试集剩余的10%。采起这种方式的缘由是,便于咱们进行屡次交叉验证,获得不一样的训练集和测试集.
(2) 咱们知道岭回归中会选取多个不一样的λ值,来找到预测偏差最小的模型;此外,算法中采用交叉验证的方法,因此对于每个λ对应着多个测试偏差值,因此在分析预测效果最好的λ以前,须要先对每一个λ对应的多个偏差求取均值。
(3) 咱们呢知道岭回归算法须要对训练集数据的每一维特征进行标准化处理,那么为保证结果的准确性,也须要对测试集进行和训练集相同的标准化操做,即测试集数据特征减去训练集该维度特征均值,再除以训练集该维度特征方差
(4) 由于采用岭回归算法时,对数据进行了标准化处理,而标准的回归算法则没有,因此在代码最后咱们仍是须要将数据进行还原,这样便于分析比较两者的真实数据的预测偏差。
由实验结果获得回归方程:
上面模型可能仍是不易理解,再看一下具体的缩减过程当中系数的变化状况:
最后获得的回归系数是通过不一样程度的衰减获得的。好比上图中第一行第四项比第二项系数大5倍,比第一项大57倍。依次来看,第四特征能够看作是最重要特征,在预测时起最主要做用,其次就是第二特征。也便是,特征对应的系数值越大,那么其对预测的决定做用也就越大;若是某一维度系数值为0,则代表该特征在预测结果中不起做用,能够被视为不重要特征。
因此,这种缩减的分析方法仍是比较有用的,由于运算这些算法能够帮助咱们充分理解和挖掘大量数据中的内在规律。当特征数较少时可能效果不够明显,而当特征数至关大时,咱们就能够据此了解特征中哪些特征是关键的,哪些是不重要的,这就为咱们节省很多成本和损耗。
7,总结
(1) 回归与分类的区别,前者预测连续型变量,后者预测离散型变量;回归中求最佳系数的方法经常使用的是最小化偏差的平方和;若是xTx可逆,那么回归算法可使用;能够经过预测值和原始值的相关系数来度量回归方程的好坏
(2) 当特征数大于样本总数时,为解决xTx不可逆的问题,咱们能够经过引入岭回归来保证可以求得回归系数
(3) 另一种缩减算法是,前向逐步回归算法,它是一种贪心算法,每一步经过修改某一维度特征方法来减少预测偏差,最后经过屡次迭代的方法找到最小偏差对应的模型
(4) 缩减法能够看作是对一个模型增长误差的同时减小方差,经过误差方差折中的方法,能够帮助咱们理解模型并进行改进,从而获得更好的预测结果