回归分析及实际案例：预测鲍鱼年龄

heda3

已于 2022-11-07 20:22:08 修改

阅读量1.1w

点赞数 20

分类专栏：机器学习文章标签：回归监督学习

于 2020-01-28 22:42:58 首次发布

本文链接：https://blog.csdn.net/heda3/article/details/104102732

版权

机器学习专栏收录该内容

28 篇文章 5 订阅

订阅专栏

上一篇文章：线性回归（Linear regression）算法

引入：

1、线性回归：

算法的优点：

结果易于理解，计算不复杂

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

目标：平方误差和最小

求解（对参数w求导等于0）的回归系数：

模型预测：

$f(x)=W^{T}X$

"""
函数说明:标准回归

Parameters:
    xArr - 特征矩阵
    yArr -响应值


Returns:
     ws- 回归系数

Author:
    heda3
Blog:
    https://blog.csdn.net/heda3
Modify:
    2020-01-10
"""   
def standRegres(xArr,yArr):
    xMat = mat(xArr); yMat = mat(yArr).T
    xTx = xMat.T*xMat#计算xTx
    if linalg.det(xTx) == 0.0:#判断行列式是否为0
        print("This matrix is singular, cannot do inverse")
        return
    ws = xTx.I * (xMat.T*yMat)#计算回归系数
    return ws

注意行列式不为零才可以计算逆矩阵

#加载测试数据

from numpy import *

def loadDataSet(fileName):      #general function to parse tab -delimited floats
    numFeat = len(open(fileName).readline().split('\t')) - 1 #get number of fields 
    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

数据描述：

特征：2，响应:数值型

使用线性回归预测数据并绘制散点图+计算模型的效果（相关系数计算）

 ##测试线性回归
 xArr,yArr=loadDataSet('ex0.txt')
 xArr[0:2]
 ws=standRegres(xArr,yArr)#计算回归系数
 xMat=mat(xArr)
 yMat=mat(yArr)
 #绘散点图
 import matplotlib.pyplot as plt
 fig=plt.figure()
 ax=fig.add_subplot(111)
 ax.scatter(xMat[:,1].flatten().A[0],yMat.T[:,0].flatten().A[0])#.A转变为数组
 xCopy=xMat.copy()
 xCopy.sort(0)#维度，行排序
 yHat=xCopy*ws
 ax.plot(xCopy[:,1],yHat)
 plt.show
 #求预测值和真实值的相关系数
 yHat1=xMat*ws
 corrcoef(yHat1.T,yMat)

2、局部加权线性回归

Locally Weighted Linear Regression， LWLR
线性回归存在的问题是：出现欠拟合

解决方法：在估计中引入偏差，从而降低预测的均方误差

也即是通过在每一小段进行拟合，以逼近真实的数据

局部加权线性回归相比普通线性回归的问题是：每次必须在整个数据集上运行，也即是必须要保存所有的训练数据

思路：在待预测点附近的每个点赋予一定的权重

求解的回归系数：

对比线性回归：

其中的W矩阵用于给每个数据点赋予权重

权重W的选择，通过使用不同的核：

例如高斯核

需要调节的参数：一个 k

#对单点估计

"""
函数说明:局部加权线性回归

Parameters:
    testPoint x空间的任意一点
    xArr - 特征矩阵
    yArr -响应值
    k 和权重有关，当k越小则使用的越少的局部数据集进行训练，k=1相当于标准线性回归

Returns:
    某个点的预测结果  testPoint * ws - 

Author:
    heda3
Blog:
    https://blog.csdn.net/heda3
Modify:
    2020-01-10
"""   
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,:]     #x-x[i]
        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

#对多点（数据集）估计

"""
函数说明:为数据集中的每个点调用lwlr
  
Parameters:
    testArr 测试的数据集
    xArr - 特征矩阵
    yArr -响应值
    k

Returns:
     testPoint * ws - 

Author:
    heda3
Blog:
    https://blog.csdn.net/heda3
Modify:
    2020-01-10
"""  
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

