Python实现最小角度回归

柔光的暖阳◎ 2022-03-20 10:59 155阅读 0赞

关于LARS的详细介绍，可以参考Bradley Efron, Trevor Hastie, Iain Johnstone and Robert Tibshirani等人的那篇经典论文Least Angle Regression，文末附有论文地址。本文仅涉及LARS的Python实现。

实现代码：

# 导入要用到的相关包和函数
    %matplotlib inline
    import numpy as np
    import pandas as pd
    from math import sqrt
    import matplotlib.pyplot as plt

# 导入要用到的红酒品质数据集
    wine_data = pd.read_csv('wine.csv')
    wine_data.head()

wine\_data的表结构如下所示：

![20190213110426806.PNG][]

# 对wine_data中的数据做标准化处理：
    wine_data = wine_data.apply(lambda x:(x - x.mean()) / x.std(),axis=0)
    wine_data.head()

wine\_data中的数据经标准化处理后，其表结构如下所示：

![20190213110807134.PNG][]

# 查看经标准化处理后的wine_data的统计信息
    wine_data.describe()

wine\_data描述性统计信息如下所示(部分)：

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQxMDgwODUw_size_16_color_FFFFFF_t_70][]

xNormalized = []  # 构造用于存放标准化处理后的属性数据的列表
    
    # 提取出wine_data中标准化处理后的标签数据并放入列表中
    labelNormalized = [float(label) for label in wine_data.iloc[:,-1].tolist()]  
    
    names = wine_data.columns.tolist()    # 提取出wine_data中所有属性的名称并放入列表中
    
    # 列表xNormalized中的每个元素对应着标准化处理后的wine_data数据中除去标签列的每一行
    for i in range(len(wine_data)):
        xNormalized.append(wine_data.iloc[i,0:-1].tolist()) 
       
    
    nrows = len(xNormalized)        # 总的样本数量
    ncols = len(xNormalized[0])     # 属性(特征)数量
    
    # 初始化相关参数
    beta = [0.0] * ncols  # beta列表为参数向量
    betaMat = []    # betaMat列表存放每次迭代后的beta参数列表
    
    nSteps = 350         # 迭代次数
    stepSize = 0.004     # 迭代步长
    nzList = []          # 存放属性对应的索引
    residualsMat = []    # residualsMat列表存放每次迭代后的residuals列表
    for i in range(nSteps):
        
        # 计算xNormalized中的每个样本对应的残差
        residuals = [0.0] * nrows
        for j in range(nrows):
            labelsHat = sum([xNormalized[j][k] * beta[k] for k in range(ncols)])
            residuals[j] = labelNormalized[j] - labelsHat
            
        residualsMat.append(list(residuals))
    
        # 分别计算每个属性与残差的关联值的平均值(两个变量的关联值等于标准化(减去均值后再除以标准差)值的乘积)
        corr = [0.0] * ncols
        for j in range(ncols):
            corr[j] = sum([xNormalized[k][j] * residuals[k] for k in range(nrows)]) / nrows
    
        iStar = 0
        corrStar = corr[0]
        
        # 找出绝对值最大的关联值及其对应的索引
        for j in range(1, (ncols)):
            if abs(corrStar) < abs(corr[j]):
                iStar = j
                corrStar = corr[j]
    
        # 利用属性与残差的最大关联值(绝对值最大)更新最大关联值对应的beta参数
        beta[iStar] += stepSize * corrStar / abs(corrStar)    # 关联值为正，beta参数对应增加，关联值为负，beta参数对应减小
        betaMat.append(list(beta))    # 将只有一个beta值被更新的beta参数列表存入列表betaMat中
    
        # 找出beta参数列表中所有不等于零的beta参数对应的索引
        nzBeta = [index for index in range(ncols) if beta[index] != 0.0]
        
        # 将在nzBeta列表中但不在nzList列表中的索引放入nzList列表中
        for q in nzBeta:
            if (q in nzList) == False:
                nzList.append(q)
    
    # 列表nameList中的属性已经按照属性对预测目标的重要性排序
    nameList = [names[nzList[i]] for i in range(len(nzList))]
    
    print(nameList)
    for i in range(ncols):
        # 绘制各个属性的系数与迭代次数的关系图像
        coefCurve = [betaMat[k][i] for k in range(nSteps)]
        xaxis = range(nSteps)
        plt.plot(xaxis, coefCurve)
    
    plt.xlabel("Steps Taken")
    plt.ylabel(("Coefficient Values"))

