京等人。(2010)，“惩罚回归、标准误差和贝叶斯套索”，贝叶斯分析，5、2 ，表明对于计算套索预测的标准误差的统计有效方法可能没有达成共识。Tibshirani似乎同意（幻灯片 43）标准错误仍然是一个未解决的问题。

在可能有帮助的相关说明中，Tibshirani 及其同事提出了对套索的显着性检验。该论文是可用的，标题为“Asignificance test for the lasso”。可以在此处找到该论文的免费版本

贝叶斯 LASSO 是计算标准误差问题的唯一替代方案。标准误差在贝叶斯 LASSO 中自动计算...您可以使用吉布斯采样方案非常轻松地实现贝叶斯 LASSO...

贝叶斯 LASSO 需要将先验分布分配给模型的参数。在 LASSO 模型中，我们有目标函数和作为正则化参数。在这里，因为我们有 -norm 用于，所以需要一种特殊类型的先验分布，LAPLACE 分布是正态分布与指数分布的比例混合作为混合密度。基于每个参数的完整条件后验将被推导出来。$||\mathbf{y}-\mathbf{X}\boldsymbol{\beta}||_2^2 + \lambda||\boldsymbol{\beta}||_1$$\lambda$$\ell_1$$\boldsymbol{\beta}$

然后可以使用 Gibbs Sampling 来模拟链。参见 Park & Cassella (2008)，“贝叶斯套索” ，JASA，103，482。

常客 LASSO 存在三个固有的缺点：

1. 必须通过交叉验证或其他方式选择$\lambda$

2. 标准误差很难计算，因为 LARS 和其他算法会产生的点估计。$\boldsymbol{\beta}$

3. 手头问题的层次结构不能使用频率模型进行编码，这在贝叶斯框架中很容易。

Sandipan Karmakar 的回答告诉你该怎么做，这应该可以帮助你了解“如何”：

> library(monomvn)
>
> ## following the lars diabetes example
> data(diabetes)
> str(diabetes)
'data.frame':   442 obs. of  3 variables:
 $ x : AsIs [1:442, 1:10] 0.038075.... -0.00188.... 0.085298.... -0.08906.... 0.005383.... ...
      ..- attr(*, "dimnames")=List of 2
      .. ..$ : NULL
  .. ..$ : chr  "age" "sex" "bmi" "map" ...

 $ y : num  151 75 141 206 135 97 138 63 110 310 ...

[...]

> ## Bayesian Lasso regression
> reg_blas <- with(diabetes, blasso(x, y))
t=100, m=8
t=200, m=5
t=300, m=8
t=400, m=8
t=500, m=7
t=600, m=8
t=700, m=8
t=800, m=8
t=900, m=5
> 
> ## posterior mean beta (setting those with >50% mass at zero to exactly zero)
> (beta <- colMeans(reg_blas$beta) * (colMeans(reg_blas$beta != 0)  > 0.5))
      b.1       b.2       b.3       b.4       b.5       b.6       b.7       b.8 
   0.0000 -195.9795  532.7136  309.1673 -101.1288    0.0000 -196.4315    0.0000 
      b.9      b.10 
 505.4726    0.0000 
> 
> ## n x nsims matrix of realizations from the posterior predictive:
> post_pred_y <- with(reg_blas, X %*% t(beta))
> 
> ## predictions:
> y_pred <- rowMeans(post_pred_y)
> head(y_pred)
[1]  52.772443 -78.690610  24.234753   9.717777 -23.360369 -45.477199
> 
> ## sd of y:
> sd_y <- apply(post_pred_y, 1, sd)
> head(sd_y)
[1] 6.331673 6.756569 6.031290 5.236101 5.657265 6.150473
> 
> ## 90% credible intervals
> ci_y <- t(apply(post_pred_y, 1, quantile, probs=c(0.05, 0.95)))
> head(ci_y)
             5%       95%
[1,]  42.842535  62.56743
[2,] -88.877760 -68.47159
[3,]  14.933617  33.85679
[4,]   1.297094  18.01523
[5,] -32.709132 -14.13260
[6,] -55.533807 -35.77809