对数据进行预测+绘制散点图

 ##测试局部加权线性回归
 xArr,yArr=loadDataSet('ex0.txt')
 #对单点估计
 yArr[0]
 lwlr(xArr[0],xArr,yArr,1.0)
 lwlr(xArr[0],xArr,yArr,0.001)
 #数据集中所有点的估计
 yHat=lwlrTest(xArr,xArr,yArr,0.01)
 #绘制散点图
 xMat=mat(xArr)
 srtInd=xMat[:,1].argsort(0)
 xSort=xMat[srtInd][:,0,:] #等价于xMat[srtInd.flatten().A[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')#.A转变为数组
 plt.show()

设置参数k值为0.01

问题：

数据的特征比样本点多时，计算矩阵的(XTX)的逆出错，也即是输入数据（特征）矩阵不是满秩的矩阵

或者是数据特征之间是高度相关时也不能计算

如何减少特征数，如何减少不重要的特征？

解决方法：缩减系数

1）岭回归

2）lasso

3）LAR

4）PCA回归

5）子集选择

1）岭回归

$X^{T}X$ 的基础上加 $\lambda I$ 使得矩阵非奇异，从而能对 $X^{T}X+\lambda I$ 求逆

求解的回归系数：

对比标准线性回归：

选择的参数： $\lambda$ 可通过交叉验证确定

注意：数据需要先标准化处理（X-mean）/var

"""
函数说明：岭回归
  
Parameters:
    xMat- 数据的特征 假设有n 样本个数有m
    yMat- 响应值
    lam-- 调节的参数
Returns:
     ws-- 计算出的回归系数
Author:
    heda3
Blog:
    https://blog.csdn.net/heda3
Modify:
    2020-01-10
"""  
def ridgeRegres(xMat,yMat,lam=0.2):
    xTx = xMat.T*xMat#2*2   n*n
    denom = xTx + eye(shape(xMat)[1])*lam#n*n
    if linalg.det(denom) == 0.0:
        print("This matrix is singular, cannot do inverse")
        return
    ws = denom.I * (xMat.T*yMat)
    return ws

"""
函数说明：岭回归参数lambda调节
  
Parameters:
    xArr- 特征矩阵
    yArr- 响应值

Returns:
    wMat - 返回一组w（维数和特征数对应）系数
Author:
    heda3
Blog:
    https://blog.csdn.net/heda3
Modify:
    2020-01-10
"""  
def ridgeTest(xArr,yArr):
    xMat = mat(xArr); yMat=mat(yArr).T
    #数据标准化处理
    yMean = mean(yMat,0)
    yMat = yMat - yMean     #to eliminate X0 take mean off of Y
    #regularize X's
    xMeans = mean(xMat,0)   #calc mean then subtract it off
    xVar = var(xMat,0)      #calc variance of Xi then divide by it
    xMat = (xMat - xMeans)/xVar
    numTestPts = 30#设置lambda参数迭代次数
    wMat = zeros((numTestPts,shape(xMat)[1]))#30*2的矩阵
    for i in range(numTestPts):
        ws = ridgeRegres(xMat,yMat,exp(i-10))
        wMat[i,:]=ws.T
    return wMat

2)lasso

对回归系数的约束：

3）前向逐步回归

每一步都尽可能的减少误差，通过设置初始权重为1，每一步所做的决策是对某个权重增加或减少一个很小的值

算法步骤：

通过多次迭代后得到趋于稳定的回归参数！

可调节的参数：步长和迭代次数

"""
函数说明：前向逐步线性回归
  
Parameters:
    xArr- 特征矩阵
    yArr- 响应值
    eps=0.01  每次迭代需要调整的步长
    numIt=100 迭代次数
Returns:
    returnMat - 返回一组w（维数和特征数对应）系数
Author:
    heda3
Blog:
    https://blog.csdn.net/heda3
Modify:
    2020-01-11
"""  
def stageWise(xArr,yArr,eps=0.01,numIt=100):
    xMat = mat(xArr); yMat=mat(yArr).T
    #数据的标准化
    yMean = mean(yMat,0)
    yMat = yMat - yMean     #can also regularize ys but will get smaller coef
    xMat = regularize(xMat)
    m,n=shape(xMat)
    returnMat = zeros((numIt,n)) #testing code remove
    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):#n个特征也即是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()#迭代numIt次找到最小误差的回归系数
        returnMat[i,:]=ws.T
    return returnMat

