主要内容

贝叶斯套索回归

打开脚本

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

套索回归是一种线性回归技术，它结合了正则化和变量选择。正则化通过降低回归系数的大小有助于防止过度拟合。套索回归的频率论观点不同于其他正则化技术，如岭回归，因为套索属性s值正好等于0，回归系数对应于不重要或冗余的预测值。

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

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

  • $y$是向量N响应。

  • $X$是一个矩阵P相应的观测预测变量。

  • $\beta$是向量P回归系数。

  • \ beta_0美元是拦截。

  • 一美元$是一个长度N一的列向量。

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

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

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

\λ美元是惩罚(或收缩)项，与其他参数不同，它不适合数据。必须指定的值\λ美元在估计之前，最好的做法是调整它。在指定值之后，$ f $相对于回归系数是最小的。得到的系数是套索估计值。有关拉索回归的频率观的更多细节，请参阅[191]．

对于本例，考虑为美国失业率创建一个预测线性模型。你需要一个一般化的模型。换句话说，你想通过去除所有冗余的预测器和所有与失业率无关的预测器来最小化模型的复杂性。

加载和预处理数据

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

负载数据模型

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

在同一图形中绘制所有系列，但在单独的子地块中绘制。

图;为j=1:size（DataTable，2）子图（4,4，j）plot（DataTable.Time，DataTable{：，j}）；title（DataTable.Properties.VariableNames（j））；结束

所有系列除了联邦基金,GS10,TB3MS,UNRATE呈指数趋势

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

haseExpoTrend=~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；title(“变量盒阴谋”）;

变量有相当相似的尺度。

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

  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项的线性公式]估计系数：估计当前值uuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuu（截距）3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.56 0.0 0.0 0.0.0.3.3.3.3.3.3.3.3.3.3.3.3.4.4.4.3.4.4.3.3.3.3.3.3.3.3.3.3.3.3.3.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7.7 7 7 7 7 7 7 7_Hoans-32.861 1.4507-22.652.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.1433.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-957：0.954 F-统计量与常数模型：292，p-值=2.26e-109 isInSig=1x4单元数组{'log_CPIAUCSL'}{'log_GS10'}{'log_M2SL'}{'FEDFUNDS'}

SLRMdlFull是一个LinearModel模型对象P-值表明，在5%显著性水平下，CPIAUCSL,GS10,M2SL,联邦基金这些都是无关紧要的变量。

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

idxreduced = ~ ismember (predictornames isInSig);XReduced = X (:, idxreduced);SLRMdlReduced = fitlm (XReduced y“VarNames”[predictornames (idxreduced) {“UNRATE”}));yFSLRF =预测(SLRMdlFull XF);yFSLRR =预测(SLRMdlReduced XF (:, idxreduced));fmsslrf =√(mean((yF - yFSLRF).^2)))
fmseSLRF = 0.6674

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

适合频率拉索回归模型

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

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

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

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

yFLasso = FitInfo。拦截+ XF * LassoBetaEstimates;mso =根号(yF - yFLasso .^2,1));

绘制回归系数相对于收缩值的大小。

hax = lassoPlot (LassoBetaEstimates FitInfo);L1Vals = hax.Children.XData;yyaxis正当h=绘图（L1VAL，Fmso，“线宽”2,“线型”,“——”）;传奇(h,“FMSE”,“位置”,“西南”)；伊莱贝尔(“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表格最佳套索截距61.154 log\u COE 0.042061 log\u cpaucsl 0 log\u GCE 0 log\u GDP 0 log\u GDPDEF 8.524 log\u GPDI 0 log\u GS10 1.6957 log\u hoans-18.937 log\u M1SL-1.3365 log\u M2SL 0.17948 log\u PCEC 0 FEDFUNDS 0 TB3MS-0.14459

频率拉索分析表明，变量CPIAUCSL,全球教育运动,国内生产总值,GPDI,PCEC,联邦基金要么无关紧要，要么多余。

拟合贝叶斯套索回归模型

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

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

在MATLAB®中实现贝叶斯套索回归时，请注意统计和机器学习工具箱之间的一些差异™ 作用套索和Econometrics Toolbox™对象lassoblm及其相关功能。

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

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

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

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

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

PriorMdl=bayeslm（p，“模型类型”,“套索”,“VarNames”，预测器名称）；表（PriorMdl.Lambda，“RowNames”，PriorMdl.VarNames）
ans=14x1表变量1截距0.01对数COE 1对数CPIUCSL 1对数GCE 1对数GDP 1对数GDPDEF 1对数GPDI 1对数GS10 1对数Hoans 1对数M1SL 1对数M2SL 1对数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);%有效样本量lambda=FitInfo.lambda*n./sqrt（FitInfo.MSE）；

