主要内容

贝叶斯套索回归

打开脚本

这个例子展示了如何使用贝叶斯套索回归来执行变量选择。

套索回归是一种结合了正则化和变量选择的线性回归技术。正则化通过降低回归系数的大小来帮助防止过拟合。套索回归的频率论观点不同于其他正则化技术，如岭回归，因为套索将一个恰好为0的值归为对应于无关紧要或冗余的预测因子的回归系数。

考虑这个多元线性回归模型:

$$y = \beta_0 1_n + X\beta + \varepsilon.$$

  • y美元的向量n响应。

  • X美元是一个矩阵p对应观测到的预测变量。

  • β\美元的向量p回归系数。

  • \ beta_0美元是截距。

  • 1) _n”美元是一个长度n1的列向量。

  • \ varepsilon美元是iid高斯扰动的随机向量，其均值为0，方差为\σ^ 2美元．

套索回归频率论观点的目标函数为

$ $ f \离开({β\ \ beta_0; y、X} \右)= 0.5 \ \ | {y -β\ beta_0 1 _n”+ X \} \右\ | _2 ^ 2 +λ\ \左右\ | \β\ \ | _1 ^ 2。$ $

\λ美元是惩罚(或收缩)项，与其他参数不同，它不适合数据。必须为指定一个值\λ美元在评估之前，最好的实践是调优它。指定一个值之后，$ f $相对于回归系数是最小的。得到的系数就是套索估计。有关套索回归的频率论视图的更多细节，请参见[193]．

对于这个例子，考虑为美国失业率创建一个预测线性模型。你需要一个能很好概括的模型。换句话说，您希望通过删除所有冗余预测因子和所有与失业率不相关的预测因子来最小化模型的复杂性。

加载和预处理数据

加载美国宏观经济数据集Data_USEconModel.mat．

负载Data_USEconModel

数据集包括MATLAB®时间表DataTimeTable，其中包含从1947年第一季度到2009年第一季度测量的14个变量;UNRATE是美国的失业率。请输入描述在命令行。

将所有系列画在同一个图中，但在不同的子图中。

图;为j = 1:size(DataTimeTable,2) subplot(4,4,j) plot(DataTimeTable. time,DataTimeTable{:，j});标题(DataTimeTable.Properties.VariableNames (j));结束

所有系列除外FEDFUNDS，GS10，TB3MS,UNRATE似乎呈指数增长趋势。

对具有指数趋势的变量应用对数变换。

hasexpotrend = ~ismember(datatitable . properties . variablenames，.．.［“FEDFUNDS”“GD10”“TB3MS”“UNRATE”]);DataTimeTableLog = varfun(@日志，DataTimeTable，“数据源”，.．.DataTimeTable.Properties.VariableNames (hasexpotrend));DataTimeTableLog = [DataTimeTableLog DataTimeTable(:， datatitable . properties . variablenames (~hasexpotrend))];

DataTimeTableLog时间表就像DataTimeTable，但那些具有指数趋势的变量是对数尺度。

在套索回归目标函数中，具有较大幅度的系数往往占惩罚的主导地位。因此，在实现套索回归时，变量具有相似的规模是很重要的。比较变量的尺度DataTimeTableLog通过在同一轴上绘制它们的箱形图。

图;箱线图(DataTimeTableLog。变量,“标签”, DataTimeTableLog.Properties.VariableNames);H = gca;h.XTickLabelRotation = 45;标题(“可变箱形图”）;

这些变量有相当相似的尺度。

要确定时间序列模型的泛化效果如何，请遵循以下步骤:

  1. 将数据划分为估计样本和预测样本。

  2. 将模型拟合到估计样本。

  3. 使用拟合模型来预测预测范围内的响应。

  4. 估计每个模型的预测均方误差(FMSE)。

  5. 选择FMSE最低的模型。

为响应和预测数据创建估计和预测样本变量。指定4年(16个季度)的预测期限。

Fh = 16;y = DataTimeTableLog。UNRATE（1:(end - fh)); yF = DataTimeTableLog.UNRATE((end - fh + 1):end); isresponse = DataTimeTable.Properties.VariableNames ==“UNRATE”；X = DataTimeTableLog{1:(end - fh)，~isresponse};XF = DataTimeTableLog{(end - fh + 1):end，~isresponse};p = size(X,2);%预测因子的数量predictornames = DataTimeTableLog.Properties.VariableNames(~isresponse);