案例1：预测鲍鱼年龄

问题出发点是：鲍鱼的年龄是通过贝壳的年轮计数确定，此方法耗时费力，如何依据其它的一些参数来推测鲍鱼的年龄？

数据集来源：UCI Machine Learning Repository: Abalone Data Set

数据集的描述：

上述的数据集给出的8个特征属性

Name / Data Type / Measurement Unit / Description
-----------------------------
Sex / nominal / -- / M, F, and I (infant)
Length / continuous / mm / Longest shell measurement
Diameter / continuous / mm / perpendicular to length
Height / continuous / mm / with meat in shell
Whole weight / continuous / grams / whole abalone
Shucked weight / continuous / grams / weight of meat
Viscera weight / continuous / grams / gut weight (after bleeding)
Shell weight / continuous / grams / after being dried

响应值：

构建模型预测

1）使用标准线性回归模型预测鲍鱼年龄

 ##使用标准线性回归模型
 ws=standRegres(abX[0:99],abY[0:99])
 yHat=mat(abX[100:199])*ws
 rssError(abY[100:199],yHat.T.A)

新数据上表现平方误差=518.6

2）使用局部线性加权线性回归模型预测鲍鱼年龄

 ##使用局部加权线性回归
 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)
 #测试集表现
 yHat01=lwlrTest(abX[100:199],abX[100:199],abY[100:199],0.1)
 yHat1=lwlrTest(abX[100:199],abX[100:199],abY[100:199],1)
 yHat10=lwlrTest(abX[100:199],abX[100:199],abY[100:199],10)
 x1=rssError(abY[100:199],yHat01.T)
 print(x1)
 x2=rssError(abY[100:199],yHat1.T)
 print(x2)
 x3=rssError(abY[100:199],yHat10.T)
 print(x3)

在不同的参数k下的效果，以及在训练集和测试集上的表现：

平方误差结果：

训练集下：

测试集下：

在未知数据上比较效果才可以很好的选择模型，10折交叉验证

3）岭回归模型预测鲍鱼年龄

 ##使用岭回归
 abX,abY=loadDataSet('abalone.txt')
 ridgeWeights=ridgeTest(abX,abY)#得到不同lambda下计算出的回归系数
 #绘制回归系数
 import matplotlib.pyplot as plt
 fig=plt.figure()
 ax=fig.add_subplot(111)
 ax.plot(ridgeWeights)
 plt.show()

横坐标为lambda值，y轴为各回归系数

通过上图看看出哪些变量对结果的预测具有影响力

为定量找到最佳参数值lambda，需要进行交叉验证获得误差最小的lambda

使用10折交叉验证计算得出最佳的岭回归系数，参与预测新的数据

在岭回归中要求数据要标准化再参与计算，那么在训练完成后新的数据如何进行预测？这个新的数据怎么利用训练的数据进行标准化？

解决方法是：利用在训练数据中得出的回归参数，通过变换实现变相的在新数据预测时的标准化

新的数据一般预测过程：

数据标准化：XT=（XTest-mean（XTrain））/Var(XTrain)

预测：Ytest=XT*Ws+mean(YTrain)

将上述的公式变换：

Ytest=(（XTest-mean（XTrain））/Var(XTrain))*Ws+mean(YTrain)

设UnReg=Ws/Var(XTrain)

constantTerm=-mean（XTrain）*Ws/Var(XTrain)+mean(YTrain)

则Ytest=XTest*UnReg+constantTerm（现在的新变换后的预测过程）

###交叉验证--岭回归
ridgeWs,ridgeunReg,ridgeConstantTerm=crossValidation(abX[0:99],abY[0:99],10)#目的是找出最佳的岭回归系数
##测试均方误差
####和标准线性回归的比较
xMat=mat(abX[100:199])
yMat=mat(abY[100:199]).T
ridgeyHat=xMat*ridgeunReg.T+ridgeConstantTerm#岭回归预测
rssError(abY[100:199],ridgeyHat.T.A)#误差计算

换种写法：

