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

 

课程的所有数据和代码在我的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'))-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的预测值
    #返回归系数

 

 

局部加权线性回归

#局部加权线性回归
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，方差为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

 

岭回归

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-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

 

 

训练模型、交叉验证

#训练算法：建立模型
#交叉验证测试岭回归
#@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))

 

预测玩具价格

#收集数据
#添加时间函数库
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=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)