预测数值型数据：回归

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

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

缺点：对非线性的数据拟合不好

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

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

回归的一般方法：

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

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

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

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

（5）测试算法：使用R2或者预测值和数据的拟合度，来分析模型的效果。

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

输入数据存放在矩阵X中,找出使误差最小的w。这里的误差是指预测y值和真实y值之间的差值，使用该误差的简单累加将使得正差值和负差值相互抵消，所以我们采用平方误差。

标准回归函数和数据导入函数：

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, cannot do inverse")
        return
    ws = xTx.I * (xMat.T*yMat)
    return ws

结果：

import regression
from numpy import *

xArr, yArr = regression.loadDataSet('ex0.txt')
print('xArr=', xArr[:2])
ws = regression.standRegres(xArr, yArr)
print('回归系数=', ws)
xMat = mat(xArr)
yMat = mat(yArr)
yHat = xMat * ws
print("预测值=", yHat)

import matplotlib.pyplot as plt

fig = plt.figure()
ax = plt.subplot(111)
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.show()

判断模型好坏：计算预测值yHat序列和真实值y序列的匹配程度，那就是计算这两个序列的相关系数。

yHat = xMat * ws
cor = corrcoef(yHat.T, yMat)
print('相关系数：', cor)

1.2 局部加权线性回归

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

其中的一个方法是局部加权线性回归（Locally Weighted Linear Regression，为LWLR）。在该算法中，我们给待预测点附近的每个点赋予一定的权重；在这个子集上基于最小均方差来进行普通的回归。与kNN一样，这种算法每次预测均需要事先选取出对应的数据子集。该算法解出回归系数w的形式如下：

其中w是一个矩阵，用来给每个数据点赋予权重。LWLR使用“核”（与支持向量机中的核类似）来对附近的点赋予更高的权重核的类型可以自由选择，最常用的核就是高斯核，高斯核对应的权重如下：

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

局部加权线性回归：

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')
# yArr[0]

# 数据集中的所有点的估计
yHat = np.zeros((len(xArr), 1))
for i in range(len(xArr)):
    # yHat[i] =  regression.lwlr(xArr[0], xArr, yArr, 1.0)
    # yHat[i] = regression.lwlr(xArr[0], xArr, yArr, 0.001)
    yHat[i] = regression.lwlr(xArr[i], xArr, yArr, 0.003)

xMat = mat(xArr)
srtInd = xMat[:, 1].argsort(0).A.flatten()  # 确保索引是一维的
xSort = xMat[srtInd][:, 1]  # 选择第二列，并确保是一维的

# 然后用Matplotlib绘图：
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(xSort.A, yHat[srtInd], 'b')  # 确保 x 和 y 都是一维的
ax.scatter(xMat[:, 1].flatten().A[0], mat(yArr).T.flatten().A[0], s=2, c='red')
plt.show()

1.3 实例：预测鲍鱼的年龄

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)
print('yHat01=', yHat01)
yHat1 = regression.lwlrTest(abX[0:99], abX[0:99], abY[0:99], 1)
print('yHat1=', yHat1)
yHat10 = regression.lwlrTest(abX[0:99], abX[0:99], abY[0:99], 10)
print('yHat10', yHat10)

# 分析预测误差的大小，可以用函数rssError()计算出这一指标：
err01 = regression.rssError(abY[0:99], yHat01.T)
print('err01=', err01)
err1 = regression.rssError(abY[0:99], yHat1.T)
print('err1=', err1)
err10 = regression.rssError(abY[0:99], yHat10.T)
print('err10=', err10)

"""可以看到，使用较小的核将得到较低的误差。那么，为什么不在所有数据
集上都使用最小的核呢？这是因为使用最小的核将造成过拟合，对新数据
不一定能达到最好的预测效果。下面就来看看它们在新数据上的表现："""
yHat01 = regression.lwlrTest(abX[100:199], abX[0:99], abY[0:99], 0.1)
new_err01 = regression.rssError(abY[100:199], yHat01.T)
print('new_err01=', new_err01)

yHat1 = regression.lwlrTest(abX[100:199], abX[0:99], abY[0:99], 1)
new_err1 = regression.rssError(abY[100:199], yHat1.T)
print('new_err0=', new_err1)

