主要内容

贝叶斯套索回归

打开脚本

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

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

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

/ * * * * * * * * * * * * * * *

  • 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 $相对于回归系数是最小的。得到的系数是套索估计值。有关拉索回归的频率观的更多细节，请参阅[191]．

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

加载和预处理数据

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

加载Data_USEconModel

数据集包括MATLAB®时间表数据表，其中包含从1947年Q1到2009年第1季度测量的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．..［“菲德福德斯”“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 =;h.XTickLabelRotation = 45;标题('变量箱图'）;

变量具有相似相似的尺度。

要确定时间序列模型的泛化程度，请遵循以下步骤:

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

  2. 将模型适合估计样本。

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

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

  5. 选择FMSE最低的型号。

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

跳频= 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 =大小(X, 2);%预测数predictornames = DataTableLog.Properties.VariableNames (~ isresponse);

拟合简单线性回归模型

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

SLRMdlFull = fitlm (X, y,'varnames'slrPValue = slrmdlfull .系数. pvalue;isInSig = SLRMdlFull。CoefficientNames (slrPValue > 0.05)
SLRMdlFull =线性回归模型:UNRATE ~[13个预测因子中14项的线性公式]Estimate 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_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.75044.2792e-06 log_GS10 0.57467 0.29313 1.9605 0.051565 log_HOANBS -32.861 1.4507 -22.652 2.3160 e-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.038799F-statistic vs. constant model: 292, p-value = 2.26e-109 isig = 1x4 cell array {'log_CPIAUCSL'} {'log_GS10'} {'log_M2SL'} {'FEDFUNDS'}

SLRMdlFull是一个LinearModel模型对象。的p- 值得提示，以5％的意义，cpiasucsl.，GS10，M2SL， 和FEDFUNDS是无关紧要的变量。

将预测数据中不含不显著变量的线性模型重新拟合。利用完整模型和简化模型，在预测范围内预测失业率，然后估计FMSEs。

idxreduced = ~ ismember (predictornames isInSig);XReduced = X (:, idxreduced);SLRMdlReduced = fitlm (XReduced y'varnames'，[预测archames（Idxreduced）{“UNRATE”}));yFSLRF =预测(SLRMdlFull XF);yFSLRR =预测(SLRMdlReduced XF (:, idxreduced));fmsslrf =√(mean((yF - yFSLRF).^2)))
FMSESLRF = 0.6674 FMSSSLRR = 0.6105

减少模型的FMSE小于完整模型的FMSE，这表明减少的模型更好地推广。

适合频率拉索回归模型

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

[LassoBetaEstimates, FitInfo] =套索(X, y,“标准化”，错误的，．..'predictornames', predictornames);

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

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

yflasso = fitinfo.intercept + xf * lassobetaaestimates;fmselasso = sqrt（平均值（（yf  -  yflasso）。^ 2,1））;

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

hax = lassoPlot (LassoBetaEstimates FitInfo);L1Vals = hax.Children.XData;yyaxis正确的h =情节(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）;LassobetaSimates（：，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

频率拉索分析表明，变量cpiasucsl.，全球教育运动，国内生产总值，GPDI，PCEC， 和FEDFUNDS要么无关紧要，要么多余。

拟合贝叶斯拉索回归模型

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

和拉索回归的频率观一样，如果你增加\λ美元，则所得模型的稀疏性单调增大。然而，不像频率套索，贝叶斯套索有不重要或冗余系数模式，不完全是0。相反，估计具有后验分布，如果0附近的区域密度很高，则不显著或冗余。

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

  • lassoblm是对象框架的一部分吗套索是一个函数。对象框架简化了计量计量工作流程。

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

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

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

通过使用创建贝叶斯套索回归贝叶斯林．指定预测器变量的数量和变量名称。显示存储在中的每个系数的默认收缩值lambda.模型的财产。

PriorMdl = bayeslm (p,“ModelType”，'套索'，'varnames', predictornames);表(PriorMdl。λ,“RowNames”PriorMdl.VarNames)
ans = 14x1 table Var1 ____ Intercept 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 fefunds 1 TB3MS 1

Priormdl.是一个lassoblm模型对象。lassoblm收缩属性1对于除拦截之外的每个系数，它具有收缩0．01．这些默认值在大多数应用程序中可能用处不大;最佳实践是使用多种值进行试验。

考虑返回的正则化路径套索．将收缩值膨胀一倍$ n / \ sqrt（mse_j）$,在那里MSE_j美元是套索跑的妈妈$ j $，$ j $= 1通过收缩值的数量。

ismissing =任何(isnan (X), 2);n =总和(~ ismissing);有效样本量λ= FitInfo.Lambda * n /√(FitInfo.MSE);

考虑返回的正则化路径套索．通过每次迭代的收缩值循环：

  1. 估计回归系数的后验分布和干扰方差，给定收缩和数据。因为预测器的规模是接近的，attribute相同的收缩到每个预测器，并attribute收缩0．01拦截。保存后验平均值和95%可信区间。

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

  3. 预测预测地平线。

  4. 估计FMSE。

numlambda =元素个数(λ);％preallocateBayeslassocoefficients =零（P + 1，Numlambda）;Bayeslassoci95 =零（P + 1,2，Numlambda）;FMSEBAYESLASSO =零（NUMLAMBDA，1）;blcplot =零（p + 1，numlambda）;估算和预测rng (10);%的再现性为了j = 1：numlambda primdl.lambda = 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) =√(mean((yF - yFBayesLasso).^2)));结尾

