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 ShoppingAPI收集数据。

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

1.7 总结

        当数据的样本数比特征数还少时候,矩阵x^{T}x的逆不能直接计算。即便当样本数比特征数多时,x^{T}x的逆仍有可能无法直接计算,这是因为特征有可能高度相关。这时可以考虑使用岭回归,因为当x^{T}x的逆不能计算时,它仍保证能求得回归参数。

        岭回归是缩减法的一种,相当于对回归系数的大小施加了限制。另一种很好的缩减法是lasso。Lasso难以求解,但可以使用计算简便的逐步线性回归方法来求得近似结果。

        缩减法还可以看做是对一个模型增加偏差的同时减少方差。偏差方差折中是一个重要的概念,可以帮助我们理解现有模型并做出改进,从而得到更好的模型。

Logo

有“AI”的1024 = 2048,欢迎大家加入2048 AI社区

更多推荐