机器学习实战-预测数值型数据：回归

1，用线性回归找到最佳拟合直线

回归的目的是预测数值型的目标值。

回归方程（regression equation）主要是求回归系数，一旦有了回归系数，在给定输入，做预测就是用回归系数乘以输入值，在将结果全部加在一起，就得到了预测值。（因为回归系数是一个向量，输入也是向量，这些运算就是求出二者的内积）

回归的一般方法：

1）收集数据：采用任意方法收集数据

2）准备数据：回归需要数值型数据，标称型数据将被转换为二值型数据

3）分析数据，二维图

4）训练算法：找到回归系数

5）测试算法：使用 $R^{2}$ 或者预测值和数据的拟合度，来分析模型的效果

6）使用算法：使用回归，在给定的时候预测出一个数值

求回归方程：

假定输入数据存放在矩阵X中，而回归系数存放在向量w中。那么对于给定的数据 $x_{1}$ ,预测结果将会通过 $Y_{1}=X_{1}w$ 给出。通过x,y数据找到回归系数向量w,一个常用的方法是找出误差最小的w。这里的误差是指预测值和真实值之间的差值。避免正差值和负差值之间相互抵消，因此采用平方误差。

平方误差可以写作

$\sum_{i=1}^{m}(y_{i}-x_{i}^Tw)^2$

用矩阵表示还可以写作 $\hat{w}=(Y-Xw)^{T}(Y-Xw)$ ,如果对w求导，得到 $X^{T}(Y-Xw)$ ，令其等于零，解出w

$\hat{w}=(X^TX)^{-1}X^Ty$

w上方的小标记表示，这事当前可以估计出的w的最优解，从现有数据上估计出的w可能并不是数据中的真实值，所以这里使用了一个"帽"符号来表示是w的一个最佳估计。注意：公式中包含 $(X^TX)^{-1}$ ，也就是需要对矩阵求逆，代码中需要判断矩阵的逆到底存不存在

另外还可以使用NumPy库中的矩阵方法，该方法叫OLS,意思是“普通最小二乘法”（ordinary least squares）

from numpy import *#加载数据
def loadDataSet(fileName):      #general function to parse tab -delimited floatsnumFeat = len(open(fileName).readline().split('\t')) - 1 #get number of fields dataMat = []; labelMat = []fr = open(fileName)for line in fr.readlines():lineArr =[]curLine = line.strip().split('\t')for i in range(numFeat):lineArr.append(float(curLine[i]))dataMat.append(lineArr)labelMat.append(float(curLine[-1]))return dataMat,labelMat#返回回归系数
def standRegres(xArr,yArr):xMat = mat(xArr); yMat = mat(yArr).TxTx = xMat.T*xMatif linalg.det(xTx) == 0.0:# det判断行列式print ("This matrix is singular, cannot do inverse")returnws = xTx.I * (xMat.T*yMat)return wsimport matplotlib.pyplot as plt
xmat=mat(xarr)
ymat=mat(yarr)
yhat=xmat*ws
fig=plt.figure()
ax=fig.add_subplot(111)
#a是个矩阵或者数组，a.flatten()就是把a降到一维，默认是按横的方向降
# 此时的a是个矩阵，降维后还是个矩阵，矩阵.A（等效于矩阵.getA()）变成了数组，A[0]就是数组里的第一个元素
ax.scatter(xmat[:,1].flatten().A[0],ymat.T[:,0].flatten().A[0])
xCopy=xmat.copy()
xCopy.sort(0)
yhat=xCopy*ws
ax.plot(xCopy[:,1],yhat)
plt.savefig("2.png")
plt.show()

几乎任意数据集都可以用上述方法建立模型，但是判断模型的好坏，可以通过计算着两个序列的相关系数来，计算预测值yHat序列和真实值y序列的匹配程度。

yhatT=yhat.T
coef=corrcoef(yhatT,ymat)
print(coef)[[1.         0.13653777][0.13653777 1.        ]]

2,局部加权线性回归

线性回归的一个问题是有可能出现欠拟合现象，因为它求的是具有最小均方误差的无偏估计。

