贝叶斯套索回归
这个例子展示了如何通过使用贝叶斯套索回归执行可变选择。
Lasso回归是一种结合正则化和变量选择的线性回归技术。正则化通过降低回归系数的大小来防止过拟合。lasso回归的频率主义观点与其他正则化技术(如岭回归)不同，因为lasso将一个值精确地定义为0，对应于不重要或冗余的预测因子。
考虑这个多元线性回归模型：
是的向量n响应。
是矩阵p相应的观察到预测变量。
是的向量p回归系数。
是截距。
是一个长度n者的列向量。
是独立同分布的高斯干扰用0和方差的平均的随机矢量
。
拉索回归的频率观的目标函数是
是点球（或收缩）和术语，不同于其他参数，不适合该数据。您必须指定一个值在进行评估之前，最佳实践是对其进行调优。在您指定一个值之后，
被最小化相对于该回归系数。最终的系数被套索估计。有关套索回归的频率论观点的详细信息，请参阅[182]。
在这个例子中，考虑为美国的失业率预测线性模型。你想一般化很好的典范。换句话说，你要通过删除所有多余的预测，并且是不相关的失业率所有预测，以尽量减少模型的复杂性。
加载和数据预处理
加载美国宏观经济数据集Data_USEconModel.mat。
加载Data_USEconModel
该数据集包括MATLAB®时间表数据表，其中包含从Q1 1947至2009年Q1测量14个变量;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.Properties.VariableNames（〜hasexpotrend））];
DataTableLog是一个时间表数据表，但是那些具有指数趋势的变量是在对数尺度上的。
具有相对大的幅度系数往往主宰在套索回归目标函数的惩罚。因此，重要的是变量有类似的规模，当你实现套索回归。比较变量的尺度在DataTableLog通过在同一轴线上绘制自己的箱形图。
图;箱线图(DataTableLog.Variables'标签'，DataTableLog.Properties.VariableNames）;H = GCA;h.XTickLabelRotation = 45;标题（“变量箱线图”);
变量有相当类似的规模。
为了确定时间序列模型如何概括，步骤如下：
将数据划分为估计和预测的样品。
适合车型估计样品。
使用拟合模型来预测应答到预测范围。
估计每个模型的预测均方误差（FMSE）。
选择具有最低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);
默认，套索通过使用值的一个序列符合套索回归模型100倍。因此,LassoBetaEstimates是一个13×100的回归系数估计矩阵。行对应于in中的预测变量X和列对应于值。FitInfo是包括的值的结构(FitInfo.Lambda）时，均方误差为各装配（FitInfo.MSE），以及所估计的截距（FitInfo.Intercept）。
通过计算返回的每个模型的FMSE套索。
yFLasso = FitInfo.Intercept + XF * LassoBetaEstimates;fmseLasso = SQRT（均值（（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 =分钟（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，比例
，其中是固定的收缩参数和
。有关详细信息，请参阅lassoblm。
与拉索回归的频率观一样，如果你增加，然后将得到的模型的稀疏性单调增加。然而，不同于频率论套索，贝叶斯套索具有不完全为0。相反，估计有一个后验分布和不显着的或冗余的，如果周围0区域具有高密度微不足道的或冗余的系数模式。
在MATLAB®中实现Bayesian lasso回归时，请注意统计和机器学习工具箱功能之间的几个差异套索和计量经济学工具箱™对象lassoblm和其相关的功能。
lassoblm是一个对象的框架的一部分，而套索是一个函数。对象框架简化了工作流程，计量经济学。不像
套索,lassoblm没有标准化的预测数据。但是，你可以为每个系数提供不同的收缩值。此功能可帮助维持系数估计的可解释性。lassoblm将一个收缩值以每个系数;它不接受像正规化路径套索。
由于lassoblm框架适合于每一个系数收缩值执行贝叶斯套索回归，可以使用一个for循环在正规化路径执行套索。
通过创建一个贝叶斯套索回归模型前bayeslm。指定预测变量和变量名的数量。显示默认收缩值存储在每个系数LAMBDA模型的属性。
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。这些默认值可能不是在大多数应用中;最好的做法是有许多值进行试验。
考虑返回正规化路径套索。的因数膨胀收缩值
，其中
是套索运行的MSE
,= 1通过收缩值的个数。
ISMISSING =任何（isnan（X），2）;N =总和（〜ISMISSING）;有效样本量%λ= FitInfo.Lambda * n /√(FitInfo.MSE);
考虑返回正规化路径套索。通过在每次迭代时的收缩值循环：
估计回归系数和干扰方差的给定的收缩率和数据的后验分布。因为预测值的尺度接近，属性相同的收缩每个预测，和属性的收缩
0.01到拦截。保存后的手段和95％置信区间。套索情节，如果95％的可信区间包括0，则对应的后验均值设置为零。
预测到预测范围。
估计FMSE。
numlambda =元素个数(λ);％预分配BayesLassoCoefficients =零（P + 1，numlambda）;BayesLassoCI95 =零（P + 1,2，numlambda）;fmseBayesLasso =零（numlambda，1）;BLCPlot =零（P + 1，numlambda）;%估计及预测RNG（10）;％用于重现对于J = 1：numlambda PriorMdl.Lambda =拉姆达（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 =查找（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
你选择了一个模型后，用整个数据集重新估计它。例如，以产生预测贝叶斯套索回归模型，创建一个先验模型和指定的收缩产生具有最小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实证