###交叉验证--岭回归
ridgeWs,ridgeunReg,ridgeConstantTerm=crossValidation(abX[0:99],abY[0:99],10)#目的是找出最佳的岭回归系数
##测试均方误差
####和标准线性回归的比较
xMat=mat(abX[100:199])
yMat=mat(abY[100:199])
ridgeyHat=xMat*ridgeunReg.T+ridgeConstantTerm#岭回归预测
rssError(yMat.A,ridgeyHat.T.A)#误差计算
xxx=yMat.A
yyy=ridgeyHat.T.A

"""
函数说明：交叉验证测试岭回归
  
Parameters:
    xArr - 特征
    yArr - 标签
    numVal=10 - 交叉验证的次数

Returns: 
    bestWeights 最佳的岭回归参数
    为了和标准线性回归比较
    unReg,constantTerm 数据标准化还原后的特征参数和常量参数
Author:
    heda3
Blog:
    https://blog.csdn.net/heda3
Modify:
    2020-01-28
""" 
def crossValidation(xArr,yArr,numVal=10):
    m = len(yArr)#样本点个数                           
    indexList = list(range(m))
    errorMat = zeros((numVal,30))#create error mat 30columns numVal rows
    for i in range(numVal):#交叉验证
        trainX=[]; trainY=[]
        testX = []; testY = []
        random.shuffle(indexList)#随机打乱样本索引
        #训练集和测试集的划分 90%训练  10%测试
        for j in range(m):#create training set based on first 90% of values in indexList
            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)    #30*特征数 get 30 weight vectors from ridge
        #30组回归系数
        for k in range(30):#loop over all of the ridge estimates
            matTestX = mat(testX); matTrainX=mat(trainX)
            meanTrain = mean(matTrainX,0)
            varTrain = var(matTrainX,0)
            
            matTestX = (matTestX-meanTrain)/varTrain #regularize test with training params
            yEst = matTestX * mat(wMat[k,:]).T + mean(trainY)#test ridge results and store
            errorMat[i,k]=rssError(yEst.T.A,array(testY))
            #print errorMat[i,k]
    #计算所有这些误差值的均值
    meanErrors = mean(errorMat,0)#errorMat为 10*30   30个岭回归参数  10次交叉验证  按照把轴向数据求平均  得到每列数据的平均值，也即是10折交叉验证的平均   calc avg performance of the different ridge weight vectors
    
    minMean = float(min(meanErrors))#哪个岭回归参数下的误差最小
    bestWeights = wMat[nonzero(meanErrors==minMean)]#找出误差最小的回归参数
    #can unregularize to get model
    #when we regularized we wrote Xreg = (x-meanX)/var(x)
    #we can now write in terms of x not Xreg:  x*w/var(x) - meanX/var(x) +meanY
    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)
    #标准化后数据还原
    constantTerm=-1*sum(multiply(meanX,unReg)) + mean(yMat)
    print("with constant term: ",constantTerm)
    return bestWeights,unReg,constantTerm

平方误差结果：

4）前向逐步回归模型预测鲍鱼年龄

 ##使用逐步回归
 xArr,yArr=loadDataSet('abalone.txt')
 returnMat1=stageWise(xArr,yArr,0.01,200)#出现来回震荡情况，原因是步长太大？原因是系数已经饱和需要调小系数
 #对比更小的步长
 returnMat2=stageWise(xArr,yArr,0.001,200)
 #绘制回归系数
 import matplotlib.pyplot as plt
 fig=plt.figure()
 ax=fig.add_subplot(111)
 ax.plot(returnMat2)
 plt.show()

参数：步长和迭代次数

基于上图可以较好的发现重要特征

定量的选择模型参数：使用类似交叉验证方法

实验结果：

##数据训练
returnMat2=stageWise(xArr,yArr,0.001,5000)
###预测
####数据标准化
xMat = mat(xArr); yMat=mat(yArr).T
xMean=mean(xMat,0)
VarX=var(xMat,0)
yMean = mean(yMat,0)

Xtest=(mat(xArr[100:199])-xMean)/VarX#标准化后的测试数据
yMat=mat(yArr[100:199])#实际标签数据

Ytest=Xtest*mat(returnMat2[4999,:]).T+yMean
rssError(yMat.A,Ytest.T.A)#.A转变为数组

均方误差：