预测数值型数据：回归 源码分析（2）

4. 缩减系数来“理解”数据

4.1 岭回归

如果数据的特征比样本点还多，那么就不能使用线性回归，因为在计算$(X^TX)^{-1}$的时候会出错。也就是输入数据的矩阵$X$不是满秩矩阵，非满秩矩阵在求逆时会出现问题，为此有了岭回归。

简单说来，岭回归就是在矩阵$X^TX$上加一个$\lambda I$从而使得矩阵非奇异，进而能对$X^TX+\lambda I$求逆。$\lambda$是一个用户定义的数值。在这种情况下，回归系数的计算公式将变成：

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

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

这段话就是岭回归最大的特点和用处。。。。

岭回归代码的原理：
这里也是通过预测误差最小化得到$\lambda$数据获取之后,首先抽一部分数据用于测试，剩余的作为训练集用于训练参数$w$训练完毕后在测试集上测试预测性能。通过选取不同的$\lambda$来重复上述测试过程，最终得到一个使预测误差最小的$\lambda$。

# -*- coding: utf-8 -*-
"""
Created on Wed Oct 25 16:49:50 2017
"""
from numpy import *
import matplotlib.pyplot as plt

# 数据导入函数
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 ridgeRegres(xMat,yMat,lam=0.2): # lambda是关键字，此处不能用
    xTx = xMat.T*xMat  # 构建矩阵x`x
    denom = xTx + eye(shape(xMat)[1])*lam  # shape(xMat)[1]得到的是特征数
    if linalg.det(denom) == 0.0: # linalg.det()用来计算行列式，检查其是否为0
        print "This matrix is singular, cannot do inverse"
        return
    ws = denom.I * (xMat.T*yMat)  # 如果矩阵非奇异，计算回归系数并返回
    return ws  # 返回的是回归系数

# 用于在一组lambda上测试结果
def ridgeTest(xArr,yArr):
    xMat = mat(xArr); yMat=mat(yArr).T
    yMean = mean(yMat,0)  # axis = 0 压缩行，对各列求均值
    yMat = yMat - yMean     # 减去均值
    xMeans = mean(xMat,0)   # 得到每列特征的均值
    xVar = var(xMat,0)      # 得到每列的方差，此处var()是求得方差
    xMat = (xMat - xMeans)/xVar # 特征减去各自的均值并除以方差，就是标准化的过程
    numTestPts = 30
    wMat = zeros((numTestPts,shape(xMat)[1]))
    for i in range(numTestPts): # 30个不同的lambda参数下调用ridgeregre()
        ws = ridgeRegres(xMat,yMat,exp(i-10)) # 这里lambda应以指数级变化，以观察到有效的结果对比
        #print 'shape ws:',shape(ws)
        wMat[i,:]=ws.T # 把系数转置后成为行向量后赋给wmat[]
    return wMat # 返回最终得到的30次不同的lambda得到的系数


# 岭回归
abX,abY=loadDataSet('abalone.txt')
ridgeWeights=ridgeTest(abX,abY)
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(ridgeWeights) # 对于多维的画图，直接plt()即可
plt.show()

这里要注意的地方就是对特征进行标准化处理：使每维特征具有相同的重要性
具体的做法是所有特征都减去各自的均值并除以方差。

运行结果：

这里想说下结果图的横纵坐标：

机器学习实战上面的横坐标为$log(\lambda)$，而此处由结果图可以看出横坐标是$10+log(\lambda)$，因为横坐标是$0-30$，也就是以上图的纵坐标为刻度的，例如对应0点的每个值就是$\lambda$最小时的8个特征对应的系数值，而此时的$\lambda=e^{(0-10)}=e^{-10}$，也就是$10+log(\lambda)=10+log(e^{-10})=0$既是此时的横坐标，结果图的其他刻度值同样也是这样的一个对应关系。

4.2 前向逐步回归（贪心算法）

当最小二乘法回归加上约束条件时：$\sum_{k=1}^{n}W_k^2\leqslant \lambda$就可以得到和岭回归一样的公式。
为什么要这样做？
因为在使用普通最小二乘法时，在当两个或者更多的特征相关时，可能会得到一个很大的正系数和一个很大的负系数，但限制了所有回归系数的平方和不大于$\lambda$时，就可以解决这个问题。

再说下lasso：
与岭回归类似，lasso也是对回归做了限定，对应的约束：$\sum_{k=1}^{n}|W_k|\leqslant \lambda$，这里是用的绝对值，但效果却和岭回归有很大的差别：

