预测数值型数据：回归

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

2 局部加权线性回归

3 示例：预测鲍鱼的年龄

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

4.1 岭回归

4.2 lasso

4.3 前向逐步回归

5 权衡偏差与方差

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

6.1 收集数据：使用google购物的API

6.2 训练算法：建立模型

7 本章结束

本节涉及的相关代码和数据

预测数值型数据：回归

本章内容：

① 线性回归

② 局部加权线性回归

③ 岭回归和逐步线性回归

④ 预测鲍鱼年龄和玩具售价

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

线性回归的相关特点主要是

优点：结果易于理解，计算上不复杂

缺点：对非线性数据不友好

适用数据类型：数值型和标称型

回归的一般方法：

① 收集数据：采用任意方法收集数据

② 准备数据：回归需要数值型数据，标称型数据将会被转为二值型数据

③ 分析数据：绘出数据的可视化二维图将有助于对数据做出理解和分析，在采用缩减法求得新回归系数之后，可以将新拟合曲线绘在图上作为对比

④ 训练算法：找到回归系数

⑤ 测试算法：使用R2或者预测值和数据的拟合度，来分析模型的结果

⑥ 使用算法：使用回归，可以在给定输入的时候预测出一个数值，这是对分类方法的提升，因为这样可以预测连续性数据而不仅仅是离散的数据标签

采用最小平方误差

# 标准回归函数和数据导入函数
from numpy import *
# 加载数据
def loadDataSet(fileName):
    # print(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]))
            # print(lineArr)
        # print(lineArr)
        dataMat.append(lineArr)
        labelMat.append(float(curLine[-1]))
    # print(dataMat)
    return dataMat,labelMat

# 线性回归函数
def standRegres(xArr,yArr):
    # 数组转为矩阵
    xMat=mat(xArr)
    yMat=mat(yArr).T
    xTx=xMat.T*xMat
    # linalg.det()函数可以直接用来计算行列式
    if linalg.det(xTx)==0.0:
        print("this matrix is singular,cannot do inverse")
        return
    # matrix.I获得与给定矩阵相同大小的乘法逆
    ws=xTx.I*(xMat.T*yMat)
    # 返回权重
    return ws

调用上述函数：

xArr,yArr=loadDataSet('ex0.txt')
# xArr
# 回归系数
ws=standRegres(xArr,yArr)
ws

输出得到的结果为：

将得到的权重系数，画图表示为：

# 绘图
def plotData(xArr,yArr,ws):
    import matplotlib.pyplot as plt
    xMat = mat(xArr)
    yMat = mat(yArr)
    figure = plt.figure()
    ax = figure.add_subplot(111)
    # 取第二个特征绘图
    # flatten()函数转化成一维数组
    # matrix.A属性返回矩阵变成的数组，和getA()方法一样
    # 绘制点
    ax.scatter(xMat[:,1].flatten().A[0], yMat.T[:,0].flatten().A[0])
    # 返回给定数据的数组形式的拷贝
    xCopy = xMat.copy()
    # 升序排序
    xCopy.sort(0)
    print (ws.shape)
    yHat = xCopy * ws # yHat 表示拟合直线的纵坐标，用回归系数求出
    ax.plot(xCopy[:,1], yHat, c = 'green')
    plt.show()


plotData(xArr,yArr,ws)

得到的输出图像为：

可以看到该拟合，虽然在数据计算上结果较好，但可能单纯的线性拟合并不能够得到很好的结果，因为原始数据很明显有一个波浪形的分布，而不单单是线性分布。

计算该拟合结果的相关系数为：

# 计算相关系数
def calcCorrel(xArr,yArr,ws):
    xMat = mat(xArr)
    yMat = mat(yArr)
    yHat=xMat*ws
    # 转置yHat向量，以保证都是行向量
    # corrcoef函数计算皮尔逊相关系数
    correl=corrcoef(yHat.T,yMat)
    return correl