局部加权线性回归（Locally weighted linear regression,LWLR）。在该算法中，给每个待测点附近的每个点赋予一定的权重；然后基于最小均方差来进行普通的回归。

$\hat{w}=(X^TWX)^{-1}X^TWy$

其中w是一个矩阵，用来给每个数据点赋予权重

LWLR使用“核”（与支持向量机中的核类似）来对附近的点赋予更高的权重。核的类型可以自由选择，最常用的核是高斯核，高斯核对应的权重如下

$w(i,i)=exp(\frac{\left | x^{(i)}-x \right |}{-2k^2})$

这样就构建了一个只含有对角元素的权重矩阵，并且点x与x(i)越近，w(i,i)将会越大()。参数k,决定了对附近的点赋予多大的权重，这也是LWLR时唯一需要考虑的参数，

import matplotlib.pyplot as pltdef lwlr(testPoint,xArr,yArr,k=1.0):xMat = mat(xArr); yMat = mat(yArr).Tm = shape(xMat)[0]weights = mat(eye((m)))for j in range(m):                      #next 2 lines create weights matrixdiffMat = testPoint - xMat[j,:]     #weights[j,j] = exp(diffMat*diffMat.T/(-2.0*k**2))xTx = xMat.T * (weights * xMat)if linalg.det(xTx) == 0.0:print("This matrix is singular, cannot do inverse")returnws = xTx.I * (xMat.T * (weights * yMat))return testPoint * wsdef lwlrTest(testArr,xArr,yArr,k=1.0):  #loops over all the data points and applies lwlr to each onem = shape(testArr)[0]yHat = zeros(m)for i in range(m):yHat[i] = lwlr(testArr[i],xArr,yArr,k)return yHatxarr,yarr=loadDataSet('ex0.txt')
xmat=mat(xarr)
ymat=mat(yarr)
yhat=lwlrTest(xarr,xarr,yarr,0.01)
srcInd=xmat[:,1].argsort(0)
xsort=xmat[srcInd][:,0,:]
print(xsort)
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(xsort[:,1],yhat[srcInd])
ax.scatter(xmat[:,1].flatten().A[0],mat(yarr).T.flatten().A[0],s=2,c='red')
plt.show()

局部加权线性回归也存在一个问题，即增加了计算量，因为他对每个点做预测是都必须使用整个数据集

3，预测鲍鱼的年龄

使用较小的核将得到较低的误差，但是不能所有数据集都使用最小的核，因为使用最小的核将造成过拟合，对新数据不一定能达到最好的预测效果。

4，缩减发来“理解”数据，

如果数据的特征比样本点还多，就不能使用之前的方法来做预测。因为 $(X^TX)^{-1}$ 计算的时候回报错。特征比样本点还多（n>m）,也就是说输入数据的矩阵X不是满秩矩阵。非满秩矩阵在求逆向时会出现问题。一般采用岭回归和lasso回归

4.1岭回归（ridge regression）

岭回归就是在矩阵 $X^TX$ 上加一个 $\lambda$ I 求逆，其中矩阵I是mxm的单位矩阵， $\lambda$ 是用户定义的一个数值，回归系数变为

$\hat{w}=(X^TX+\lambda I)^{-1}X^Ty$

岭回归最先用来处理特征数多于样本数的情况，现在也用于在估计中加入偏差，从而得到更好的估计。通过限制 $\lambda$ 来限制所有w之和，通过引入惩罚性，能够减少不重要的参数，能更好的理解数据，这个技术叫做缩减(shrinkage)

通过预测误差最小化得到 $\lambda$ 。通过选取不同的 $\lambda$ 来重复上述测试过程，最终得到一个是测试误差最小的 $\lambda$

