主要内容

贝叶斯套索回归

打开脚本

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

套索回归是一个线性回归技术,结合正则化和变量选择。正则化有助于防止过度拟合减少回归系数的大小。套索回归的频率论的观点不同于其他正则化技术,例如,岭回归,因为套索属性的值是0到回归系数对应于无关紧要的或多余的预测因子。

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

$ $ y = \ beta_0 1 _n”+ X \β+ \ varepsilon。$ $

  • y美元是一个向量的n响应。

  • X美元是一个矩阵p相应的预测变量。

  • β\美元是一个向量的p回归系数。

  • \ beta_0美元是拦截。

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

  • \ varepsilon美元是一个随机向量的iid 0均值和方差的高斯干扰\σ^ 2美元。

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

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

\λ美元点球(或收缩)项,与其他参数不同,不适合数据。你必须指定一个值\λ美元估计之前,和一个最佳实践是调整它。在您指定一个值,$ f $最小化对回归系数。拉索产生的系数估计。更多细节的频率论的观点套索回归,明白了[207]。

对于这个示例,考虑创建一个预测线性模型对美国失业率。你想要一个模型,概括了。换句话说,你想最小化模型复杂性通过删除所有多余的预测和预测与失业率是不相关的。

加载和数据预处理

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

负载Data_USEconModel

数据集包括MATLAB®时间表DataTimeTable,其中包含14个变量测量从1947年第一季度到2009年一季度;UNRATE美国失业率。更多细节,回车描述在命令行中。

情节都在同一个图系列,但在不同的次要情节。

图;为j = 1:尺寸(DataTimeTable, 2)次要情节(4 4 j)情节(DataTimeTable.Time, DataTimeTable {: j});标题(DataTimeTable.Properties.VariableNames (j));结束

所有系列除了FEDFUNDS,GS10,TB3MS,UNRATE似乎一个指数的趋势。

对数变换应用到这些变量有一个指数的趋势。

hasexpotrend = ~ ismember (DataTimeTable.Properties.VariableNames…(“FEDFUNDS”“GD10”“TB3MS”“UNRATE”]);DataTimeTableLog = varfun (@log DataTimeTable,“数据源”,…DataTimeTable.Properties.VariableNames (hasexpotrend));DataTimeTableLog = [DataTimeTableLog DataTimeTable (:, DataTimeTable.Properties.VariableNames (~ hasexpotrend)));

DataTimeTableLog是一个时间表DataTimeTable,但这些变量与一个指数趋势在对数尺度。

系数有较大震级往往占主导地位的点球套索回归目标函数。因此,它是重要的变量有类似的规模当你实现套索回归。比较变量的尺度DataTimeTableLog通过绘制他们的箱形图在同一轴。

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

的变量都相当类似的鳞片。

确定时间序列模型推广,如何遵循这个过程:

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

  2. 符合模型来估计样本。

  3. 使用拟合模型来预测响应预测地平线。

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

  5. FMSE最低的选择模型。

创建评估和预测样本变量的响应和预测数据。指定一个预测地平线的4年(16人)。

跳频= 16;y = DataTimeTableLog。UNRATE(1:(end - fh)); yF = DataTimeTableLog.UNRATE((end - fh + 1):end); isresponse = DataTimeTable.Properties.VariableNames ==“UNRATE”;X = DataTimeTableLog{1:(结束- fh), ~ isresponse};XF = DataTimeTableLog{(结束- fh + 1):最终,~ isresponse};p =大小(X, 2);%的预测数量predictornames = DataTimeTableLog.Properties.VariableNames (~ isresponse);

符合简单线性回归模型

适合简单线性回归模型来估计样本。指定变量名(响应变量的名称必须是最后一个元素)。提取p值。识别无关紧要的系数通过使用5%水平的意义。

