主要内容

贝叶斯套索回归

打开脚本

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

拉索回归是一种结合了正则化和变量选择的线性回归技术。正则化通过降低回归系数的大小来防止过拟合。套索回归的频率主义者观点与其他正则化技术(如岭回归)不同，因为套索将一个完全为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 $相对于回归系数是最小的。得到的系数是套索估计。有关套索回归的频率论观点的更多细节，请参见[186]．

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

加载和预处理数据

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

负载Data_USEconModel

数据集包括MATLAB®时间表数据表，其中包含从1947年Q1到2009年Q1的14个变量;UNRATE是美国的失业率。有关详细信息，请输入描述在命令行。

在同一个图中绘制所有系列，但在单独的副图中。

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

所有系列除FEDFUNDS，GS10，TB3MS,UNRATE呈现指数趋势。

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

hasexpotrend = ~ismember(DataTable.Properties.VariableNames，.．.［“FEDFUNDS”“GD10”“TB3MS”“UNRATE”]);DataTableLog = varfun(@log,DataTable，“数据源”，.．.DataTable.Properties.VariableNames (hasexpotrend));DataTableLog = [DataTableLog DataTable(:，DataTable. properties . variablenames (~hasexpotrend))];

DataTableLog是一个时间表数据表，但那些具有指数趋势的变量在对数尺度上。

在套索回归目标函数中，具有较大量级的系数往往占惩罚的主导地位。因此，当您实现套索回归时，变量具有相似的尺度是很重要的。比较中变量的比例DataTableLog通过在同一轴上画箱形图。

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

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

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

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

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

  3. 使用拟合模型对预测视界内的响应进行预测。

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

  5. 选择FMSE最低的模型。

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

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

拟合简单线性回归模型

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

SLRMdlFull = fitlm(X,y，“VarNames”slrPValue = slrmdlfull . coefficient . pvalue;isinsing = SLRMdlFull。系数名称(slrPValue > 0.05)
SLRMdlFull =线性回归模型:UNRATE ~[13个预测因子中14项的线性公式]估计SE tStat pValue ________ ________ _______ __________(拦截)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_GDP 14.522 3.9406 - 4.8651 0.00030659 log_gddpdef 11.926 2.9298 4.0704 7.1598e-05 log_GPDI -2.5939 0.54603 -4.7504 4.2792e-06 log_GS10 0.7467 0.29313 1.9605 0.051565 log_HOANBS -32.861 1.4507 -22.652 2.316e -53 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, issising);XReduced = X(:，idxreduced);SLRMdlReduced = fitlm(XReduced,y，“VarNames”[predictornames (idxreduced) {“UNRATE”}));yFSLRF = predict(SLRMdlFull,XF);yFSLRR = predict(SLRMdlReduced,XF(:，idxreduced));fmseSLRF = sqrt(mean((yF - yFSLRF).^2))
fmseSLRF = 0.6674 fmseSLRR = 0.6105

简化模型的FMSE小于完整模型的FMSE，表明简化模型具有更好的泛化能力。

拟合频率套索回归模型

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

[LassoBetaEstimates,FitInfo] =套索(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. fmseLasso)。Df == 7));idx = fmseLasso == fmsebestlasso;bestLasso = [FitInfo.Intercept(idx);LassoBetaEstimates (:, idx)];表(bestLasso,“RowNames”，[“拦截”predictornames])
ans = 14x1表bestLasso _________拦截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联邦基金0 TB3MS -0.14459

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

拟合贝叶斯套索回归模型

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

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

当您在MATLAB®中实现贝叶斯套索回归时，请注意统计和机器学习工具箱™函数之间的几个区别套索和计量经济学工具箱™对象lassoblm和它的相关函数。

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

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

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

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

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

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

PriorMdl是一个lassoblm模型对象。lassoblm的收缩1对于除截距外的其他系数，截距的收缩量为0.01．这些默认值在大多数应用程序中可能没有用;最佳实践是对许多值进行试验。

考虑返回的正则化路径套索．将收缩值膨胀一个因子$ n / \√(MSE_j)美元,在那里MSE_j美元套索运行的均方误差是多少$ j $，$ j $通过收缩值的数量= 1。

ismissing = any(isnan(X)，2);N = sum(~ismissing);有效样本量%lambda = FitInfo.Lambda*n./sqrt(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] = estimate(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 =预测(EstMdl,XF);fmseBayesLasso(j) =√(均值((yF - yFBayesLasso).^2));结束

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

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

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

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表bestBayesLasso bestCI95 ______________ ______________________截获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_GPDI -2.0973 -3.2168 -0.96728 log_GDPDEF 10.281 5.1677 15.936 log_GPDI -2.0973 -3.2168 -0.96728 log_GS10 0.82485 0.22748 1.4237 log_HOANBS -32.538 - 3.656 -4.1499 -2.0066 log_M2SL -1.1007 -3.1243 0.75844 log_PCEC -3.3342 -9.6496 1.8482联邦基金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表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_GDPDEF 14.522 16.701 0 9.8839 log_gddpdef 11.926 9.1106 8.524 10.281 log_GPDI -2.5939 -2.6963 0 -2.0973 log_GS10 0.7467 0 1.6957 0.82485 log_HOANBS -32.861 -33.782 -18.937 -32.538 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 = estimate(PriorMdl，[X;XF]、[y;yF]);
方法:套索MCMC采样10000张图观察数:201预测因子数:14 |意味着性病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经验日志_gpdi | -2.5884 0.5294[-3.628， -1.553] 0.000经验日志_gs10 | 0.6910 0.2577[0.194, 1.197] 0.997经验日志_hoanbs | -32.3356 1.4941[-35.276， -29.433] 0.000经验日志_m1sl | -2.2279 0.5043[-3.202， -1.238] 0.000经验日志_m2sl | 0.3617 0.9123[-1.438， - 1.179] 0.652经验日志_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]经验Sigma2 | 0.1186 0.0122[0.097, 0.145] 1.000经验的

另请参阅

||

相关的话题

计量经济学工具箱文档

万博1manbetx

用MATLAB建模金融风险的实用指南

用MATLAB建模金融风险的实用指南

下载电子书