import matplotlib.pyplot as plt#岭回归
def ridgeRegres(xMat,yMat,lam=0.2):xTx = xMat.T*xMatdenom = xTx + eye(shape(xMat)[1])*lamif linalg.det(denom) == 0.0:print("This matrix is singular, cannot do inverse")returnws = denom.I * (xMat.T*yMat)return wsdef ridgeTest(xArr,yArr):xMat = mat(xArr); yMat=mat(yArr).TyMean = mean(yMat,0)yMat = yMat - yMean     #to eliminate X0 take mean off of Y#regularize X'sxMeans = mean(xMat,0)   #calc mean then subtract it offxVar = var(xMat,0)      #calc variance of Xi then divide by itxMat = (xMat - xMeans)/xVar #特征数据标准化处理 numTestPts = 30wMat = zeros((numTestPts,shape(xMat)[1]))for i in range(numTestPts):ws = ridgeRegres(xMat,yMat,exp(i-10))wMat[i,:]=ws.Treturn wMatabx,aby=loadDataSet('abalone.txt')
ridgeweights=ridgeTest(abx,aby)
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(ridgeweights)
plt.show()

4.2 lasso回归

增加如下约束时，普通的最小二乘法回归会得到与岭回归的一样的公式：

$\sum_{k=1}^{n}w_{k}^2\leqslant\lambda$

使用普通的最小二乘法回归在当两个或更多的特征相关时，可能会得出一个很大的正系数和一个很大的负系数。正是因为此约束条件存在，使用岭回归可以避免此问题。

与岭回归类似，另一个缩减方法lasso也对回归系数做了限定，对应的约束条件

$\sum_{k=1}^{n}\left | w_k \right |\leqslant\lambda$

此处约束条件使用绝对值取代了平方和。当 $\lambda$ 足够小的时候，一些系数会因此被迫缩减到0，这个特性可以帮我们更好的理解数据。

4.3前向逐步回归

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

'''
数据标准化，使其分布满足0均值和单位方差
在每轮跌倒过程中：设置当前最小误差lowestError为正无穷对每个特征：增大或缩小：改变一个系数得到一个新的w计算新w下的误差如果误差Error小于当前最小误差lowestError：设置wbest等于当前的w将w设置为新的wbest
'''

def regularize(xMat):#regularize by columnsinMat = xMat.copy()inMeans = mean(inMat,0)   #calc mean then subtract it offinVar = var(inMat,0)      #calc variance of Xi then divide by itinMat = (inMat - inMeans)/inVarreturn inMat#逐步线性回归算法的实现，与lasso做法相似单计算简单
def stageWise(xArr,yArr,eps=0.01,numIt=100):xMat = mat(xArr); yMat=mat(yArr).TyMean = mean(yMat,0)yMat = yMat - yMean     #can also regularize ys but will get smaller coefxMat = regularize(xMat)m,n=shape(xMat)#returnMat = zeros((numIt,n)) #testing code removews = zeros((n,1)); wsTest = ws.copy(); wsMax = ws.copy()for i in range(numIt):print(ws.T)lowestError = inf; #正无穷for j in range(n):for sign in [-1,1]:wsTest = ws.copy()wsTest[j] += eps*signyTest = xMat*wsTestrssE = rssError(yMat.A,yTest.A)if rssE < lowestError:lowestError = rssEwsMax = wsTestws = wsMax.copy()stagedata2=stageWise(abx,aby,0.001,5000)xmat=mat(abx)
ymat=mat(aby).T
xmat=regularize(xmat)
ym=mean(ymat,0)
ymat=ymat-ym
weights=standRegres(xmat,ymat.T)
print(weights.T)...
[[ 0.044 -0.011  0.12   0.022  2.023 -0.963 -0.105  0.187]]
[[ 0.0430442  -0.02274163  0.13214087  0.02075182  2.22403814 -0.99895312-0.11725427  0.16622915]]

可以看到进过5000次迭代后，逐步线性回归算法与常规的最小二乘法效果类似。该算法主要优点可以找出重要的特征，这样就可能及时停止对那些不重要特征的收集。

5，权衡偏差和方差

模型和测量之间存在差异，就出现了误差。对复杂的过程进行简化，这将导致差异的发生。

一般采用降低核大小和采用缩减法。缩减法通过吧一些特征的回归系数缩减到0，同时也就减少了模型的复杂度。

6，实例：预测乐高玩具套装的价格

用回归法预测乐高套装的价格：

1）收集数据：任何方式

2）准备数据：从返回的JSON中抽取价格