print(nameList)的输出结果为：

['alcohol', 'volatile acidity', 'sulphates', 'total sulfur dioxide', 'chlorides', 
    'fixed acidity', 'pH','free sulfur dioxide', 'citric acid', 'residual sugar', 'density']

属性值与迭代次数的关系图像如下图所示：

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQxMDgwODUw_size_16_color_FFFFFF_t_70 1][]

# 通过betaMat的值进一步查看beta参数是如何更新的
    betaMat[20:40]

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQxMDgwODUw_size_16_color_FFFFFF_t_70 2][]

利用交叉验证来确定LARS生成结果中的最佳系数组合：

# 建立交叉验证循环以确定最佳系数组合
    
    # 初始化相关参数
    nxval = 10          # nxval为交叉验证的折数
    nSteps = 350        # 迭代次数
    stepSize = 0.004    #迭代步长
            
    mseMat = []   # 用于存放MSE列表的列表
    for ixval in range(nxval):
        
        m = []    # m用于存放每次迭代后计算的的MSE
        # 划分测试和训练索引集
        idxTest = [a for a in range(nrows) if a%nxval == ixval]
        idxTrain = [a for a in range(nrows) if a%nxval != ixval]
    
        # 划分测试与训练用的属性集和标签集
        xTrain = [xNormalized[r] for r in idxTrain]
        xTest = [xNormalized[r] for r in idxTest]
        labelTrain = [labelNormalized[r] for r in idxTrain]
        labelTest = [labelNormalized[r] for r in idxTest]
    
        # 在训练数据集上训练LARS回归
        nrowsTrain = len(idxTrain)
        nrowsTest = len(idxTest)
    
        # 初始化向量参数beta
        beta = [0.0] * ncols
    
        # 初始化用于存放每次迭代后的beta参数列表的betaMat列表
        betaMat = []
        betaMat.append(list(beta))
        
        # 初始化用于存放误差的列表
        errors = []
        for i in range(nSteps):
            b = []
            errors.append(b)
            
        for iStep in range(nSteps):
            
            # 计算xTrain中的每个样本对应的残差
            residuals = [0.0] * nrowsTrain
            for j in range(nrowsTrain):
                labelsHat = sum([xTrain[j][k] * beta[k] for k in range(ncols)])
                residuals[j] = labelTrain[j] - labelsHat
    
            # 分别计算每个属性与残差的关联值的平均值(两个变量的关联值等于标准化(减去均值后再除以标准差)值的乘积)
            corr = [0.0] * ncols
            for j in range(ncols):
                corr[j] = sum([xTrain[k][j] * residuals[k] for k in range(nrowsTrain)]) / nrowsTrain
    
            # 找出绝对值最大的关联值及其对应的索引
            iStar = 0
            corrStar = corr[0]
            for j in range(1, (ncols)):
                if abs(corrStar) < abs(corr[j]):
                    iStar = j
                    corrStar = corr[j]
    
            # 利用属性与残差的最大关联值(绝对值最大)更新最大关联值对应的beta参数
            beta[iStar] += stepSize * corrStar / abs(corrStar)    # 关联值为正，beta参数对应增加，关联值为负，beta参数对应减小
            betaMat.append(list(beta))      # 将只有一个beta值被更新的beta参数列表存入列表betaMat中
    
            # 基于刚才计算的bate参数在测试数据集上计算预测值，并进一步计算样本外误差
            for j in range(nrowsTest):
                labelsHat = sum([xTest[j][k] * beta[k] for k in range(ncols)])
                err = labelTest[j] - labelsHat
                errors[iStep].append(err)  # error中的元素是列表，列表中的每个元素存储的是对应每次迭代更新后的所有预测值与对应真实值之间的误差
                
        # 计算每次迭代更新beta参数后的MSE
        for errVect in errors:
            mse = sum([x*x for x in errVect])/len(errVect)
            m.append(mse)       # 将每次迭代更新beta参数后计算的MSE存入m列表中
            
        mseMat.append(m)    # 将每次划分训练数据集和测试数据集后的m存入mseMa列表中
    
    
    mseMat = np.array(mseMat)    # 将mseMat转化为array对象
    avg_mse = np.mean(mseMat,axis=0)   # 计算350次迭代产生的350个MSE在10次交叉验证后的均值
    minMse = min(avg_mse)    # 最佳系数组合对应的最小MSE
    minPt = [i for i in range(len(avg_mse)) if avg_mse[i] == minMse ][0]    # 最小MSE对应的最佳系数组合的索引
    print("Minimum Mean Square Error:", minMse)     
    print("Index of Minimum Mean Square Error:", minPt)   
    
    # 绘制MSE与迭代次数的关系图像
    xaxis = range(len(avg_mse))
    plt.plot(xaxis, avg_mse)
    
    plt.xlabel("Steps Taken")
    plt.ylabel(("Mean Square Error"))