• 绝对值：在$\lambda$很小时，一些系数会被迫缩减到0，这样就可以洗漱数据，得到更有用的特征，也就是能更好的理解数据，类似于正则化中的范数。

前向逐步回归：可以得到和lasso差不多的效果，但更加简单，它是一种贪心算法，即每一步都尽可能的减少误差。

# -*- coding: utf-8 -*-
"""
Created on Wed Oct 25 16:49:50 2017
"""
from numpy import *
import matplotlib.pyplot as plt

# 数据导入函数
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 ridgeRegres(xMat,yMat,lam=0.2): # lambda是关键字，此处不能用
    xTx = xMat.T*xMat  # 构建矩阵x`x
    denom = xTx + eye(shape(xMat)[1])*lam  # shape(xMat)[1]得到的是特征数
    if linalg.det(denom) == 0.0: # linalg.det()用来计算行列式，检查其是否为0
        print "This matrix is singular, cannot do inverse"
        return
    ws = denom.I * (xMat.T*yMat)  # 如果矩阵非奇异，计算回归系数并返回
    return ws  # 返回的是回归系数

# 用于在一组lambda上测试结果
def ridgeTest(xArr,yArr):
    xMat = mat(xArr); yMat=mat(yArr).T
    yMean = mean(yMat,0)  # axis = 0 压缩行，对各列求均值
    yMat = yMat - yMean     # 减去均值
    xMeans = mean(xMat,0)   # 得到每列特征的均值
    xVar = var(xMat,0)      # 得到每列的方差，此处var()是求得方差
    xMat = (xMat - xMeans)/xVar # 特征减去各自的均值并除以方差，就是标准化的过程
    numTestPts = 30
    wMat = zeros((numTestPts,shape(xMat)[1]))
    for i in range(numTestPts): # 30个不同的lambda参数下调用ridgeregre()
        ws = ridgeRegres(xMat,yMat,exp(i-10)) # 这里lambda应以指数级变化，以观察到有效的结果对比
        #print 'shape ws:',shape(ws)
        wMat[i,:]=ws.T # 把系数转置后成为行向量后赋给wmat[]
    return wMat # 返回最终得到的30次不同的lambda得到的系数

# 标准化处理：方差为0，方差为1
def regularize(xMat):   
    inMat = xMat.copy()
    inMeans = mean(inMat,0)   
    inVar = var(inMat,0)      
    inMat = (inMat - inMeans)/inVar
    return inMat

# 计算预测误差的大小
def rssError(yArr,yHatArr): # yArr和yHatArr都需要是数组
    return ((yArr-yHatArr)**2).sum()


# 前向逐步线性回归(贪心算法)
def stageWise(xArr,yArr,eps=0.01,numIt=100): # eps 每次迭代的步长，numIt 迭代的次数
    xMat = mat(xArr); yMat=mat(yArr).T
    yMean = mean(yMat,0) # axis = 0 压缩行，对各列求均值
    yMat = yMat - yMean     
    xMat = regularize(xMat) # 特征标准化为均值为0，方差为1
    m,n=shape(xMat) # 特征矩阵的行和列
    returnMat = zeros((numIt,n)) # 返回的矩阵100*n
    ws = zeros((n,1)); wsTest = ws.copy(); wsMax = ws.copy() 
    for i in range(numIt): # 迭代numIt次
        #print 'ws.T:',ws.T # 打印出来分析效果
        lowestError = inf; # 每一次迭代误差初始值设为正无穷
        for j in range(n): # 在所有特征上循环
            for sign in [-1,1]: # 每个特征运行两次
                wsTest = ws.copy() # 每次的ws都是上一次最优后得到的
                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 # 最后返回的是贪心算法得到的系数矩阵


# 前向逐步回归
xArr,yArr=loadDataSet('abalone.txt')
wsMat_1=stageWise(xArr,yArr,0.01,200)
#print 'wsMat:',wsMat
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(wsMat_1) # 对于多维的画图，直接plt()即可
ax.set_title('Step length=0.01, Number of iterations:200')
plt.show()

wsMat_2=stageWise(xArr,yArr,0.001,5000)
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(wsMat_2) # 对于多维的画图，直接plt()即可
ax.set_title('Step length=0.001, Number of iterations:5000')
plt.show()

运行结果：

其中返回的系数矩阵为：

• 其中的第一列和第六列都为0，即说明这两个特征的系数为0，也就是该特征不是主要特征，从而起到降维的作用。
• 还有就是可以看出在参数步长为0.01时，一段是时间后饱和后，一些系数就在特征值之间来回震荡，这是因为步长太大的缘故。由对比不同的步长和迭代次数可以看出系数逐步稳定的过程，不再震荡。
• 逐步线性回归可以帮助理解现有的模型，当构建好一个模型后，可以运算该算法找出重要特征，这样就可以停止不重要的特征的收集。如果用于测试的话，该算法可以通过迭代后构建很多的同类模型，可以使用类似于10折交叉验证法比较这些模型，最终选择使误差最小的模型。