3）分析数据：可视化并观察数据

4）训练算法：构建不同的模型，采用逐步线性回归和直接的线性回归模型

5）测试算法：使用交叉验证来测试不同的模型，分析哪个效果最好

6）使用算法：这次练习的目标就是生产数据模型

6.1收集数据，

from time import sleep
import json
import urllib
def searchForSet(retX, retY, setNum, yr, numPce, origPrc):sleep(10) #防止短时间内有过多的api调用myAPIstr = 'AIzaSyD2cR2KFyx12hXu6PFU-wrWot3NXvko8vY'searchURL = 'https://www.googleapis.com/shopping/search/v1/public/products?key=%s&country=US&q=lego+%d&alt=json' % (myAPIstr, setNum)pg = urllib.request.urlopen(searchURL)retDict = json.loads(pg.read()) #打开和解析操作for i in range(len(retDict['items'])):try:currItem = retDict['items'][i]if currItem['product']['condition'] == 'new':newFlag = 1else: newFlag = 0listOfInv = currItem['product']['inventories']for item in listOfInv:sellingPrice = item['price']#信息过滤可以采用贝叶斯来判断，此处使用剪短的启发式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, origPrc])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, 269.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)

2，训练算法：建立模型

使用 standRegres ，在新数据集上进行回归处理拿到ws,得到公式

#交叉验证测试回归
def crossValidation(xArr,yArr,numVal=10):m = len(yArr)                           indexList = range(m)errorMat = zeros((numVal,30))#create error mat 30columns numVal rowsfor i in range(numVal):trainX=[]; trainY=[]testX = []; testY = []random.shuffle(indexList) #对数据进行混洗for j in range(m):#create training set based on first 90% of values in indexListif j < m*0.9: trainX.append(xArr[indexList[j]])trainY.append(yArr[indexList[j]])else:testX.append(xArr[indexList[j]])testY.append(yArr[indexList[j]])wMat = ridgeTest(trainX,trainY)    #get 30 weight vectors from ridgefor k in range(30):#loop over all of the ridge estimatesmatTestX = mat(testX); matTrainX=mat(trainX)meanTrain = mean(matTrainX,0)varTrain = var(matTrainX,0)matTestX = (matTestX-meanTrain)/varTrain #regularize test with training paramsyEst = matTestX * mat(wMat[k,:]).T + mean(trainY)#test ridge results and storeerrorMat[i,k]=rssError(yEst.T.A,array(testY))#保存每个λ对应的多个误差值#print errorMat[i,k]meanErrors = mean(errorMat,0)#calc avg performance of the different ridge weight vectorsminMean = float(min(meanErrors))bestWeights = wMat[nonzero(meanErrors==minMean)]#can unregularize to get model#when we regularized we wrote Xreg = (x-meanX)/var(x)#we can now write in terms of x not Xreg:  x*w/var(x) - meanX/var(x) +meanYxMat = mat(xArr); yMat=mat(yArr).TmeanX = mean(xMat,0); varX = var(xMat,0)unReg = bestWeights/varXprint("the best model from Ridge Regression is:\n",unReg)print ("with constant term: ",-1*sum(multiply(meanX,unReg)) + mean(yMat))

通过交叉验证落后得到的的回归系数，存在大小不一，如果在特征很多的情况下，可以方便的选择出那些特征是关键的，那些特征是不重要的。

7，总结

回归与分类不同，分类预测离散型变量，回归预测连续性变量。求特征对应的最佳回归系数的方法是最小化误差的平方和。前提是 $X^TX$ 的逆存在。数据集上计算的回归方程不一定是最佳的，可以使用预测yhat和原始值y，的相关性来度量回归方程的好坏。

当样本数比特征数还少的时候， $X^TX$ 的逆不能直接计算，可以考虑是用岭回归，

岭回归是缩减法的一种，相当于对回归系数的大小做了限制。另一种锁甲你发是lasso,lasso难以求解，可以考虑使用简单的逐步线性回归来求得近似结果

缩减法可以看成是对一个模型增加编差的同时减少方差。偏差方差折中是一个重要的概念