拟合简单线性回归模型

对估计样本拟合一个简单的线性回归模型。指定变量名(响应变量的名称必须是最后一个元素)。提取p值。通过使用5%的显著性水平来识别不显著系数。

SLRMdlFull = fitlm(X,y，“VarNames”，DataTimeTableLog.Properties.VariableNames) slrPValue = SLRMdlFull.Coefficients.pValue;isessig = SLRMdlFull。系数名称(slrPValue > 0.05)
SLRMdlFull =线性回归模型:UNRATE ~[13个预测因子中14项的线性公式]估计SE tStat pValue ________ ________ _______ __________ (Intercept) 88.014 5.3229 16.535 2.6568e-37 log_COE 7.1111 2.3426 3.0356 0.0027764 log_CPIAUCSL -3.4032 2.4611 -1.3828 0.16853 log_GCE -5.63 1.1581 -4.8615 2.6245e-06 log_GDPDEF 11.522 3.9406 3.6851 0.00030659 log_GPDI -2.5939 0.54603 -4.7504 4.2792e-06 log_GS10 0.57467 0.29313 1.9605 0.051565 log_HOANBS -32.861 1.4507 -22.652 2.316e -53 log_M1SL -3.3023 0.53185 -6.209 3.914e-09 log_M2SL -1.38451．0438 -1.3264 0.18647 log_PCEC -7.143 3.4302 -2.0824 0.038799 FEDFUNDS 0.044542 0.047747 0.93287 0.3522 TB3MS -0.1577 0.064421 -2.448 0.015375 Number of observations: 185, Error degrees of freedom: 171 Root Mean Squared Error: 0.312 R-squared: 0.957, Adjusted R-Squared: 0.954 F-statistic vs. constant model: 292, p-value = 2.26e-109 isInSig = 1x4 cell array {'log_CPIAUCSL'} {'log_GS10'} {'log_M2SL'} {'FEDFUNDS'}

SLRMdlFull是一个LinearModel模型对象。的p-值表明，在5%显著性水平下，CPIAUCSL，GS10，M2SL,FEDFUNDS都是不重要的变量。

对预测数据中不存在不显著变量的线性模型进行改装。利用完整模型和简化模型，将失业率预测到预测范围内，然后估计FMSEs。

idxreduced = ~ismember(predictornames, isessig);XReduced = X(:，idxreduced);SLRMdlReduced = fitlm(XReduced,y，“VarNames”[predictornames (idxreduced) {“UNRATE”}));yFSLRF = predict(SLRMdlFull,XF);yFSLRR = predict(SLRMdlReduced,XF(:，idxreduced));fmseSLRR =√(mean((yF - yFSLRF).^2))
fmseSLRF = 0.6674 fmseSLRR = 0.6105

简化模型的FMSE小于完整模型的FMSE，说明简化模型具有更好的泛化性。

拟合频率套索回归模型

用套索回归模型拟合数据。使用默认的正则化路径。套索默认情况下标准化数据。因为变量的规模相似，所以不要将它们标准化。返回关于正则化路径上拟合的信息。

[LassoBetaEstimates,FitInfo] = lasso(X,y，“标准化”假的,.．.“PredictorNames”, predictornames);