上述代码中print("Minimum Mean Square Error", minMse)和print("Index of Minimum Mean Square Error", minPt)的输出结果分别是：

Minimum Mean Square Error: 0.649156243914687和Index of Minimum Mean Square Error： 285

MSE与迭代次数的关系图像如下图所示：

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQxMDgwODUw_size_16_color_FFFFFF_t_70 3][]

由此可以确定最佳系数组合为betaMat\[285\](betaMat是以wine\_data中的所有数据作为训练数据得到的所有系数组合组成的列表)：

betaMat[285]

其值为：  
\[0.008, -0.23200000000000018, -0.012, 0.008, -0.10400000000000006, 0.036000000000000004,

\-0.11200000000000007, 0.0, -0.07600000000000004, 0.17600000000000013, 0.3800000000000003\]

其他参考：

《Python机器学习——预测分析核心算法》Michael Bowles著

[https://en.wikipedia.org/wiki/Correlation\_and\_dependence][https_en.wikipedia.org_wiki_Correlation_and_dependence]

[https://blog.csdn.net/qq\_41080850/article/details/86560645][https_blog.csdn.net_qq_41080850_article_details_86560645]

论文Least Angle Regression及本文用到的红酒数据集的地址：[https://pan.baidu.com/s/1ONfd1pg5B3FH6vkn4hOLlQ][https_pan.baidu.com_s_1ONfd1pg5B3FH6vkn4hOLlQ]

[20190213110426806.PNG]: /images/20220320/2f1d1435952e46f9a20506389df0baa2.png
[20190213110807134.PNG]: /images/20220320/84066de79a8045a096943f61e297a7a9.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQxMDgwODUw_size_16_color_FFFFFF_t_70]: /images/20220320/b0e59ed76f2f49ba80ff2338a91292ac.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQxMDgwODUw_size_16_color_FFFFFF_t_70 1]: /images/20220320/1a197ee37a2b49b9886a98a4108daaff.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQxMDgwODUw_size_16_color_FFFFFF_t_70 2]: /images/20220320/7d3b1837ca3b4b50ad2be912067149a3.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQxMDgwODUw_size_16_color_FFFFFF_t_70 3]: /images/20220320/e4d78875acef434e81569e419cd5f46e.png
[https_en.wikipedia.org_wiki_Correlation_and_dependence]: https://en.wikipedia.org/wiki/Correlation_and_dependence
[https_blog.csdn.net_qq_41080850_article_details_86560645]: https://blog.csdn.net/qq_41080850/article/details/86560645
[https_pan.baidu.com_s_1ONfd1pg5B3FH6vkn4hOLlQ]: https://pan.baidu.com/s/1ONfd1pg5B3FH6vkn4hOLlQ