在每一次迭代,估计运行吉布斯采样器，从完整的条件句(见分析易易行的后索)，并返回empiricalblm模型EstMdl，其中包含绘图和评估汇总表总结．您可以为吉布斯采样器指定参数，例如绘制或细化方案的数量。

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

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

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

当归一化L1惩罚超过0.1时，模型趋向于更好地泛化。为了平衡最小的FMSE和模型复杂性，选择最简单的模型，使FMSE接近0.68。

Idx = find（fmsebayeslasso> 0.68,1）;fmsebestbayeslasso = fmsebayeslasso（Idx）;Bestbayeslasso =贝叶斯类科科菲素（：，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-32.538 -35.589 -29.526 log_M1SL -3.0656 -4.1499 -2.0066 log_M2SL -1.1007 -3.1243 0.5844 log_PCEC -3.3342 -9.6496 1.8482

确定变量是否不重要或冗余的一种方法是检查0是否在其相应系数的95%可信区间内。使用这个标准,cpiasucsl.，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 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_GS100.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 FEDFUNDS

选择模型后，使用整个数据集重新估算它。例如，要创建预测贝叶斯套索回归模型，创建先前的模型，并指定具有最小FMSE的最简单模型的收缩，然后使用整个数据集来估计它。

bestLambda =λ(idx);PriorMdl = bayeslm (p,“ModelType”，'套索'，“λ”，Bestlambda，．..'varnames', predictornames);posteriormdl =估计（priormdl，[x; xf]，[y; yf]）;
方法:lasso MCMC sampling with 10000 draw观测数: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,0.000 Empirical log_GDPDEF | 11.8648 3.9955 [4.111, 19.640] 0.999 Empirical log_GDPDEF | 8.8777 2.4221 [4.033, 13.745] 1.000 Empirical log_GPDI | -2.5884 0.5294 [-3.628， -1.553] 0.000 Empirical log_GS10 | 0.910 0.2577 [0.194, 1.197] 0.997 Empirical log_HOANBS | -32.3356 1.4941 [-35.276，0.000 Empirical log_M1SL | -2.2279 0.5043 [-3.202， -1.238] 0.000 Empirical log_M2SL | 0.3617 0.9123 [-1.438, 2.179] 0.652 Empirical log_PCEC | -11.3018 3.0353 [-17.318， -5.252] 0.000 Empirical FEDFUNDS | -0.0132 0.0490 [-0.109, 0.082] 0.392 Empirical TB3MS | -0.1128 0.0660 [-0.244, 0.016] 0.042 Empirical Sigma2 | 0.1186 0.0122 [0.097，0.145) 1.000经验

另请参阅

||

相关的话题

OuthoMetrics工具箱文档

金宝app

用MATLAB建模财务风险的实用指南

用MATLAB建模财务风险的实用指南

下载电子书