贝叶斯套索回归

打开脚本

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

Lasso回归是一种结合正则化和变量选择的线性回归技术。正则化通过降低回归系数的大小来防止过拟合。lasso回归的频率主义观点与其他正则化技术(如岭回归)不同，因为lasso将一个值精确地定义为0，对应于不重要或冗余的预测因子。

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

$$ Y = \ beta_0 1_n + X \的β+ \ varepsilon。$$

  • $ Y $是n响应。

  • $ X $是矩阵p相应的观察到预测变量。

  • β\美元是p回归系数。

  • $ \ $ beta_0是截距。

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

  • $ \ $ varepsilon是独立同分布的高斯干扰用0和方差的平均的随机矢量$ \西格马^ 2 $。

拉索回归的频率观的目标函数是

$$˚F\左（{\ beta_0，\β; Y，X} \右）= 0.5 \左\ |{Y - \ beta_0 1_n + X \的β} \右\ | _2 ^ 2 + \拉姆达\左\ |\测试\右\ |。_1 ^ 2 $$

$ \ $拉姆达是惩罚(或收缩)项，与其他参数不同，它不适合数据。您必须指定一个值$ \ $拉姆达在进行评估之前，最佳实践是对其进行调优。在您指定一个值之后，$ F $被最小化相对于该回归系数。最终的系数被套索估计。有关套索回归的频率论观点的详细信息，请参阅[182]。

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

加载和预处理数据

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

负载Data_USEconModel

该数据集包括MATLAB®时间表数据表，包括14个由1947年第一季至2009年第一季所量度的变量;UNRATE是美国的失业率。更多详情，请输入描述在命令行。

以相同的图形绘制所有的系列，但在单独的副情节。

图;为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（@日志，数据表，“数据源”,...DataTable.Properties.VariableNames (hasexpotrend));DataTableLog = [DataTableLog DataTable(:，DataTable. properties . variablenames (~hasexpotrend)];

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

在lasso回归目标函数中，具有较大量值的系数往往占主导地位。因此，重要的是变量有一个类似的规模，当你实施lasso回归。比较变量in的比例DataTableLog通过在同一轴上绘制它们的箱形图。

图;箱线图(DataTableLog.Variables'标签',DataTableLog.Properties.VariableNames);甘氨胆酸h =;h。XTickLabelRotation = 45; title(“变量盒阴谋”);

这些变量的规模相当相似。

为了确定时间序列模型如何概括，步骤如下：

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

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

  3. 使用拟合模型来预测应答到预测范围。

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

  5. 选择具有最低FMSE模型。

创建估计，为响应和预测数据的预测样本变量。指定的4年（16个季度）预测地平线。

FH = 16;Y = DataTableLog.UNRATE（1：（结束 -  FH））;YF = DataTableLog.UNRATE（（端 -  FH + 1）：结束）;isresponse = DataTable.Properties.VariableNames ==“UNRATE”;X = {DataTableLog 1：（结束 -  FH），〜isresponse};XF = DataTableLog {（端 -  FH + 1）：端，〜isresponse};P =尺寸（X，2）;预测的数量％predictornames = DataTableLog.Properties.VariableNames (~ isresponse);

飞度简单线性回归模型

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

SLRMdlFull = fitlm (X, y,“VarNames”，DataTableLog.Properties.VariableNames)， slrPValue = slrmdlfull .系数;isInSig = SLRMdlFull。CoefficientNames (slrPValue > 0.05)
SLRMdlFull =线性回归模型：UNRATE〜[线性公式在13个预测14项]估计系数：估计SE TSTAT p值________ _______ __________（截距）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.38280.16853 log_GCE -5.63 1.1581 -4.8615 2.6245e-06 log_GDP 14.522 3.9406 3.6851 0.00030659 log_GDPDEF 11.926 2.9298 4.0704 7.1598e-05 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.3116e-53 log_M1SL-3.3023 0.53185 -6.209 3.914e-09 log_M2SL -1.3845 1.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观察数：185，错误自由度：171均方根错误：0.312 R平方：0.957，调整R平方：0.954 F统计与常数模型：292，p值= 2.26e-109 isInSig = 1×4单元阵列{ 'log_CPIAUCSL'} {'升og_GS10 '} {' log_M2SL '} {' FEDFUNDS'}

SLRMdlFull是线性模型模型对象。的p-值表明，在5%的显著性水平下，CPIAUCSL,GS10,M2SL,FEDFUNDS是无关紧要的变量。

重新装上线性模型没有预测数据的显着变量。使用完全和简化模型，预测失业率到预测范围，则估计FMSEs。

idxreduced =〜ismember（predictornames，isInSig）;XReduced = X（：，idxreduced）;SLRMdlReduced = fitlm（XReduced，Y，“VarNames”[predictornames (idxreduced) {'UNRATE'}]）;yFSLRF =预测（SLRMdlFull，XF）;yFSLRR =预测（SLRMdlReduced，XF（：，idxreduced））;fmseSLRF = SQRT（均值（（YF  - 。yFSLRF）^ 2））fmseSLRR = SQRT（均值（（YF  - 。yFSLRR）^ 2））
fmseSLRF = 0.6674 fmseSLRR = 0.6105

简化模型的FMSE小于全模型的FMSE，说明简化模型的泛化效果更好。

适合频率论套索回归模型

根据数据拟合lasso回归模型。使用默认的正则化路径。套索默认情况下，标准化的数据。由于变量是类似的尺度，不规范他们。返回有关沿着正规化路径的千篇一律的信息。

[LassoBetaEstimates，FitInfo] =套索（X，Y，“标准化”假的,...“PredictorNames”,predictornames);

