机器学习笔记9-各种线性回归

最新推荐文章于 2023-04-26 11:12:38 发布

Tobesix

最新推荐文章于 2023-04-26 11:12:38 发布

阅读量1.6k

点赞数 1

分类专栏：机器学习篇文章标签：机器学习

本文链接：https://blog.csdn.net/intricate_/article/details/105064763

版权

机器学习篇专栏收录该内容

21 篇文章 7 订阅

订阅专栏

本章介绍多种线性回归：

1.普通线性回归

2.局部加权线性回归

3.岭回归

4.逐步线性回归

1.普通线性回归(LR)

对于给定的x，预测结果y等于：

$y_{i} = x_{i}^{T}\omega$

平方误差损失函数：

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

使用最小二乘法可以得到最佳的w:

$\hat{\omega }=(X^{T}X)^{-1}X^{T}y$

缺点:线性回归有可能出现欠拟合现象，因为它求的时具有最小均方方差的无偏估计。

2.局部加权线性回归(LWLR)

我们可以给待预测点附近的每一个点赋予一定的权重，然后在这个子集上基于最小方差进行普通的回归。

该算法解出回归系数w的形式如下：

$\hat{\omega }=(X^{T}WX)^{-1}X^{T}Wy$

其中W时一个矩阵，用来给每一个数据赋予权重。

局部加权线性回归使用"核"（与支持向量机中的核类似）来对附近的点赋予更高的权重。核的类型可以自由选择，最常用的核就是高斯核，高斯核对应权重如下：

$w{(i,i)}=exp\left ( \frac{\left | x^{(i)}-x \right |}{-2k^{2}} \right )$

缺点：对于每一个数据点，都要运行整个数据集。

3.岭回归

如果数据的特征比样本点还多时，以上两个方法就失效了，因为在计算 $(X^{T}X)^{-1}$ 的时候会出错。

非满秩矩阵在求逆时会出现问题，统计学家引入了岭回归的概念。

岭回归简单来说就是在矩阵 $X^{T}X$ 上加上一个 $\lambda I$ 从而使得矩阵非奇异，其中矩阵 $I$ 是一个m*m的矩阵，对角线上元素全为1，其他元素全为0。而 $\lambda$ 是一个用户定义的数值。在这种情况下，回归系数的计算公式变为：

$\hat{\omega }=(X^{T}X+\lambda I)^{-1}X^{T}y$

岭回归最先用来处理特征数多于样本数的情况，现在也用于在估计中加入偏差，从而得到更好的估计，这里引入 $\lambda$ 来限制了所有w之和，通过引入该惩罚项，能够减少不重要的参数，类似于正则化参数。

通过岭回归的回归系数变化图，在最左边即λ最小时，可以得到所有系数的原始值(与普通线性回归相同) 而在右边，系数全部缩减成0 ;在中间部分的某值将可以取得最好的预测效果，另外，可以根据图中系数大小来判断系数对预测结果的影响力。

4.前向逐步回归

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

逐步线性回归通过多次迭代得到的图像和岭回归得到的图像类似(关于Y轴对称)，和岭回归一样逐步线性回归的优点在于它可以帮助人们理解现有的模型并做出改进，构建模型后，可以运行该算法找出重要的特征，这样就有可能及时停止对那些不重要特征的收集。当应用缩减方法（如逐步线性回归和岭回归）时，模型也就增加了偏差，与此同时减少了模型的方差。

Python代码：

import numpy as np
import matplotlib.pyplot as plt

# 读取文件，文件的第一列值为1.0，代表x0（偏移量），第二列为横坐标x，第三列为纵坐标y
def loadDataSet(fileName):
    fr = open(fileName)
    dataMat = []; labelMat = []
    for line in fr.readlines():
        lineArr =[float(example) for example in line.strip().split('\t')]
        dataMat.append(lineArr[:-1])
        labelMat.append(float(lineArr[-1]))
    return np.array(dataMat), np.array(labelMat)