5. 方差和偏差的简单理解

方差指的是模型之间的差异，而偏差指的是模型预测值和数据之间的差异！！！

当应用缩减方法时，模型增加了偏差，同时却减少了模型的方差。

可以看出将一些系数缩减到很小的值或直接缩减为0，这是一个减少模型复杂度的过程，但同时也是一个增大模型偏差的过程。

方差是可以度量的。如果从鲍鱼数据中取一个随机样本集（例如取其中100个数据）并用线性模型拟合，将会得到一组回归系数。同理，再取出另一组随机样本集并拟合，将会得到另一组回归系数。这些系数间的差异大小也就是模型方差大小的反映

6. 预测乐高玩具套装的价格

理论应用：首先从拍卖站点抽取一些数据，再使用一些回归法进行实验来为数据找到最佳的岭回归模型。这样就可以通过实际效果来看看偏差和方差间的折中效果。

算法流程：

1 收集数据：用google shopping的api收集数据
2 准备数据：从返回的json数据中抽取价格
3 分析数据：可视化并观察数据
4 训练算法：构建不同的模型，采用岭回归和普通线性回归训练模型
5 测试算法：使用交叉验证来测试不同的模型，选择效果最好的模型

• 收集数据：使用Google 购物的API来获取玩具套装的相关信息和价格，可以通过urllib2发送http请求，API将以JSON格式返回需要的产品信息，python的JSON解析模块可以帮助我们从JSON格式中解析出所需要的数据。收集数据的代码如下：

由于对爬数据不太会和其他原因，这里没有运行出来，但我觉得还是有必要对源码好好分析下，这里仅给出代码解析：

# -*- coding: utf-8 -*-
from time import sleep
import json
import urllib2

# 购物信息的获取函数
def searchForSet(retX, retY, setNum, yr, numPce, origPrc): # 调用API并数据抽取
    sleep(10) # 休息10秒钟，防止短时间内有过多的API调用
    # 拼接查询的url字符串，添加API的Key和待查询的套装信息
    searchURL = 'https://www.googleapis.com/shopping/search/v1/public/products?key=%s&country=US&q=lego+%d&alt=json' % (myAPIstr, setNum)
    pg = urllib2.urlopen(searchURL)   # 打开 URL 等待返回数据
    retDict = json.loads(pg.read())   # 利用json打开和解析url获得的数据，数据信息存入字典中
    print ('i am here!')

    # 遍历数据的每一个条目
    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'] # 得到该条目玩具商品的价格
                if  sellingPrice > origPrc * 0.5: # 价格低于原价的50%视为不完整套装
                    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)

lgx=[];lgy=[]
print ('setDataCollect:',setDataCollect(lgx,lgy))

其中不完整的套装的检索是用的启发式！！

• 训练算法，建立模型

要达到的目的是构建的模型可以对售价做出预测，并帮助理解现有数据。

现用岭回归进行模型的建立，上一篇博客讲过如何对系数进行缩减，下面的代码将是如何用缩减法确定最佳回归系数。