calcCorrel(xArr,yArr,ws)

得到的结果为：

2 局部加权线性回归

线性回归的一个问题就是有可能出现欠拟合的现象，因为他求的是具有最小均方误差的无偏估计。显而易见，如果模型欠拟合将不能取得最好的预测效果。所以有些方法允许在估计中引入一些偏差，从而降低预测的均方误差。

其中一个方法就是局部加权线性回归（LWLR）。我们给待预测点附近的每个点赋予一定的权重，然后在这个基础上基于最小均方差来进行普通的回归。与KNN一样，这种算法每次预测均需事先选取处对应的数据子集。LWLR使用‘核’来对附近的点赋予更高的权重。核的类型可以自由选择，常用的有高斯核。

def lwlr(testPoint,xArr,yArr,k=1.0):
    xMat=mat(xArr)
    yMat=mat(yArr).T
    m=shape(xMat)[0]
    # eye()函数返回一个二维数组，对角线为1，其余地方为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

调用上述函数：

lwlr(xArr[0],xArr,yArr)

得到的输出结果为：

画出该拟合的结果曲线为：

def plotDatalwlr(xArr,yArr,k=0.01):
    yHat=lwlrTest(xArr,xArr,yArr,k)
    xMat=mat(xArr)
    # 对样本x升序排序并返回索引
    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')
    plt.show()
xArr,yArr=loadDataSet('ex0.txt')
plotDatalwlr(xArr,yArr)

得到的拟合图像为：

局部加权线性回归也存在一个问题，即增加了计算量，因为他对每个点做预测时都必须使用整个数据集。如果避免这些计算将可以减少程序运行时间，从而减缓因计算量增加带来的问题。

3 示例：预测鲍鱼的年龄

# 计算误差
def rssError(yArr,yHatArr):
    # **表示幂运算，前面为底数，后面为指数
    return ((yArr-yHatArr)**2).sum()

abX,abY=loadDataSet('abalone.txt')

plotDatalwlr(abX[0:99],abY[0:99],0.1)
plotDatalwlr(abX[0:99],abY[0:99],1)
plotDatalwlr(abX[0:99],abY[0:99],10)

# 训练集
yHat01_=lwlrTest(abX[0:99],abX[0:99],abY[0:99],0.1)
yHat1_=lwlrTest(abX[0:99],abX[0:99],abY[0:99],1)
yHat10_=lwlrTest(abX[0:99],abX[0:99],abY[0:99],10)

err01_=rssError(abY[0:99],yHat01_.T)
err1_=rssError(abY[0:99],yHat1_.T)
err10_=rssError(abY[0:99],yHat10_.T)

print(err01_,err1_,err10_)

# 测试集
yHat01=lwlrTest(abX[100:199],abX[0:99],abY[0:99],0.1)
yHat1=lwlrTest(abX[100:199],abX[0:99],abY[0:99],1)
yHat10=lwlrTest(abX[100:199],abX[0:99],abY[0:99],10)

err01=rssError(abY[100:199],yHat01.T)
err1=rssError(abY[100:199],yHat1.T)
err10=rssError(abY[100:199],yHat10.T)

print(err01,err1,err10)

得到的输出结果为：

可以看到，使用较小的核可以获得较低的误差，但是使用最小的核可能会造成过拟合，对新数据不一定能够达到最好的效果

接下来再和简单的线性回归做个比较

ws=standRegres(abX[0:99],abY[0:99])
yHat=mat(abX[100:199])*ws
rssError(abY[100:199],yHat.T.A)

得到结果为：

简单线性回归达到了与局部加权线性回归类似的效果，这也表明：必须在未知数据上比较效果才能够选取到最佳模型

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