'''
普通线性回归-最小二乘法
求解回归系数 w = ((X^T)X)^(-1)(X^T)y
@param xArr：数据集（包括x0=1）
@param yArr：目标值
'''
def standRegres(xArr,yArr):
    xMat = np.mat(xArr); yMat = np.mat(yArr).T
    xTx = xMat.T*xMat
    if np.linalg.det(xTx) == 0.0:  # NumPy线性代数库linalg中的函数det()计算行列式，如果Mat.T*xMat的行列式为0，则不存在X的逆，将会出现错误
        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 = np.mat(xArr)
# yMat = np.mat(yArr)
# 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() # 绘制拟合曲线
# print(np.corrcoef(yHat.T,yMat)) # 计算相关系数



'''
局部加权线性回归函数
@param testPoint：待预测点
@param xArr：数据集
@param yArr：目标值
@param k：高斯核中参数k，k越大，用于训练回归模型的数据越多
注意：局部加权线性回归函数的缺点在于，对于每一个数据点，都要运行整个数据集
'''
def lwlr(testPoint,xArr,yArr,k=1.0):
    xMat = np.mat(xArr); yMat = np.mat(yArr).T
    m = xMat.shape[0]
    weights = np.mat(np.eye((m)))  # 创建对角阵，主对角线均为1
    for j in range(m):       # 权重值大小以指数级衰减
        diffMat = testPoint - xMat[j,:]
        weights[j,j] = np.exp(diffMat*diffMat.T/(-2.0*k**2))
    xTx = xMat.T * (weights * xMat)
    if np.linalg.det(xTx) == 0.0:
        print("This matrix is singular, cannot do inverse")
        return
    ws = xTx.I * (xMat.T * (weights * yMat)) # 回归系数w
    return testPoint * ws  # 返回预测值

# 循环计算每一个数据点的预测值
# 缺点：对于每一个数据点，都要运行整个数据集
def lwlrTest(testArr,xArr,yArr,k=1.0):
    m = testArr.shape[0]
    yHat = np.zeros(m)
    for i in range(m):
        yHat[i] = lwlr(testArr[i],xArr,yArr,k)
    return yHat


# 测试局部加权线性回归函数
# xArr,yArr = loadDataSet('ex0.txt')
# yHat = lwlrTest(xArr, xArr, yArr, 0.003) # k=0.003，会出现过拟合
# xMat = np.mat(xArr)
# srtInd = xMat[:,1].argsort(0)  # 索引按照升序排列
# xSort = xMat[srtInd][:,0,:]  # 将xMat按照升序排列
# fig = plt.figure()
# ax = fig.add_subplot(111)
# ax.plot(xSort[:,1], yHat[srtInd]) # 绘制回归直线
# ax.scatter(xMat[:,1].flatten().A[0], np.mat(yArr).T.flatten().A[0], s=2, c='red')
# plt.show()
# print(np.corrcoef(yHat.T, yArr))


'''
残差平方和rss
@param yArr：目标值
@param yHatArr：预测值
'''
def rssError(yArr,yHatArr): # yArr和yHatArr都需要是数组
    return ((yArr-yHatArr)**2).sum()


'''
岭回归（解决特征比样本点还多的问题（n>m），防止不存在矩阵X的逆。现在也用于在估计中加入偏差，从而得到更好的估计）
@param xMat；数据集
@param yMat：目标值
@param lam：λ的值，用于限制所有w之和，通过引入该惩罚项，减少不重要的参数，在统计学中叫做缩减
'''
def ridgeRegres(xMat, yMat, lam=0.2):
    xTx = xMat.T * xMat
    denom = xTx + np.eye(xMat.shape[1]) * lam
    if np.linalg.det(denom) == 0.0:  # 还需要检查行列式，防止lam设置为0
        print("This matrix is singular, cannot do inverse")
        return
    ws = denom.I * (xMat.T * yMat)
    return ws

