机器学习实战之回归
一，引言
前面讲到的基本都是分类问题，分类问题的目标变量是标称型数据，或者离散型数据。而回归的目标变量为连续型，也即是回归对连续型变量做出预测，最直接的办法是依据输入写出一个目标值的计算公式，这样，对于给定的输入，利用该公式可以计算出相应的预测输出。这个公式称为回归方程，而求回归方程显然就是求该方程的回归系数，而一旦有了这些回归系数，再给定输入，就可以将这些回归系数乘以输入值，就得到了预测值。
二，线性回归
线性回归，简单而言，就是将输入项分别乘以一些常量，再将结果加起来得到输出。假设输入数据存放在矩阵x中，而回归系数存放在向量w中，那么对于给定的数据x1，预测结果将会通过y=xTw给出。那么，如何才能够找出最佳的回归系数向量w呢？
很容易想到使用最小化误差的w，但是这里的误差为预测y值和真实y值的差值，使用该误差的简单累加将会出现正差值和负差值的相互抵消，所以，我们可以采用平方误差来进行度量。即：
Σ(yi-xiTW)2,i=1,2,3......N(N为样本总数)
这样，用矩阵表示可以写成(y-Xw)T(y-Xw).因为要求函数的极小值，再对w求导得，xT(y-Xw)，则令其等于0，即可得到w的最优解：
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) 缩减法可以看做是对一个模型增加偏差的同时减少方差，通过偏差方差折中的方法，可以帮助我们理解模型并进行改进，从而得到更好的预测结果