SLRMdlFull = fitlm (X, y,“VarNames”DataTimeTableLog.Properties.VariableNames) slrPValue = SLRMdlFull.Coefficients.pValue;isInSig = SLRMdlFull。CoefficientNames (slrPValue > 0.05)
SLRMdlFull =线性回归模型:UNRATE ~(与14条款13预测线性公式)估计系数:估计SE tStat pValue说____ __________(拦截)88.014 5.3229 16.535 2.6568 e-37 log_COE log_GCE log_CPIAUCSL 0.0027764 7.1111 2.3426 3.0356 -3.4032 2.4611 -1.3828 0.16853 -5.63 1.1581 -4.8615 2.6245 e-06 log_GDP log_GDPDEF 0.00030659 14.522 3.9406 3.6851 11.926 2.9298 4.0704 7.1598 e-05 log_GPDI -2.5939 0.54603 -4.7504 4.2792 e-06 log_GS10 log_HOANBS 0.051565 0.57467 0.29313 1.9605 -32.861 1.4507 -22.652 2.3116 e-53 log_M1SL -3.3023 0.53185 -6.209 3.914 e-09 log_M2SL FEDFUNDS log_PCEC 0.18647 -1.3845 1.0438 -1.3264 -7.143 3.4302 -2.0824 0.038799 0.044542 0.047747 0.93287 0.3522 TB3MS -0.1577 0.064421 -2.448 0.015375的观测数量:185年,错误自由度:171根均方误差:0.312平方:0.957,调整平方:0.954 f统计量与常数模型:292年,假定值= 2.26 e - 109 isInSig = 1 x4单元阵列{‘log_CPIAUCSL} {‘log_GS10} {‘log_M2SL} {' FEDFUNDS '}

SLRMdlFull是一个LinearModel模型对象。的p值表明,5%水平的意义,CPIAUCSL,GS10,M2SL,FEDFUNDS是无关紧要的变量。

改装没有无关紧要的变量的线性模型的预测数据。充分利用和减少模型,预测失业率预测地平线,然后估计fms。

idxreduced = ~ ismember (predictornames isInSig);XReduced = X (:, idxreduced);SLRMdlReduced = fitlm (XReduced y“VarNames”[predictornames (idxreduced) {“UNRATE”}));yFSLRF =预测(SLRMdlFull XF);yFSLRR =预测(SLRMdlReduced XF (:, idxreduced));fmseSLRF =√意味着(yF - yFSLRF) ^ 2)) fmseSLRR(12(平均((yF - yFSLRR) ^ 2)。)
fmseSLRF fmseSLRR = 0.6105 = 0.6674

减少模型的FMSE小于FMSE完整的模型,这表明减少模型概括更好。

适合拉索频率论的回归模型

适合套索对数据的回归模型。使用默认的正规化的道路。套索标准化数据默认情况下。因为类似尺度变量,不规范。返回的信息符合沿着正规化的道路。

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

默认情况下,套索符合套索回归模型通过使用一系列的值的100倍\λ美元。因此,LassoBetaEstimates是一个13 -到- 100矩阵的回归系数的估计。行对应于预测变量X和列对应的值\λ美元。FitInfo是一个结构,它包括的值\λ美元(FitInfo.Lambda),每个合适的均方误差(FitInfo.MSE)和估计的拦截(FitInfo.Intercept)。

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

yFLasso = FitInfo。拦截+ XF * LassoBetaEstimates;fmseLasso =√意味着(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);LassoBetaEstimates (:, idx)];表(bestLasso,“RowNames”,(“拦截”predictornames])
ans = 14 x1表bestLasso _____拦截61.154 log_COE 0.042061 log_CPIAUCSL log_GCE 0 log_GDP log_GDPDEF 8.524 log_GPDI 0 0.17948 -1.3365 -18.937 1.6957 log_GS10 log_HOANBS log_M1SL log_M2SL log_PCEC FEDFUNDS 0 TB3MS -0.14459

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

符合贝叶斯套索回归模型

贝叶斯观点的套索回归,回归系数的先验分布的拉普拉斯(双指数),意思是0和规模σ\ / \λ美元,在那里\λ美元是固定的收缩参数和美元\σ^ 2 \ sim搞笑(A, B)美元。更多细节,请参阅lassoblm。

与拉索的频率论的观点回归,如果你增加\λ美元,那么由此产生的模型的稀疏单调增加。但是,与频率论的套索,贝叶斯套索微不足道或冗余模式系数不是0。相反,估计后验分布和无关紧要的或冗余如果0周边地区具有较高的密度。