def ridgeTest(xArr, yArr):
    xMat = np.mat(xArr);
    yMat = np.mat(yArr).T
    yMean = np.mean(yMat, 0)  # 计算每一列的均值
    yMat = yMat - yMean
    # 对特征矩阵X归一化
    xMeans = np.mean(xMat, 0)  # 计算每一列的均值
    xVar = np.var(xMat, 0)  # 按照列计算方差
    xMat = (xMat - xMeans) / xVar  # 归一化，广播
    numTestPts = 30  # 计算30个不同的λ下的权值
    wMat = np.zeros((numTestPts, xMat.shape[1]))  # 30个λ下的权值组成的权值矩阵
    for i in range(numTestPts):
        ws = ridgeRegres(xMat, yMat, np.exp(i - 10))
        wMat[i, :] = ws.T
    return wMat  # 返回权值矩阵

'''
在最左边即λ最小时，可以得到所有系数的原始值(与普通线性回归相同) 
而在右边，系数全部缩减成0 ;在中间部分的某值将可以取得最好的预测效果。
'''
# 测试岭回归函数
# abX,abY = loadDataSet('abalone.txt')
# ridgeWeights = ridgeTest(abX,abY)
# print(ridgeWeights)  # 输出权值矩阵
# fig = plt.figure()
# ax = fig.add_subplot(111)
# ax.plot(ridgeWeights)  # 横坐标为λ，纵坐标(一列8个点)为权值w[0]~w[7]
# plt.show()

# 按照列进行归一化（均值为0，方差为1）
def regularize(xMat):
    inMat = xMat.copy()
    inMeans = np.mean(inMat,0)   # 按列求均值
    inVar = np.var(inMat,0)  # 按列求方差
    inMat = (inMat - inMeans)/inVar
    return inMat

'''
前向逐步回归：容易找出重要特征，即使停止对那些不重要特征的收集。
如果用于测试，可以使用类似于10折交叉验证的方法比较这些模型，选择rss最小的模型
@param xArr：数据集
@param yArr：目标值
@param eps：步长
@param numIt：迭代次数
'''
def stageWise(xArr,yArr,eps=0.01,numIt=100):
    xMat = np.mat(xArr); yMat = np.mat(yArr).T
    yMean = np.mean(yMat,0)
    yMat = yMat - yMean     # 也可以归一化y，但会得到更小的相关系数
    xMat = regularize(xMat)
    m,n = xMat.shape
    returnMat = np.zeros((numIt,n)) # 记录每次迭代的权值向量，构成一个矩阵 (迭代次数，特征数)
    ws = np.zeros((n,1)); wsMax = ws.copy()  # 权重初始化
    for i in range(numIt):
        lowestError = np.inf  # 设置当前最小误差为+∞
        for j in range(n):  # 对每个特征
            for sign in [-1,1]:  # 增大或减小
                wsTest = ws.copy()
                wsTest[j] += eps*sign # 按照一定步长，改变系数，得到一个新的W
                yTest = xMat*wsTest
                rssE = rssError(yMat.A,yTest.A) # 新W下的误差
                if rssE < lowestError:  # 如果误差小于当前误差，更新误差值
                    lowestError = rssE
                    wsMax = wsTest
        ws = wsMax.copy()
        returnMat[i,:] = ws.T
    return returnMat  # 返回迭代矩阵

# 测试前向逐步回归方法
# xArr,yArr = loadDataSet('abalone.txt')
# returnMat5000 = stageWise(xArr,yArr,0.001,5000) 
# print(returnMat5000[-1,:])
# fig = plt.figure()  # 绘图显示
# ax = fig.add_subplot(111)
# ax.plot(returnMat5000) # 横坐标为迭代次数，纵坐标(一列8个点)为权值w[0]~w[7]
# plt.show()
# # 与最小二乘法比较，发现和5000次迭代，步长0.001的结果类似
# xMat = np.mat(xArr)
# yMat = np.mat(yArr).T
# xMat = regularize(xMat)
# yMat = yMat - np.mean(yMat,0)
# weights = standRegres(xMat,yMat.T)
# print(weights.T)