默认情况下，套索通过使用一个值序列对lasso回归模型进行100次拟合$ \ $拉姆达。因此,LassoBetaEstimates是一个13×100的回归系数估计矩阵。行对应于in中的预测变量X和列对应于值$ \ $拉姆达。FitInfo是包括的值的结构$ \ $拉姆达(FitInfo.Lambda）时，均方误差为各装配（FitInfo.MSE），以及所估计的截距（FitInfo.Intercept)。

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

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

绘制回归系数的幅度相对于所述收缩值。

HAX = lassoPlot（LassoBetaEstimates，FitInfo）;L1Vals = hax.Children.XData;yyaxis正确的h =情节(L1Vals fmseLasso,'行宽'，2，'的LineStyle',“——”);图例（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表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 FEDFUNDS 0 TB3MS -0.14459

该频率论套索分析表明，变量CPIAUCSL,GCE,GDP,GPDI,PCEC,FEDFUNDS不是无关紧要就是多余。

拟合贝叶斯Lasso回归模型

在套索回归的贝叶斯视图，回归系数的先验分布是拉普拉斯（双指数），均值为0，比例σ\ / \λ美元，其中$ \ $拉姆达固定的收缩参数和$ \西格马^ 2 \ SIM IG（A，B）$。有关详细信息，请参阅lassoblm。

与拉索回归的频率观一样，如果你增加$ \ $拉姆达，然后将得到的模型的稀疏性单调增加。然而，不同于频率论套索，贝叶斯套索具有不完全为0。相反，估计有一个后验分布和不显着的或冗余的，如果周围0区域具有高密度微不足道的或冗余的系数模式。

在MATLAB®中实现Bayesian lasso回归时，请注意统计和机器学习工具箱功能之间的几个差异套索和计量经济学工具箱™对象lassoblm和其相关的功能。

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

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

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

由于lassoblm框架是适合执行贝叶斯拉索回归收缩值每系数，你可以使用一个for循环执行拉索正则化路径。

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