在MATLAB®实现贝叶斯套索回归时,要注意几个差异的统计和机器学习工具箱™函数套索和计量经济学工具箱™对象lassoblm及其相关功能。

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

  • 不像套索,lassoblm不规范的预测数据。然而,您可以提供不同的收缩值对于每个系数。该功能有助于维护系数估计的可解释性。

  • lassoblm一个收缩值适用于每个系数;它不接受正规化的道路套索。

因为lassoblm框架是适合进行贝叶斯套索每系数回归为一个收缩值,您可以使用一个for循环执行套索正规化的道路。

通过使用之前创建一个贝叶斯套索回归模型bayeslm。指定数量的预测变量和变量名。显示每个系数存储在默认的收缩值λ模型的属性。

PriorMdl = bayeslm (p,“ModelType”,“套索”,“VarNames”,predictornames);表(PriorMdl.Lambda,“RowNames”PriorMdl.VarNames)
__拦截ans Var1 = 14 x1表0.01 log_COE 1 log_CPIAUCSL log_GCE 1 log_GDP log_GDPDEF 1 log_GPDI log_GS10 1 log_HOANBS log_M1SL 1 log_M2SL 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) =√意味着(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 =找到(fmseBayesLasso > 0.68, 1);fmsebestbayeslasso = fmseBayesLasso (idx);bestBayesLasso = BayesLassoCoefficients (:, idx);bestCI95 = BayesLassoCI95 (:,:, idx);表(bestBayesLasso bestCI95,“RowNames”,(“拦截”predictornames])
______________________拦截ans = 14 x2表bestBayesLasso bestCI95 * * * 90.587 79.819 101.62 log_COE 6.5591 1.8093 11.239 log_CPIAUCSL log_GCE -2.2971 -6.9019 1.7057 -5.1707 -7.4902 -2.8583 log_GDP log_GDPDEF 9.8839 2.3848 17.672 10.281 5.1677 15.936 log_GPDI log_GS10 -2.0973 -3.2168 -0.96728 0.82485 0.22748 1.4237 log_HOANBS log_M1SL -32.538 -35.589 -29.526 -3.0656 -4.1499 -2.0066 log_M2SL log_PCEC -1.1007 -3.1243 0.75844 -3.3342 -9.6496 1.8482 FEDFUNDS TB3MS -0.15682 -0.29385 -0.022145 0.043672 -0.056192 0.14254

有一种方法可以确定一个变量无关紧要或冗余是通过检查是否在其相应的系数0 95%可信区间。使用这个标准,CPIAUCSL,M2SL,PCEC,FEDFUNDS是无关紧要的。

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

SLRR = 0 (p + 1);SLRR([真idxreduced]) = SLRMdlReduced.Coefficients.Estimate;表([SLRMdlFull.Coefficients.Estimate;fmseSLRR),…[SLRR;fmseSLRR),…[bestLasso;fmsebestlasso),…[bestBayesLasso;fmsebestbayeslasso),“VariableNames”,…{“SLRFull”“SLRReduced”“套索”“BayesLasso”},…“RowNames”,(PriorMdl.VarNames;MSE的])
ans = 15 x4表SLRFull SLRReduced套索BayesLasso ________ __________ ________ __________拦截log_COE 90.587 88.014 88.327 61.154 7.1111 10.854 0.042061 6.5591 log_CPIAUCSL -3.4032 0 0 -2.2971 log_GCE log_GDP 14.522 16.701 -5.1707 -5.63 - -6.1958 0 0 9.8839 log_GPDI log_GDPDEF 11.926 9.1106 8.524 10.281 -2.5939 -2.6963 -2.0973 log_GS10 0.57467 0 1.6957 0.82485 log_HOANBS log_M2SL log_M1SL -32.538 -32.861 -33.782 -18.937 -3.3023 -2.8099 -1.3365 -3.0656 -1.3845 -3.3342 -7.143 - -14.207 0.17948 - -1.1007 log_PCEC 0 0 0 0.043672 0.044542 FEDFUNDS TB3MS MSE 0.61048 0.61048 0.79425 -0.15682 -0.1577 -0.078449 -0.14459 0.69639

在你选择了一个模型,使用整个数据集重新评估它。例如,创建一个预测贝叶斯套索回归模型,创建一个先验模型和指定收缩产生与最小FMSE最简单的模型,然后估计使用整个数据集。

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

另请参阅

||

相关的话题