- 已编辑
使用Horseshoe 先验的Bayes回归及代码解析
最近的学习笔记,感觉质量比较差不足以给主站投稿,故发在这里让大家提提意见,也希望能有所帮助。
Horseshoe prior是一种稀疏bayes监督学习的方法。通过对模型参数的先验分布中加入稀疏特征,从而得到稀疏的估计。
horseshoe prior属于multivariate scale mixtures of normals的分布族。所以和其他常用的稀疏bayes学习方法,Laplacian prior, (Lasso), Student-t prior,非常类似。
在回归问题下的模型假设:
可以额外加上以下超参数模型
其中 叫做local shrinkage parameter,局部压缩参数,对于不同的可以有不同的压缩系数, 叫做global shrinkage parameter
叫horseshoe,马蹄的原因大概是shrinkage weights的图像呈现马蹄形。
Carvalho, Polson, and Scott (2010) 给出了几个horseshoe estimator的相关性质以及和传统其他压缩先验方法的比较。
这篇文章给出了一个用于判断各种压缩先验性质的方法:
在这种情况下, 就是压缩率。
而horseshoe prior在这个指标下就相对很优秀:对于很小的参数,压缩效果明显,对于相对比较大的参数,几乎不压缩。
R包horseshoe提供了horseshoe estimator的计算函数,而本文主要目的就是对horseshoe.R中的算法进行整合和解析。
horseshoe = function(y,X, method.tau = c("fixed", "truncatedCauchy","halfCauchy"), tau = 1,
method.sigma = c("fixed", "Jeffreys"), Sigma2 = 1,
burn = 1000, nmc = 5000, thin = 1, alpha = 0.05)
函数的定义,对于 可以选择常数,截尾柯西分布和半柯西分布,默认的tau是1,同样 也可以选择固定值和Jeffreys prior () 。
if(p>n)
algo=1
else
algo=2
在函数介绍中就有,当p>n时,使用的是Bhattacharya et.al. (2016)的算法,当n>p时,使用的是 Rue (2001)的算法对 的full conditional posterior distribution进行抽样。
因为是Normal-Normal case,conditional on , 是Normal prior加Normal likelihood,所以对于Bayes回归问题,beta的全条件后验分布基本上可以写成:
这种形式,所以就是用Gibbs算法对 抽样。
现在介绍在p>n情况下的算法。
来源于Bhattacharya. et.al. (2016),
假设我们要从中抽样,有
算法如下
- 抽 ,
- 令
- 解线性方程 得到w
- 令
算法的推导过程是:
我们目标是从 中进行抽样,由于
那么第一块D,则可以通过抽 直接得到,也就是算法第4步中的. 问题就是怎么得到第二块。
第二块等于
因为
所以相当于只需要对 加上再做线性变换 就能得到第二块的分布,也就对应算法过程3和4。
而为了得到,则只需要对 乘以 加, 对应算法过程2。
这样的算法过程就对应代码:
lambda_star=tau*lambda
U=as.numeric(lambda_star^2)*t(X)
## step 1 ##
u=stats::rnorm(l2,l0,lambda_star)
v=X%*%u + stats::rnorm(n)
## step 2 ##
v_star=solve((X%*%U+I_n),((y/sqrt(sigma_sq))-v))
Beta=sqrt(sigma_sq)*(u+U%*%v_star)
当p<n,也就是algo==2时则是用经典的Cholesky分解进行抽样,故不再赘述。
下一块则是抽取 的过程,这块算法使用了 Polson. et. al. (2014)的supplement material中,使用Damien et.al. (1999)算法构造的Slice sampler。
首先,我们可以写出 的后验分布:
做变量代换: , . 那么 的full-conditional distribution转变为的 full conditional distribution,由变量代换公式:
然后通过Damien et. al. (1999)的Slice sampler:
如果要生成密度函数为 的随机数, 可以写成如下形式:
其中 是已知形式的密度函数, 是非负可逆函数,也就是如果有 ,那么可以反推出。那么就可以通过添加辅助变量的方法生成服从 的随机数。
令辅助变量 , 那么 和 的联合分布就是
这样,关于X的边际分布就是。我们构造Gibbs sampler,其中 的full conditional distribution就是够早的uniform, 而X的full conditional distribution,则是分布, 但是局限在区域.
因为是已知好抽的分布,是可逆函数,所以就把问题变成求区域和从好抽的分布中抽的问题。
应用在如上例子中,就是 , 是指数分布。
那算法:
从中抽 .
求区域 , 也就是
在限制于 的exp分布上抽。
这样就能得到。然后再逆变化回 就能得到 的样本,又因为 之间是独立的,所以这一步能向量化操作同时对所有的p个抽样。
其中第三步的算法实现是:
由于局限在了 上,而指数分布是单调的。就可以用一般的分布抽样的方法,从定义在的uniform抽,F是累积分布函数, 在这里就是 , 然后再逆回去,得到目标随机数。
算法:
- 从 中抽
对应代码:
## update lambda_j's in a block using slice sampling ##
eta = 1/(lambda^2)
upsi = stats::runif(p,0,1/(1+eta))
tempps = Beta^2/(2*sigma_sq*tau^2)
ub = (1-upsi)/upsi
# now sample eta from exp(tempv) truncated between 0 & upsi/(1-upsi)
Fub = 1 - exp(-tempps*ub) # exp cdf at ub
Fub[Fub < (1e-4)] = 1e-4; # for numerical stability
up = stats::runif(p,0,Fub)
eta = -log(1-up)/tempps
lambda = 1/sqrt(eta);
剩下的更新 的部分类似。但是由于是全局的参数,所以单个式子的形式就变成了叠乘的形式,所以就从exp分布变味了gamma分布。
The density of
Gamma distribution
参考文献:
Bhattacharya, Anirban, Antik Chakraborty, and Bani K. Mallick. 2016. “Miscellanea fast sampling with Gaussian scale mixture priors in high-dimensional regression.” Biometrika 103 (4): 985–91. https://doi.org/10.1093/biomet/asw042.
Carvalho, Carlos M, Nicholas G Polson, and James G Scott. 2010. “The horseshoe estimator for sparse signals.” Biometrika 97 (2): 465–80. https://doi.org/10.1093/biomet/asq017.
Damien, Paul, Jon Wakefield, and Stephen Walker. 1999. “Gibbs sampling for Bayesian non-conjugate and hierarchical models by using auxiliary variables.” Journal of the Royal Statistical Society. Series B: Statistical Methodology 61 (2): 331–44. https://doi.org/10.1111/1467-9868.00179.
Polson, Nicholas G., James G. Scott, and Jesse Windle. 2014. “The Bayesian bridge.” Journal of the Royal Statistical Society. Series B: Statistical Methodology 76 (4): 713–33. https://doi.org/10.1111/rssb.12042.