yHat10 = regression.lwlrTest(abX[100:199], abX[0:99], abY[0:99], 10)
new_err10 = regression.rssError(abY[100:199], yHat10.T)
print('new_err10=', new_err10)

"""从上述结果可以看到，在上面的三个参数中，核大小等于10时的测试误差
最小，但它在训练集上的误差却是最大的。接下来再来和简单的线性回归
做个比较："""
ws = regression.standRegres(abX[0:99], abY[0:99])
yHat = mat(abX[100:199])*ws
com = regression.rssError(abY[100:199], yHat.T.A)
print('com=', com)

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

如果输入的矩阵不是满秩矩阵，其在求逆时会出现问题，因此需要引入岭回归进行缩减。其次是lasso方法，效果很好但计算复杂。前向逐步回归，效果和lasso方法差不多且更容易实现。

1.4.1 岭回归

岭回归就是在矩阵 $x^{^{T}}x$ 上加一个λI从而使得矩阵非奇异，进而能对 $x^{^{T}}x$ + λI求逆。其中矩阵I是一个m×m的单位矩阵，对角线上元素全为1，其他元素全为0。而λ是一个用户定义的数值，回归系数的计算公式将变成：

岭回归：

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

结果：

abX, abY = regression.loadDataSet('abalone.txt')
ridgeWeights = regression.ridgeTest(abX, abY)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(ridgeWeights)
plt.show()

1.4.2 lasso

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

与岭回归相似，lasso也对回归系数做了限定，对应的约束条件：

1.4.3 前向逐步回归

它属于一种贪心算法，即每一步都尽可能减少误差。一开始，所有的权重都设为1，

然后每一步所做的决策是对某个权重增加或减少一个很小的值。

该算法的伪代码如下所示：

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

在每轮迭代过程中：

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

对每个特征：

增大或缩小：

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

计算新W下的误差

如果误差Error小于当前最小误差lowestError：设置Wbest等于当前的W

将W设置为新的Wbest

前向逐步线性回归：

def regularize(xMat):
    # 假设 regularize 函数对数据进行标准化处理
    inMeans = np.mean(xMat, 0)
    inVar = np.var(xMat, 0)
    xMat = (xMat - inMeans) / inVar
    return xMat

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

结果：

xArr, yArr = regression.loadDataSet('abalone.txt')
sta = regression.stageWise(xArr, yArr, 0.01, 200)
print('sta:', sta)
# 更小的步长和更多的步数
sta1 = regression.stageWise(xArr, yArr, 0.001, 5000)
print('sta1:', sta1)
# 最小二乘法
xMat=mat(xArr)
yMat=mat(yArr).T
xMat=regression.regularize(xMat)
yM = mean(yMat, 0)
yMat = yMat - yM
weights = regression.standRegres(xMat, yMat.T)
print(weights.T)

1.5 权衡偏差与方差

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

偏差方差折中与测试误差及训练误差的关系。上面的曲线就是测试误差，在中间部分最低。为了做出最好的预测，我们应该调整模型复杂度来达到测试误差的最小值。

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

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

（1）收集数据：用Google Shopping的API收集数据。

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

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

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

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

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

1.6.1 收集数据：使用Google购物的API

from time import sleep
import json
import urllib.request

def searchForSet(retX, retY, setNum, yr, numPce, origPrc):
    sleep(10)
    myAPIstr = 'get from code.google.com'  # 请确保你有一个有效的 API 密钥
    searchURL = 'https://www.googleapis.com/shopping/search/v1/public/products?key=%s&country=US&q=lego+%d&alt=json' % (myAPIstr, setNum)
    try:
        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 = 1
                else:
                    newFlag = 0
                listOfInv = 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)
    except Exception as e:
        print('Failed to retrieve data:', e)
        print('Please check the URL and your network connection, and try again.')

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(regression.setDataCollect(lgX, lgY))

1.6.2 训练算法：建立模型

交叉验证测试岭回归：

def crossValidation(xArr, yArr, numVal=10):
    m = len(yArr)
    indexList = 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, 0)
        varTrain = var(matTrainX, 0)
        matTestX = (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))

结果：

laX1 = mat(ones((shape(lgX))))
ws = regression.standRegres(laX1, lgY)
print(ws)

print(regression.crossValidation(lgX, lgY, 10))
# 缩减
print(regression.ridgeTest(lgX, lgY))