如果数据的特征比样本点还多，也就是说输入数据的矩阵X不是满秩矩阵。非满秩矩阵在求逆时会出现问题，因此，引入了“岭回归”的概念，这就是第一种缩减方法，接着是lasso方法，效果很好但复杂，最后介绍了第二种缩减方法，称为前向逐步回归，可以得到与lasso差不多的效果，也更容易实现。

4.1 岭回归

岭回归就是在矩阵XTX上加入一个λI从而使矩阵非奇异，进而就能对矩阵求逆，其中矩阵I时一个单位矩阵，λ由用户自定义。

缩减方法可以去掉不重要的参数，因此能够更好的理解数据，此外，与简单的线性回归相比，缩减法能够取得更好的效果。

# 用于计算回归系数
# 实现给定系数下的岭回归
def ridgeRegres(xMat,yMat,lam=0.2):
    xTx=xMat.T*xMat
    denom=xTx+eye(shape(xMat)[1])*lam
    if linalg.det(denom)==0.0:
        print("this matrix is singlar,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)
    # var()函数求样本方差的无偏估计值，如果参数是1，就是有偏估计
    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

abX,abY=loadDataSet('abalone.txt')
ridgeWeights=ridgeTest(abX,abY)

# 绘制出回归系数与log（λ）的关系
import matplotlib.pyplot as plt
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(ridgeWeights)
plt.show()

得到的结果为：

4.2 lasso

lasso对回归系数做出了限定，使用绝对值取代了平方和，在λ足够小的时候，一些系数会因此被迫缩减到0，这个特性将有助于理解数据

4.3 前向逐步回归

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

伪代码：

数据标准化，使其分布满足0均值和单位方差

在每轮迭代中：

设置当前最小误差lowestError为正无穷

对每个特征：

增大或减小：

改变一个系数得到一个新的w

计算新w下的误差

如果误差Error小于当前最小误差lowestError：

设置Wbest等于当前w

将w设置为新的Wbest

# 实现矩阵归一化
def regularize(xMat):
    inMat = xMat.copy()
    # 得到平均数，压缩行，对每一列求平均值
    inMeans = mean(inMat,0)  
    
    inVar = var(inMat,0)    
    inMat = (inMat - inMeans)/inVar
    return inMat

# 前向逐步线性回归
def stageWise(xArr,yArr,eps=0.01,numIt=100):
    xMat=mat(xArr)
    yMat=mat(yArr).T
    # mean()函数求平均值
    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

调用函数：

xArr,yArr=loadDataSet('abalone.txt')
stageWise(xArr,yArr,0.01,200)

得到的输出结果为：

与最小二乘法的结果进行对比得到

# 最小二乘法的结果
xMat=mat(xArr)
yMat=mat(yArr).T
xMat=regularize(xMat)
yM=mean(yMat,0)
yMat=yMat-yM
weights=standRegres(xMat,yMat.T)
weights.T

结果为：

5 权衡偏差与方差

任何时候，一旦发现模型和测量值直接存在差异，就说明出现了误差。当考虑模型中的“噪声”或者说误差时，必须考虑其来源。可能会因为对复杂的模型进行简化，这将导致在模型和测量值之间出现“噪声”或者“误差”，若无法理解数据的真实生成过程，也会导致差异的发生。

一般认为，误差一般由三个部分组成：偏差、测量误差和随机噪声

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

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

①收集数据：用google Shopping的APi收集数据

②准备数据：从返回的json数据中抽取价格

③分析数据：可视化并观察数据

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

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

⑥使用算法：这次联系的目的就是生成数据模型

6.1 收集数据：使用google购物的API

from bs4 import BeautifulSoup
def scrapePage(retX, retY, inFile, yr, numPce, origPrc):
     # 打开并读取HTML文件
    with open(inFile, encoding='utf-8') as f:
        html = f.read()
    # beautiful函数将复杂HTML文档转换成一个复杂的树形结构
    soup = BeautifulSoup(html)
    i = 1
    # 根据HTML页面结构进行解析
    currentRow = soup.find_all('table', r = "%d" % i)
    while(len(currentRow) != 0):
        currentRow = soup.find_all('table', r = "%d" % i)
        title = currentRow[0].find_all('a')[1].text
        lwrTitle = title.lower()
        # 查找是否有全新标签
        if (lwrTitle.find('new') > -1) or (lwrTitle.find('nisb') > -1):
            newFlag = 1.0
        else:
            newFlag = 0.0
        # 查找是否已经标志出售，我们只收集已出售的数据
        soldUnicde = currentRow[0].find_all('td')[3].find_all('span')
        if len(soldUnicde) == 0:
            print("商品 #%d 没有出售" % i)
        else:
            # 解析页面获取当前价格
            soldPrice = currentRow[0].find_all('td')[4]
            priceStr = soldPrice.text
            priceStr = priceStr.replace('$','')
            priceStr = priceStr.replace(',','')
            if len(soldPrice) > 1:
                priceStr = priceStr.replace('Free shipping', '')
            sellingPrice = float(priceStr)
            # 去掉不完整的套装价格
            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)
        i += 1
        currentRow = soup.find_all('table', r = "%d" % i)
 