# 交叉验证测试岭回归
# xArr:从网站中获得的玩具套装样本数据，yArr：样本对应的出售价格，numVal：交叉验证的次数
def crossValidation(xArr,yArr,numVal=10): 
    m = len(yArr)     #  获取样本数                       
    indexList = range(m)
    errorMat = zeros((numVal,30)) # 10次测试  每次有30组回归系数 可以得到误差  
    for i in range(numVal): # 测试次数，默认为10次
        trainX=[]; trainY=[] # 训练数据集和标签 
        testX = []; testY = [] # 测试数据集和标签
        random.shuffle(indexList) # 混洗索引列表，以实现训练集或测试集数据点的随机选取
        for j in range(m): # 遍历每个样本
            if j < m*0.9:  # 数据集90%作为训练集
                trainX.append(xArr[indexList[j]])
                trainY.append(yArr[indexList[j]])
            else: # 剩余10%作为测试集
                testX.append(xArr[indexList[j]])
                testY.append(yArr[indexList[j]])
        wMat = ridgeTest(trainX,trainY)   # 得到岭回归的所有回归系数，得到了30组不同的回归系数
        for k in range(30):  # 对于30组不同的岭回归得到的回归系数进行测试，计算误差，选取最好的
            matTestX = mat(testX); matTrainX=mat(trainX) 
            meanTrain = mean(matTrainX,0) # 训练数据 均值
            varTrain = var(matTrainX,0)   # 训练数据 方差
            # 特征减去各自的均值并除以方差，就是标准化的过程，岭回归需要使用标准化的数据，
            # 因此数据也需要使用与训练集相同的参数来标准化
            matTestX = (matTestX-meanTrain)/varTrain  # 用训练集的参数将测试数据标准化 
            # 之所以加上训练标签的均值也是为了标准化一致，从而得到的测试集预测值是其真实预测值
            yEst = matTestX * mat(wMat[k,:]).T + mean(trainY) # yEst对应的是所有测试样本在每个lambda下的预测值
            # 计算误差,errorMat默认是10*30，保存的是10行交叉验证，每行的30个元素分别是对应的lambda所有测试样本预测值的误差和
            errorMat[i,k]=rssError(yEst.T.A,array(testY)) 
    meanErrors = mean(errorMat,0)  # 计算误差估计值的均值，此处为10折，每个lambda对应10个误差
    minMean = 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) +meanY
    xMat = mat(xArr); yMat=mat(yArr).T
    meanX = mean(xMat,0); varX = var(xMat,0)
    # 数据标准化还原操作
    unReg = bestWeights/varX
    print "the best model from Ridge Regression is:\n",unReg
    print "with constant term: ",-1*sum(multiply(meanX,unReg)) + mean(yMat) # 还原计算预测结果

这里采用岭回归来训练模型，并且采用交叉验证的方法来求出每个λ对应的测试误差的均值，最后分析选出预测误差最小的回归模型。

注意：

• 这里对于数据集采用随机的方式（random.shffle()）选取训练集和测试集，训练集占数据总数的90%，测试集剩余的10%。采取这种方式的原因是，便于我们进行多次交叉验证，得到不同的训练集和测试集.
• 我们知道岭回归中会选取多个不同的λ值，来找到预测误差最小的模型；此外，算法中采用交叉验证的方法，所以对于每一个λ对应着多个测试误差值，所以在分析预测效果最好的λ之前，需要先对每个λ对应的多个误差求取均值。
• 岭回归算法需要对训练集数据的每一维特征进行标准化处理，那么为保证结果的准确性，也需要对测试集进行和训练集相同的标准化操作，即测试集数据特征减去训练集该维度特征均值，再除以训练集该维度特征方差
• 因为采用岭回归算法时，对数据进行了标准化处理，而标准的回归算法则没有，所以在代码最后我们还是需要将数据进行还原，这样便于分析比较二者的真实数据的预测误差。

以上是很重要的几点要注意的地方！！！

从机器学习实战的运行结果可以看出具体的缩减过程中系数的变化情况：最后得到的回归系数是经过不同程度的衰减得到的，大的特征系数可以看做是最重要特征，在预测时起最主要作用。特征对应的系数值越大，那么其对预测的决定作用也就越大。如果某一维度系数值为0，则表明该特征在预测结果中不起作用，可以被视为不重要特征。

　　所以，这种缩减的分析方法还是比较有用的，因为运算这些算法可以帮助我们充分理解和挖掘大量数据中的内在规律。当特征数较少时可能效果不够明显，而当特征数相当大时，我们就可以据此了解特征中哪些特征是关键的，哪些是不重要的，这就为我们节省不少成本和损耗。

总结：

(1) 回归与分类的区别，前者预测连续型变量，后者预测离散型变量；回归中求最佳系数的方法常用的是最小化误差的平方和；如果$x^Tx$可逆，那么回归算法可以使用；可以通过预测值和原始值的相关系数来度量回归方程的好坏

(2) 当特征数大于样本总数时，$x^Tx$不可逆，即便当样本总数大于特征数，$x^Tx$的逆仍有可能无法计算，因为特征可能高度相关，我们可以通过引入岭回归来保证能够求得回归系数。

(3) 另外一种缩减算法是，前向逐步回归算法，它是一种贪心算法，每一步通过修改某一维度特征方法来减小预测误差，最后通过多次迭代的方法找到最小误差对应的模型

(4) 缩减法可以看做是对一个模型增加偏差的同时减少方差，通过偏差方差折中的方法，可以帮助我们理解模型并进行改进，从而得到更好的预测结果

（5）当预测值和特征之间是非线性的关系时，这时线性的模型就难以拟合，但可以使用树结构来预测。

shuffle()的示例：

In [2]: a=range(10)

In [3]: random.shuffle(a)

In [4]: a
Out[4]: [0, 9, 1, 7, 3, 8, 2, 5, 4, 6]