机器学习实战系列（七）：数值回归与预测

课程的所有数据和代码在我的Github：Machine learning in Action，目前刚开始做，有不对的欢迎指正，也欢迎大家star。除了版本差异，代码里的部分函数以及代码范式也和原书不一样（因为作者的代码实在让人看的别扭，我改过后看起来舒服多了）。在这个系列之后，我还会写一个scikit-learn机器学习系列，因为在实现了源码之后，带大家看看SKT框架如何使用也是非常重要的。

1、线性回归

现有一数据集，其分布如下图所示，

通过观察发现可以通过一个线性方程去拟合这些数据点。可设直线方程为 y=wx. 其中w称为回归系数。那么现在的问题是，如何从一堆x和对应的y中确定w？一个常用的方法就是找出使误差最小的w。这里的误差是指预测y值和真实y值之间的差值，我们采用平方误差，写作：

用矩阵还可以写作：，如果对w求导，得到，令其等于零，解出w为：

注意此处公式包含对矩阵求逆，所以求解时需要先对矩阵是否可逆做出判断。以上求解w的过程也称为“普通最小二乘法”。

几乎任一数据集都可以用上述方法建立模型，只是需要判断模型的好坏，计算预测值yHat和实际值yMat这两个序列的相关系数，可以查看它们的匹配程度。

2、局部加权线性回归

局部加权线性回归给待预测点附近的每个点赋予一定的权重，用于解决线性回归可能出现的欠拟合现象。与ｋＮＮ法类似，这种算法每次预测均需要事先选取出对应的数据子集，然后在这个子集上基于最小均分差来进行普通的回归。该算法解出回归系数的形式如下：

其中w是一个权重矩阵，通常采用核函数来对附近的点赋予权重，最常用的核函数是高斯核，如下：

这样就构建了一个只含对角元素的权重矩阵W并且点x与x(i)越近，w(i,i)将会越大，k值控制衰减速度，且k值越小被选用于训练回归模型的数据集越小。

3、岭回归

如果数据的特征比样本点多（n>m），也就是说输入数据的矩阵x不是满秩矩阵。而非满秩矩阵在求逆时会出错，所以此时不能使用之前的线性回归方法。为解决这个问题，统计学家引入了岭回归的概念。

简单来说，岭回归就是在矩阵xTx上加一个λI从而使得矩阵非奇异，进而能对 xTx+λI 求逆，其中I是一个mxm的单位矩阵。在这种情况下，回归系数的计算公式将变成：

这里通过引入λ来限制了所有w之和，通过引入该惩罚项，能减少不重要的参数，这个技术在统计学中也叫缩减。运行之后得到下图，横轴表示第i组数据，纵轴表示该组数据对应的回归系数值。从程序中可以看出lambda的取值为 exp(i-10) 其中i=0~29。所以结果图的最左边，即λ最小时，可以得到所有系数的原始值（与线性回归一致）；而在右边，系数全部缩减为0；在中间部分的某些值可以取得最好的预测效果。

4、前向逐步回归

前向逐步回归算法属于一种贪心算法，即每一步尽可能减少误差。一开始，所有的权重都设为1，然后每一步所做的决策是对某个权重增加或减少一个很小的值。

该算法伪代码如下所示：

标准线性回归

"""
Created on Sat Jul 14 15:02:40 2018@author: Administrator
"""from numpy import *#解析文件中的数据为适合机器处理的形式
def loadDataSet(filename):numFeat=len(open(filename).readline().split('\t'))-1dataMat=[];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()计算矩阵行列式#计算矩阵行列式是否为0if linalg.det(xTx)==0.0:print('This matrix is singular,cannot do inverse')return #如果可逆，根据公式计算回归系数ws=xTx.I*(xMat.T*yMat)#可以用yHat=xMat*ws计算实际值y的预测值#返回归系数

局部加权线性回归

#局部加权线性回归
def rssError(yArr,yHatArr):#返回平方误差和return ((yArr-yHatArr)**2).sum()def lwlr(testPoint,xArr,yArr,k=1.0):'''#每个测试点赋予权重系数@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

前向逐步回归

#前向逐步回归
#@eps：每次迭代需要调整的步长
def stageWise(xArr,yArr,eps=0.01,numIt=100):xMat=mat(xArr);yMat=mat(yArr)yMean=mean(yMat,0)yMat=yMat-yMean#将特征标准化处理为均值为0，方差为1xMat=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=rssEwsMat=wsTestws=wsMat.copy()returnMat[i,:]=wsreturn returnMat

岭回归

def ridgeRegres(xMat,yMat,lam=0.2):'''#岭回归@xMat:样本数据@yMat：样本对应的原始值@lam：惩罚项系数lamda，默认值为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-yMeanxMeans=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.Treturn wMat

训练模型、交叉验证

#训练算法：建立模型
#交叉验证测试岭回归
#@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).TmeanX=mean(xMat,0);valX=val(xMat,0)#数据标准化还原操作unReg=bestWeigths/valXprint('the best model from Ridge Regression is :\n',unReg)print('with constant term :',-1*sum(multiply(meanX,unReg))+mean(yMat))

预测玩具价格

#收集数据
#添加时间函数库
from time import sleep
#添加json库
import json
#添加urllib2库
import requests
#@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=requests.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=1else: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)