《机器学习实战》学习笔记———预测数值型数据：回归

---

引言

前面所学习的基本都是分类问题，分类问题的目标变量是标称型数据，或者离散型数据。而回归的目标变量为连续型，也即是回归对连续型变量做出预测，最直接的办法是依据输入写出一个目标值的计算公式，这样，对于给定的输入，利用该公式可以计算出相应的预测输出。这个公式称为回归方程，而求回归方程显然就是求该方程的回归系数，而一旦有了这些回归系数，再给定输入，就可以将这些回归系数乘以输入值，就得到了预测值。

---

用线性回归拟合

回归方程（regression equation），回归系数（regression weights），求回归系数的过程就是回归。说到回归，一般都是指线性回归（linear regression），还存在非线性回归模型。

假定输入数据存放在矩阵$X$中，而回归系数存放在向量$w$中，那么对于给定的数据$X1$，预测结果将会通过$Y1=X_{1}^{T}w$。问题是，自己有一些$x$和对应的$y$，如何找到$w$呢？常用的方法是找到使误差（一般采用平方误差，否则正负误差会相互抵消）最小的$w$。平方误差可写作：

$\sum_{i=1}^{m}\left ( y_{i}-x_{i}^{T}w \right )^{2}$

用矩阵表示还可以写成$\left ( Y-Xw \right )^{T}\left ( Y-Xw \right )$，如果对$w$求导，得到$X^{T}\left ( Y-Xw \right )$，令其等于零，解出$w$如下（这里的w^w^是ww的一个最优解、最佳估计）：

$\hat{w}=\left ( X^{T}X \right )^{-1}X^{T}Y$
公式中的$\left ( X^{T}X \right )^{-1}$是对矩阵求逆，需要在代码中判断逆矩阵是否存在，可以计算行列式，如果行列式为零，那么计算逆矩阵的时候会出现错误。这种求解最佳$w$的方法称为OLS“普通最小二乘法”（ordinary least squares）。

求解完$w$，可以计算$\hat{Y}=X\hat{w}$，然后可以来计算预测值和真实值的相关性。

实例分析

对于下图的散点图

可用如下代码完成回归过程

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

再运行如下代码

xArr,yArr=loadDataSet('ex0.txt')
ws = standRegres(xArr,yArr)
print(ws)

可知回归系数为

将其图像化出即可得

xMat=mat(xArr)
yMat=mat(yArr)
yHat = xMat*ws
fig=plt.figure()
ax = fig.add_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()

最佳拟合直线方法将数据视为直线进行建模。但数据中似乎还存在其他的潜在模式。可根据数据来局部调整预测。

---

局部加权线性回归

线性回归会出现欠拟合现象，因为它求的是最小均方误差的无偏估计。可以在估计中引入一些偏差，从而降低预测的均方误差。其中一个方法是局部加权线性回归（Locally Weighted Linear Regression，LWLR），该算法中给待测点附近的每个点赋予一定的权重，然后在这个子集上基于最小均方差来进行普通的回归。与kNN一样，此算法每次预测均需事先选取出对应的数据子集。该算法解出的回归系数的形式如下：

$\hat{w}=\left ( X^{T}WX \right )^{-1}X^{T}WY$
其中$w$是一个矩阵，用来给每个数据点赋予权重。LWLR使用“核”（与支持向量机中的核类似）来对附近的点赋予更高的权重。核的类型可以自由选择，最常用的是高斯核，高斯核对应的权重如下：

$w\left ( i,i \right )=exp\left ( \frac{\left | x^{\left ( i \right )}-x \right |}{-2k^{2}} \right )$
这样就构建了一个只含对角元素的权重矩阵$W$，并且点$x$与$w\left ( i,i \right )$越近，$w\left ( i,i \right )$将会越大，公式中的$k$需要用户指定，它决定了对附近的点赋予多大的权重。

要注意区分权重$W$和回归系数$w$，与kNN一样，该加权模型认为样本点距离越近，越可能符合同一个线性模型。

实例分析

针对上面的数据，运用如下代码

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):                      #next 2 lines create weights matrix
        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):  #loops over all the data points and applies lwlr to each one
    m = shape(testArr)[0]
    yHat = zeros(m)
    for i in range(m):
        yHat[i] = lwlr(testArr[i],xArr,yArr,k)
    return yHat

运用取不同的$k$值，得到不同的回归效果

xArr,yArr=loadDataSet('ex0.txt')
ws = standRegres(xArr,yArr)
yHat=lwlrTest(xArr,xArr,yArr,0.003)#or1,0.01
xMat=mat(xArr)
srtInd=xMat[:,1].argsort(0)
xSort=xMat[srtInd][:,0,:]
fig=plt.figure()
ax = fig.add_subplot(111)
ax.plot(xSort[:,1],yHat[srtInd])
ax.scatter(xMat[:,1].flatten().A[0],mat(yArr).T[:,0].flatten().A[0],s=2,c='red')
plt.show()

$k$=1.0时

$k$=0.01时

$k$=0.03时

$k$的选取有事能挖出数据的潜在规律，有时考虑太多的噪声，进而导致了过拟合现象，如$k=0.03$时。

---

实例：预测鲍鱼的年龄

将回归用于真实数据，有一份来自UCI数据聚合的数据，记录了鲍鱼（一种介壳类水生动物）的年龄。鲍鱼年龄可以从鲍鱼壳的层数推算得到。

加入如下代码

def rssError(yArr,yHatArr): 
    return ((yArr-yHatArr)**2).sum()

再添加：

abX,abY=loadDataSet('abalone.txt')
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)
rssError(abY[0:99],yHat01.T)
rssError(abY[0:99],yHat1.T)
rssError(abY[0:99],yHat10.T)

可得到三种不同 k k <script type="math/tex" id="MathJax-Element-70">k</script>值得误差

可以看到。使用较小的核将得到较低的误差。但有时不用最小的核的原因是会造成过拟合，对新数据不一定能达到最好的预测效果。比如将上述预测用于在新数据上

abX,abY=loadDataSet('abalone.txt')
yHat01=lwlrTest(abX[100:199],abX[0:99],abY[0:99],0.1)
print(rssError(abY[100:199],yHat01.T))
yHat1=lwlrTest(abX[100:199],abX[0:99],abY[0:99],1)
print(rssError(abY[100:199],yHat1.T))
yHat10=lwlrTest(abX[100:199],abX[0:99],abY[0:99],10)
print(rssError(abY[100:199],yHat10.T))

可看出，在上面三个参数，核大于10时的测试误差最小，但在训练集上的误差是最大的。

---

代码

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

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):                      #next 2 lines create weights matrix
        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):  #loops over all the data points and applies lwlr to each one
    m = shape(testArr)[0]
    yHat = zeros(m)
    for i in range(m):
        yHat[i] = lwlr(testArr[i],xArr,yArr,k)
    return yHat

def lwlrTestPlot(xArr,yArr,k=1.0):  #same thing as lwlrTest except it sorts X first
    yHat = zeros(shape(yArr))       #easier for plotting
    xCopy = mat(xArr)
    xCopy.sort(0)
    for i in range(shape(xArr)[0]):
        yHat[i] = lwlr(xCopy[i],xArr,yArr,k)
    return yHat,xCopy

def rssError(yArr,yHatArr): #yArr and yHatArr both need to be arrays
    return ((yArr-yHatArr)**2).sum()

---

参考文献

《机器学习实战》
https://blog.csdn.net/namelessml/article/details/52554754

---