本章首先介绍线性回归，包括其名称的由来和实现。接下来本章将讨论回归在“欠拟合”的情况下的缩减技术。最后将融合所有技术预测鲍鱼年龄和玩具售价。
线性回归
优点：结果易于理解，计算上不复杂
缺点：对非线性的数据拟合不好
适用数据类型：数值型和标称型数据
下面开始计算标准回归函数和数据导入函数
from 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.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).T
xTx = xMat.T*xMat
if linalg.det(xTx) == 0.0:#计算行列式
print "This matrix is singular, connot do inverse"
return
ws = xTx.I * (xMat.T*yMat)#.I是逆
ws = linalg.solve(xTx,xMat.T*yMat)#一上一句同义
return ws
下面开始测试数据
In [30]: import regression
...: xArr,yArr = regression.loadDataSet('ex0.txt')
...: ws = regression.standRegres(xArr,yArr)
In [34]: xArr[0:2]
Out[34]: [[1.0, 0.067732], [1.0, 0.42781]]
第一个值为1，即x0。第二个值x1，就是横坐标
In [35]: ws
Out[35]:
matrix([[ 3.00774324],
[ 1.69532264]])
ws存放的是回归系数。下面开始使用新的ws计算yHat：
xMat = mat(xArr)
yMat = mat(yArr)
yHat = xMat*ws
然后绘制散点图
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(xMat[:,1].flatten().A[0],yMat.T[:,0])
xCopy = xMat.copy()
xCopy.sort(0)#先排序
yHat = xCopy*ws
ax.plot(xCopy[:,1],yHat)
plt.show()
最佳拟合直线
在Python中，NumPy库提供了相关系数的计算方法：可以通过命令corrcoef(yEstimate,yActual)来计算相关性。
In [41]: yHat = xMat*ws
In [42]: corrcoef(yHat.T,yMat)
Out[42]:
array([[ 1. , 0.98647356],
[ 0.98647356, 1. ]])
线性回归有时会出现欠拟合的现象，因为她求的是具有最小均方误差的无偏估计。所以有些时候允许在估计中引入一些偏差，从而降低预测的均分误差。
#局部加权线性回归函数
def lwlr(testPoint,xArr,yArr,k=1.0):
xMat = mat(xArr);yMat = mat(yArr).T
m = shape(xMat)[0]
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)
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 = zeros(m)
for i in range(m):
yHat[i] = lwlr(testArr[i],xArr,yArr,k)
return yHat
并输入如下命令
xArr,yArr = regression.loadDataSet('ex0.txt')
#对单点估计
In [53]: yArr[0]
Out[53]: 3.176513
In [58]: regression.lwlr(xArr[0],xArr,yArr,1.0)
Out[58]: matrix([[ 3.12204471]])
In [59]: regression.lwlr(xArr[0],xArr,yArr,0.001)
Out[59]: matrix([[ 3.20175729]])
yHat = regression.lwlrTest(xArr,xArr,yArr,0.003)#为了得到数据集所有的点
然后开始绘图，首先对xArr排序：
xMat=mat(xArr)
srtInd=xMat[:,1].argsort(0)
xSort=xMat[srtInd][:,0,:]
然后开始绘制
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(xSort[:,1],yHat[srtInd])
ax.scatter(xMat[:,1].flatten().A[0],mat(yArr).T.flatten().A[0], s=2,c='red')#c:散点的颜色，s：散点的大小 ，alpha:是透明程度，make:散点样式
plt.show()
下面是三种k的取值下的结果图：
k=0.003
k=0.01
k=1.0
k=0.003的时候显然是过拟合，k=1.0的时候和最小二乘法差不多。但局部加权线性回归也存在计算量的问题，因为他对每个点做预测时都必须使用整个数据集。
接下来，我们将回归用于真实数据。
#加入regression,py
def rssError(yArr,yHatArr):
return ((yArr-yHatArr)**2).sum()
先计算数据
abX,abY = regression.loadDataSet('abalone.txt')
yHat01 = regression.lwlrTest(abX[0:99],abX[0:99],abY[0:99],0.1)
yHat1 = regression.lwlrTest(abX[0:99],abX[0:99],abY[0:99],1)
yHat10 = regression.lwlrTest(abX[0:99],abX[0:99],abY[0:99],10)
然后计算错误指标
In [92]: regression.rssError(abY[0:99],yHat01.T)
Out[92]: 56.820227823583345
In [93]: regression.rssError(abY[0:99],yHat1.T)
Out[93]: 429.89056187016047
In [94]: regression.rssError(abY[0:99],yHat10.T)
Out[94]: 549.11817088250825
但是k=0.1显然是过拟合现象
我们使用新的数据来看预测效果
In [102]: yHat01 = regression.lwlrTest(abX[100:199],abX[0:99],abY[0:99],0.1)
...: regression.rssError(abY[100:199],yHat01.T)
...:
Out[102]: 23989.306318956587
In [103]: yHat1 = regression.lwlrTest(abX[100:199],abX[0:99],abY[0:99],1)
...: regression.rssError(abY[100:199],yHat1.T)
Out[103]: 573.52614418975463
In [104]: yHat10 = regression.lwlrTest(abX[100:199],abX[0:99],abY[0:99],10)
...: regression.rssError(abY[100:199],yHat10.T)
上述结果可以看出，核大小等于10的时候误差反而最小。
接下来和简单线性回归作比较：
In [105]: ws = regression.standRegres(abX[0:99],abY[0:99])
In [106]: yHat = mat(abX[100:199])*ws
In [107]: regression.rssError(abY[100:199],yHat.T.A)
Out[107]: 518.636315324674
简单线性回归达到了与局部加权线性回归类似的效果。这说明：必须在未知数据上比较效果才能选取最佳模型。
局部加权线性回归的问题在于，每次必须在整个数据集上运行。
下面使用岭回归处理数据
def ridgeRegres(xMat,yMat,lam=0.2):#计算回归函数，lam即Python保留关键字lambda
xTx = xMat.T*xMat
denom = xTx + eye(shape(xMat)[1])*lam#岭
if linalg.det(denom) == 0.0:
print "This martrix 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
numTestPts = 30
wMat = zeros((numTestPts,shape(xMat)[1]))
for i in range(numTestPts):
ws = ridgeRegres(xMat,yMat,exp(i-10))
wMat[i,:]=ws.T
return wMat
然后开始绘图
import regression
abX,abY = regression.loadDataSet('abalone.txt')
ridgeWeights = regression.ridgeTest(abX,abY)
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(ridgeWeights)
plt.show()
log(lambda)
在最左边，lambda最小的时候，可以得到所有系数的原始值（与线性回归一致）；在最右边，系数全部缩减为0；在中间部分的某值将可以取得最好的预测效果。为了定量地找到最佳参数值，还需要进行交叉验证。
前后逐步回归可以得到和lasso差不多的效果，但更加简单。它属于一种贪心算法，即每一步都尽可能减少误差。一开始，所有的权重都设为1，然后每一步所做的决策是对某个权重增加或者减少一个很小的值。
伪代码如下：
数据标准化，使其分布满足0均值和单位方差：
在每轮迭代过程中：
设置当前最小误差Lowe'sError为正无穷
对每个特征：
增大或缩小：
改变一个系数得到一个新的w
计算新W下的误差
如果误差Error小于当前最小误差Lowe'sError：设置Wbest等于当前的W
将W设置为新的Wbest
#前向逐步线性回归
def stageWise(xArr,yArr,eps=0.01,numIt=100):#逐步线性回归实现。输入数据，预测变量，步长，迭代次数
xMat = mat(xArr); yMat = mat(yArr).T
yMean = mean(yMat,0)#均值
yMat = yMat - yMean
xMat = regularize(xMat)#标准化
m,n = shape(xMat)
returnMat = zeros((numIt,n))
ws = 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*sign
yTest = xMat*wsTest
rssE = rssError(yMat.A,yTest.A)
if rssE < lowestError:
lowestError = rssE
wsMax = wsTest
ws = wsMax.copy()
returnMat[i,:] = ws.T
return returnMat
def regularize(xMat):
inMat = xMat.copy()
inMeans = mean(inMat,0)
inVar = var(inMat,0)
inMat = (inMat - inMeans)/inVar
return inMat
并且测试这段代码
In [22]: import regression
...: xArr,yArr = regression.loadDataSet('abalone.txt')
...: regression.stageWise(xArr,yArr,0.01,200)
...:
[[ 0. 0. 0. 0. 0. 0. 0. 0.]]
[[ 0. 0. 0. 0.01 0. 0. 0. 0. ]]
[[ 0. 0. 0. 0.02 0. 0. 0. 0. ]]
[[ 0. 0. 0. 0.03 0. 0. 0. 0. ]]
[[ 0. 0. 0. 0.04 0. 0. 0. 0. ]]
[[ 0. 0. 0. 0.05 0. 0. 0. 0. ]]
[[ 0. 0. 0. 0.06 0. 0. 0. 0. ]]
[[ 0. 0. 0.01 0.06 0. 0. 0. 0. ]]
[[ 0. 0. 0.01 0.06 0. 0. 0. 0.01]]
[[ 0. 0. 0.01 0.06 0. 0. 0. 0.02]]
#省略部分
[[ 0.05 0. 0.09 0.03 0.31 -0.64 0. 0.36]]
[[ 0.04 0. 0.09 0.03 0.31 -0.64 0. 0.36]]
[[ 0.05 0. 0.09 0.03 0.31 -0.64 0. 0.36]]
[[ 0.04 0. 0.09 0.03 0.31 -0.64 0. 0.36]]
[[ 0.05 0. 0.09 0.03 0.31 -0.64 0. 0.36]]
[[ 0.04 0. 0.09 0.03 0.31 -0.64 0. 0.36]]
[[ 0.05 0. 0.09 0.03 0.31 -0.64 0. 0.36]]
[[ 0.04 0. 0.09 0.03 0.31 -0.64 0. 0.36]]
可以发现w1和w6都是0，这表明他们不对目标值产生影响。在eps设为0.01的情况下，一段时间后系数就已经饱和并在特定值之间来回震荡，这是因为步长太大的缘故。
下面试试更小的步长和更多的步数：
regression.stageWise(xArr,yArr,0.001,5000)
array([[ 0. , 0. , 0. , ..., 0. , 0. , 0. ],
[ 0. , 0. , 0. , ..., 0. , 0. , 0. ],
[ 0. , 0. , 0. , ..., 0. , 0. , 0. ],
...,
[ 0.043, -0.011, 0.12 , ..., -0.963, -0.105, 0.187],
[ 0.044, -0.011, 0.12 , ..., -0.963, -0.105, 0.187],
[ 0.043, -0.011, 0.12 , ..., -0.963, -0.105, 0.187]])
接下来把这些结果和最小二乘法进行比较，后者的结果可以通过如下代码获得：
In [25]: xMat = mat(xArr)
...: yMat = mat(yArr).T
...: xMat = regression.regularize(xMat)
...: yM = mean(yMat,0)
...: yMat = yMat - yM
...: weights = regression.standRegres(xMat,yMat.T)
In [27]: weights.T
Out[27]:
matrix([[ 0.0430442 , -0.02274163, 0.13214087, 0.02075182, 2.22403814,
-0.99895312, -0.11725427, 0.16622915]])
逐步线性回国算法主要的优点在于他可以帮助人们理解现有的模型并做出改进。当应用缩减方法（如逐步回归或岭回归）时，模型也就增加了偏差，与此同时却减小了模型的方差。
下面的部分涉及谷歌api，但是原文太旧，找不到新版怎么操作，暂时搁浅...