考虑返回的正则化路径套索．在每个迭代中循环收缩值:

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

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

  3. 预测进入预测范围。

  4. 估计FMSE。

numlambda=numel（λ）；% 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,“显示”，false）；Bayeslasscoefficients（：，j）=Summary.Mean（1:，end-1））；BLCPlot（：，j）=Summary.Mean（1:，end-1））；Bayeslasscocients（：，j）=Summary.CI95（1:，end-1），；idx=Bayeslasscoci95（：，end-1，j）<=0；BLCPlot（idx，j）=0；yFBayesLasso=forecast（EstMdl，XF）；Fmsebayeslassaso（j）=sqrt（Mean）（Mean）（yF-yF））；Bayeso；结束

在每一次迭代,估计运行Gibbs采样器从完整条件中采样（请参见解析可处理后验概率)并返回经验性LBLM模型EstMdl，其中包含绘图和估算汇总表总结。您可以指定Gibbs采样器的参数，例如绘制数或细化方案。

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

绘制关于标准化L1惩罚的系数后均值（除截距外系数大小之和）。在同一个图上，但在右y轴上，绘制FMES。

L1Vals =总和(abs (BLCPlot(2:最终,:)),1)/ max(总和(abs (BLCPlot(2:最终,:)),1));图;情节(L1Vals BLCPlot(2:最终,:))包含(“L1”)；伊莱贝尔(“系数估计”)；YY轴正当h =情节(L1Vals fmseBayesLasso,“线宽”2,“线型”,“——”）;传奇(h,“FMSE”,“位置”,“西南”)；伊莱贝尔(“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])
2.2 2.2.2表数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计数据统计6 log_GPDI-2.0973-3.2168-0.96728 log_GS10 0.82485 0.22748 1.4237 log_Hoans-32。538-35.589-29.526 log_M1SL-3.0656-4.1499-2.0066 log_M2SL-1.1007-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,联邦基金是无关紧要的。

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

SLRR=0（p+1,1）；SLRR（[true idxreduced]）=SLRMdlReduced.coverties.Estimate；table（[SLRMdlFull.coverties.Estimate；fmseslr]，...[SLRR;fmseSLRR),...[最佳套索；最佳套索]，...[bestBayesLasso;fmsebestbayeslasso),“变化无常”,...{“SLRFull”“减少”“套索”“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=lambda（idx）；PriorMdl=bayeslm（p，“模型类型”,“套索”,“拉姆达”bestLambda,...“VarNames”, predictornames);PosteriorMdl =估计(PriorMdl [X;XF]、[y;yF]);
方法：lasso MCMC抽样10000次，观察次数：201个预测数：14 |平均Std CI95正分布----------------------------------------------------------------------截距| 85.9643 5.2710[75.544,96.407]1.000经验对数| 12.3574 2.1817[8.001,16.651]1.000经验对数| 0.14981.9150[-3.733，4.005]0.525经验对数|-7.1850 1.0556[-9.208，-5.101]0.000经验对数| 11.8648 3.9955[4.111，19.640]0.999经验对数| GDPDEF 8.8777 2.4221[4.033，13.745]1.000经验对数1240GPDI-2.84-10.690][0.194,1.197]0.997经验对数|-32.3356 1.4941[-35.276，-29.433]0.000经验对数|-2.2279 0.5043[-3.202，-1.238]0.000经验对数| 0.3617 0.9123[-1.438,2.179]0.652经验对数|-11.3018 3.0353[-17.318，-1.252 0.000经验对数1240.040]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建模财务风险的实用指南

下载电子书