默认情况下,套索的值序列拟合套索回归模型100次\λ美元．因此,LassoBetaEstimates是一个13 × 100的回归系数估计值矩阵。中的行对应预测变量X，列对应的值\λ美元．FitInfo结构是否包含的值\λ美元（FitInfo。λ)，每次拟合的均方误差(FitInfo。均方误差)，以及估计截距(FitInfo。拦截）.

计算返回的每个模型的FMSE套索．

yFLasso = FitInfo。拦截+ XF*LassoBetaEstimates; fmseLasso = sqrt(mean((yF - yFLasso).^2,1));

画出回归系数与收缩值的关系。

hax = lassoPlot(LassoBetaEstimates,FitInfo);L1Vals = hax.Children.XData;yyaxis正确的h = plot(L1Vals,fmseLasso，“线宽”2,“线型”，“——”）;传奇(h,“FMSE”，“位置”，“西南”）;ylabel (“FMSE”）;标题(“拉索频率论者”）

具有7个预测因子(df = 7)的模型似乎很好地平衡了最小FMSE和模型复杂性。得到包含7个预测因子且FMSE最小的模型对应的系数。

fmsebestlasso = min(fmseLasso(FitInfo。Df == 7));idx = fmseLasso == fmsebestlasso;bestLasso = [FitInfo.Intercept(idx);LassoBetaEstimates (:, idx)];表(bestLasso,“RowNames”, (“拦截”predictornames])
ans = 14x1 table bestLasso _________ Intercept 61.154 log_COE 0.042061 log_CPIAUCSL 0 log_GCE 0 log_GDP 0 log_GDPDEF 8.524 log_GPDI 0 log_GS10 1.6957 log_HOANBS -18.937 log_M1SL -1.3365 log_M2SL 0.17948 log_PCEC 0 FEDFUNDS 0 TB3MS -0.14459

频率论套索分析表明，变量CPIAUCSL，全球教育运动，国内生产总值，GPDI，PCEC,FEDFUNDS不是无关紧要就是多余的。

拟合贝叶斯拉索回归模型

在套索回归的贝叶斯观点中，回归系数的先验分布是拉普拉斯(双指数)，均值为0，有尺度σ\ / \λ美元,在那里\λ美元固定收缩参数和美元\σ^ 2 \ sim搞笑(A, B)美元．详情请参见lassoblm．

就像套索回归的频率论观点一样，如果你增加\λ美元，则所得模型的稀疏性单调递增。然而，与频率套索不同，贝叶斯套索具有不重要或冗余的系数模式，这些模式并不完全为0。相反，估计具有后验分布，如果0周围的区域具有高密度，则估计是不重要或冗余的。

在MATLAB®中实现贝叶斯套索回归时，请注意统计数据和机器学习工具箱™函数之间的几个差异套索和Econometrics Toolbox™对象lassoblm以及相关的函数。

  • lassoblm是对象框架的一部分，而套索是一个函数。目标框架简化了计量经济学工作流程。

  • 不像套索，lassoblm没有标准化预测器数据。但是，您可以为每个系数提供不同的收缩值。这一特性有助于保持系数估计值的可解释性。

  • lassoblm对每个系数应用一个收缩值;它不接受这样的正则化路径套索．

因为lassoblm框架适用于对每个系数的收缩值执行贝叶斯套索回归，您可以使用for循环在正则化路径上执行套索。

建立贝叶斯套索回归先验模型bayeslm．指定预测变量的数量和变量名称。控件中存储的每个系数的默认收缩值λ模型的属性。

PriorMdl = bayeslm(p，“ModelType”，“套索”，“VarNames”, predictornames);表(PriorMdl。λ,“RowNames”PriorMdl.VarNames)
ans = 14x1表Var1 ____拦截0.01 log_COE 1 log_CPIAUCSL 1 log_GCE 1 log_GDP 1 log_GDPDEF 1 log_GPDI 1 log_GS10 1 log_HOANBS 1 log_M1SL 1 log_M2SL 1 log_PCEC 1 FEDFUNDS 1 TB3MS 1

PriorMdl是一个lassoblm模型对象。lassoblm收缩的属性1对于每个系数，除了截距，它有收缩0.01．这些默认值在大多数应用程序中可能没有用处;最好的做法是尝试多种价值。

考虑返回的正则化路径套索．将收缩率膨胀到×的倍数$ n / \√(MSE_j)美元,在那里MSE_j美元套索的MSE是否运行$ j $，$ j $= 1到收缩值的数目。

ismissing = any(isnan(X)，2);N = sum(~ismissing);有效样本量%lambda = FitInfo.Lambda*n./√(FitInfo.MSE);

考虑返回的正则化路径套索．循环遍历每次迭代的收缩值:

  1. 估计回归系数和扰动方差的后验分布，给定收缩和数据。由于预测因子的尺度接近，将相同的收缩归因于每个预测因子，并将收缩归因于0.01去拦截。存储后验均值和95%可信区间。

  2. 对于套索图，如果95%可信区间包含0，则将相应的后验均值设置为0。

  3. 预测进入预测范围。

  4. 估计FMSE。

Numlambda = numel(lambda);% PreallocateBayesLassoCoefficients = 0 (p+1,numlambda);BayesLassoCI95 = 0 (p+1,2,numlambda);fmseBayesLasso = 0 (numlambda,1);BLCPlot = 0 (p+1,numlambda);%估计和预测rng (10);%用于再现性为j = 1:numlambda PriorMdl。Lambda = Lambda (j);[EstMdl,Summary] =估计(PriorMdl,X,y，“显示”、假);BayesLassoCoefficients(:，j) =摘要。均值(1:(end - 1));BLCPlot(:，j) =摘要。均值(1:(end - 1));BayesLassoCI95(:，:，j) =摘要。CI95(1:(end - 1)，:);idx = BayesLassoCI95(:，2,j) > 0 & BayesLassoCI95(:，1,j) <= 0;BLCPlot(idx,j) = 0;yFBayesLasso = forecast(EstMdl,XF);fmseBayesLasso(j) =√(mean((yF - yFBayesLasso).^2));结束

在每次迭代中，估计运行吉布斯采样器对完整条件进行采样(参见易于分析处理的后验)，并返回empiricalblm模型EstMdl，其中包含图纸和估算汇总表总结．您可以为吉布斯采样器指定参数，例如拉伸次数或细化方案。

一个很好的做法是通过检查痕迹图来诊断后验样本。但是，因为绘制了100个分布，所以本例不执行该步骤。

绘制关于归一化L1惩罚的系数的后验均值(除截距外的系数的大小之和)。在同一张图上，但在右边y轴上，画出fmse。

L1Vals =总和(abs (BLCPlot(2:最终,:)),1)/ max(总和(abs (BLCPlot(2:最终,:)),1));图;情节(L1Vals BLCPlot(2:最终,:))包含(“L1”）;ylabel (的系数估计）;yyaxis正确的h = plot(L1Vals,fmseBayesLasso，“线宽”2,“线型”，“——”）;传奇(h,“FMSE”，“位置”，“西南”）;ylabel (“FMSE”）;标题(“贝叶斯套索”）

当标准化L1惩罚增加到超过0.1时，模型倾向于更好地泛化。为了平衡最小的FMSE和模型复杂性，选择FMSE最接近0.68的最简单的模型。

idx = find(fmseBayesLasso > 0.68,1);fmsebestbayeslasso = fmseBayesLasso(idx);bestBayesLasso = BayesLassoCoefficients(:，idx);bestCI95 = BayesLassoCI95(:，:，idx);表(bestBayesLasso bestCI95,“RowNames”, (“拦截”predictornames])
ans = 14x2 table bestBayesLasso bestCI95 ______________ ______________________ Intercept 90.587 79.819 101.62 log_COE 6.5591 1.8093 11.239 log_CPIAUCSL -2.2971 -6.9019 1.7057 log_GCE -5.1707 -7.4902 -2.8583 log_GDP 9.8839 2.3848 17.672 log_GDPDEF 10.281 5.1677 15.936 log_GPDI -2.0973 -3.2168 -0.96728 log_GS10 0.82485 - 3.22748 1.4237 log_HOANBS -32.538 -4.1499 - 2.0026 log_M1SL -3.0656 - 4.1299 - 0.4844 log_PCEC -3.3342 -9.6496 1.8482 FEDFUNDS 0.043672 -0.0561920.14254 tb3ms -0.15682 -0.29385 -0.022145

确定变量是否不重要或冗余的一种方法是检查0是否在其对应的系数95%可信区间内。使用这个标准，CPIAUCSL，M2SL，PCEC,FEDFUNDS是无关紧要的。

在表格中显示所有估算值，以便进行比较。

SLRR = 0 (p + 1,1);SLRR([true idxreduced]) = SLRMdlReduced.Coefficients.Estimate;表([SLRMdlFull.Coefficients.Estimate;fmseSLRR),.．.[SLRR;fmseSLRR),.．.[bestLasso;fmsebestlasso),.．.[bestBayesLasso;fmsebestbayeslasso),“VariableNames”，.．.｛“SLRFull”“SLRReduced”“套索”“BayesLasso”}，.．.“RowNames”, (PriorMdl.VarNames;MSE的])
ans = 15x4 table SLRFull SLRReduced Lasso BayesLasso ________ __________ ________ __________ Intercept 88.014 88.327 61.154 90.587 log_COE 7.1111 10.854 0.042061 6.5591 log_CPIAUCSL -3.4032 00 -2.2971 log_GCE -5.63 -6.1958 0 -5.1707 log_GDP 14.522 16.701 0 9.8839 log_GDPDEF 11.926 9.1106 8.524 10.281 log_GPDI -2.5939 -2.6963 0 -2.0973 log_GS10 0.57467 0 1.6957 0.82485 log_HOANBS -32.861 -33.782 -18.937 - 32.585 log_M1SL -3.3023 -2.8099 -1.3365 -3.0656 log_M2SL -1.3845 0 0.17948 -1.1007 log_PCEC-7.143 -14.207 0 -3.3342联邦基金0.044542 0 0 0.043672 tb3ms -0.1577 -0.078449 -0.14459 -0.15682 mse 0.61048 0.61048 0.79425 0.69639

在你选择了一个模型之后，用整个数据集重新估计它。例如，要创建一个预测贝叶斯套索回归模型，请创建一个先验模型，并指定收缩量，生成具有最小FMSE的最简单模型，然后使用整个数据集对其进行估计。

bestLambda = lambda(idx);PriorMdl = bayeslm(p，“ModelType”，“套索”，“λ”bestLambda,.．.“VarNames”, predictornames);PosteriorMdl =估计(PriorMdl，[X;XF]、[y;yF]);
方法:套索MCMC抽样10000张，观测数:20114 |意味着性病CI95积极的分布  ----------------------------------------------------------------------------- 拦截| 85.9643 - 5.2710[75.544,96.407]1.000经验log_COE | 12.3574 - 2.1817[8.001, 16.651] 1.000经验log_CPIAUCSL | 0.1498 - 1.9150[-3.733, 4.005] 0.525经验log_GCE | -7.1850 - 1.0556[-9.208, -5.101] 0.000经验log_GDP | 11.8648 - 3.9955[4.111, 19.640] 0.999经验log_GDPDEF | 8.8777 - 2.4221 (4.033,13.745] 1.000经验log_GPDI | -2.5884 0.5294[-3.628， -1.553] 0.000经验log_GS10 | 0.6910 0.2577[0.194, 1.197] 0.997经验log_HOANBS | -32.3356 1.4941[-35.276， -29.433] 0.000经验log_M1SL | -2.2279 0.5043[-3.202， -1.238] 0.000经验log_M2SL | 0.3617 0.9123[-1.438, 2.179] 0.652经验log_PCEC | -11.3018 3.0353[-17.318， -5.252] 0.000经验联邦基金| -0.0132 0.0490[-0.109,0.082]0.392经验TB3MS | -0.1128 0.0660 [-0.244，0.016] 0.042经验Sigma2 | 0.1186 0.0122[0.097, 0.145] 1.000经验

另请参阅

||

相关的话题