def setDataCollect(retX, retY):
    scrapePage(retX, retY, './lego/lego8288.html', 2006, 800, 49.99)  # 2006年的乐高8288,部件数目800,原价49.99
    scrapePage(retX, retY, './lego/lego10030.html', 2002, 3096, 269.99)  # 2002年的乐高10030,部件数目3096,原价269.99
    scrapePage(retX, retY, './lego/lego10179.html', 2007, 5195, 499.99)  # 2007年的乐高10179,部件数目5195,原价499.99
    scrapePage(retX, retY, './lego/lego10181.html', 2007, 3428, 199.99)  # 2007年的乐高10181,部件数目3428,原价199.99
    scrapePage(retX, retY, './lego/lego10189.html', 2008, 5922, 299.99)  # 2008年的乐高10189,部件数目5922,原价299.99
    scrapePage(retX, retY, './lego/lego10196.html', 2009, 3263, 249.99)  # 2009年的乐高10196,部件数目3263,原价249.99

调用上面函数：

lgX=[]
lgY=[]
setDataCollect(lgX,lgY)

得到结果：

6.2 训练算法：建立模型

为上面收集到的数据建立模型。构建出来的模型可以对售价做出预测

# 创建一个全为1 的矩阵
lgX1=mat(ones((63,5)))
print(lgX1[0])
# 将原数据矩阵lgX复制到新数据矩阵lgX1的第一列到第五列
lgX1[:,1:5]=mat(lgX)
# 确认数据复制的正确性
print(lgX[0])
print(lgX1[0])
ws=standRegres(lgX1,lgY)
ws

得到的输出结果为：

交叉测试验证岭回归

# 交叉测试验证岭回归
                            #  交叉验证的次数
def crossValidation(xArr,yArr,numVal=10):
    m=len(yArr)
    indexList=list(range(m))
    errorMat=zeros((numVal,30))
    for i in range(numVal):
        trainX=[]
        trainY=[]
        testX=[]
        testY=[]
        random.shuffle(indexList)
        # 将数据分为测试集和训练集
        for j in range(m):
            if 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)
        for k in range(30):
            matTestX=mat(testX)
            matTrainX=mat(trainX)
            meanTrain=mean(matTrainX)
            varTrain=var(matTrainX,0)
            matTeatX=(matTestX-meanTrain)/varTrain
            yEst=matTestX*mat(wMat[k,:]).T+mean(trainY)
            errorMat[i,k]=rssError(yEst.T.A,array(testY))
    meanErrors=mean(errorMat,0)
    minMean=float(min(meanErrors))
    bestWeights=wMat[nonzero(meanErrors==minMean)]
    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))

crossValidation(lgX,lgY,10)

得到的输出结果为：