PriorMdl = bayeslm (p,“ModelType”,“套索”,“VarNames”,predictornames);表(PriorMdl.Lambda,'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 =任何（isnan（X），2）;N =总和（〜ISMISSING）;有效样本量%λ= FitInfo.Lambda * n /√(FitInfo.MSE);

考虑正则化路径套索。循环通过每次迭代的收缩值:

  1. 估计回归系数和干扰方差的给定的收缩率和数据的后验分布。因为预测值的尺度接近，属性相同的收缩每个预测，和属性的收缩0.01到拦截。保存后的手段和95％置信区间。

  2. 套索情节，如果95％的可信区间包括0，则对应的后验均值设置为零。

  3. 预测进入预测范围。

  4. 估计FMSE。

numlambda =元素个数(λ);% 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。λ=λ(j);[EstMdl,总结]=估计(PriorMdl, X, y,'显示'、假);BayesLassoCoefficients (:, j) =总结。意思是(1:结束(- 1));BLCPlot (:, j) =总结。意思是(1:结束(- 1));BayesLassoCI95 (:,:, j) =总结。CI95 (1: (- 1),:);idx = BayesLassoCI95(:，2,j) > 0 & BayesLassoCI95(:，1,j) <= 0;BLCPlot (idx j) = 0;yFBayesLasso =预测(EstMdl XF);fmseBayesLasso(j) = sqrt(mean(((yF - yFBayesLasso).^2));结束

在每次迭代中，估计运行吉布斯采样从全样本条件表达式（见易于分析后验）并返回empiricalblm模型EstMdl，其中包含绘图和评估汇总表摘要。您可以为Gibbs sampler指定参数，例如绘制的数量或细化方案。

一个好的做法是通过检查跟踪地块诊断后的样品。然而，由于100个分布绘制，本实施例情况下，不执行该步骤。

积系数的后部手段相对于归一化的L1处罚（系数的除截距幅度的总和）。在同样的情节，但在右边的Y轴，绘制FMSEs。

L1Vals =总和（绝对值（BLCPlot（2：端，:)），1）/最高（总和（ABS（BLCPlot（2：端，:)），1））;图;情节（L1Vals，BLCPlot（2：端，:)）xlabel（“L1”);ylabel（“系数估计”);yyaxis正确的H =情节（L1Vals，fmseBayesLasso，'行宽'，2，'的LineStyle',“——”);图例（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 -5.1707 log_GCE -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 0.22748 1.4237 log_HOANBS  -32.538 -35.589 -29.526 log_M1SL -3.0656 -4.1499 -2.0066 -1.1007 log_M2SL -3.1243 0.75844 log_PCEC -3.3342 -9.6496 1.8482 FEDFUNDS 0.043672 -0.056192 0.14254 TB3MS -0.15682 -0.29385 -0.022145

以确定是否变量是微不足道的或冗余的一种方式是通过检查是否0处于其相应系数95％的可信区间。使用该标准，CPIAUCSL,M2SL,PCEC,FEDFUNDS是微不足道的。

在表中显示所有的估计进行比较。

SLRR = 0 (p + 1,1);回归系数;回归系数表([SLRMdlFull.Coefficients.Estimate;fmseSLRR),...[SLRR;fmseSLRR),...[bestLasso;fmsebestlasso),...[bestBayesLasso;fmsebestbayeslasso),'VariableNames',...{'SLRFull''SLRReduced''套索'“BayesLasso”}，...'RowNames',(PriorMdl.VarNames;'MSE']）
ANS = 15x4表SLRFull SLRReduced套索BayesLasso ________ __________ ________ __________截取88.014 88.327 61.154 90.587 log_COE 7.1111 10.854 0.042061 6.5591 log_CPIAUCSL -3.4032 0 0 -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.538 log_M1SL -3.3023 -2.8099 -1.3365 -3.0656 -1.3845 log_M2SL 0 0.17948 -1.1007 log_PCEC -7.143 -14.207 0 -3.3342 FEDFUNDS 0.044542 0 0 0.043672 TB3MS  -0.1577 -0.078449 -0.14459 -0.15682 MSE 0.61048 0.61048 0.79425 0.69639

选择一个模型后，使用整个数据集重新估计它。例如，创建一个预测贝叶斯lasso回归模型，创建一个先验模型，并指定产生最小FMSE的最简单模型的收缩，然后使用整个数据集估计它。

bestLambda =λ(idx);PriorMdl = bayeslm (p,“ModelType”,“套索”,“λ”bestLambda,...“VarNames”,predictornames);PosteriorMdl =估计(PriorMdl [X;XF]、[y;yF]);
方法:lasso 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经验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经验FEDFUNDS | -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实证

另请参阅

||

相关话题

计量经济学工具文档
金宝app
实用指南建模金融风险与MATLAB

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

下载电子书