使用Horseshoe 先验的Bayes回归及代码解析

最近的学习笔记，感觉质量比较差不足以给主站投稿，故发在这里让大家提提意见，也希望能有所帮助。

Horseshoe prior是一种稀疏bayes监督学习的方法。通过对模型参数的先验分布中加入稀疏特征，从而得到稀疏的估计。

horseshoe prior属于multivariate scale mixtures of normals的分布族。所以和其他常用的稀疏bayes学习方法，Laplacian prior, (Lasso), Student-t prior，非常类似。

在回归问题下的模型假设：
$y=X\beta+\epsilon, \epsilon\sim N(0,\sigma^2)$
$\beta_j \sim N(0,\sigma^2 \lambda_j^2 \tau^2)$
$\lambda_j \sim Half-Cauchy(0,1)$
可以额外加上以下超参数模型
$\tau \sim Half-Cauchy (0,1)$
$\sigma^2 \sim 1/\sigma^2$

其中 $\lambda_j$ 叫做local shrinkage parameter，局部压缩参数，对于不同的 $\theta_j$ 可以有不同的压缩系数， $\tau$ 叫做global shrinkage parameter

叫horseshoe，马蹄的原因大概是shrinkage weights的图像呈现马蹄形。

Carvalho, Polson, and Scott (2010) 给出了几个horseshoe estimator的相关性质以及和传统其他压缩先验方法的比较。

这篇文章给出了一个用于判断各种压缩先验性质的方法:
$E\left(\theta_{i} \mid y\right)=\int_{0}^{1}\left(1-\kappa_{i}\right) y_{i} p\left(\kappa_{i} \mid y\right) \mathrm{d} \kappa_{i}=\left\{1-E\left(\kappa_{i} \mid y\right)\right\} y_{i}$
在这种情况下, $E(\kappa_i|y)$ 就是压缩率。

而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$ 可以选择常数，截尾柯西分布和半柯西分布，默认的tau是1，同样 $\sigma$ 也可以选择固定值和Jeffreys prior ( $\sim 1/\sigma^2$ ) 。

if(p>n)
    algo=1
  else
    algo=2

在函数介绍中就有，当p>n时，使用的是Bhattacharya et.al. (2016)的算法，当n>p时，使用的是 Rue (2001)的算法对 $\beta$ 的full conditional posterior distribution进行抽样。

因为是Normal-Normal case，conditional on $\sigma,\tau,\lambda$ , $\beta$ 是Normal prior加Normal likelihood，所以对于Bayes回归问题，beta的全条件后验分布基本上可以写成:
$\beta \mid y, \lambda, \tau, \sigma \sim N\left(A^{-1} X^{\mathrm{T}} y, \sigma^{2} A^{-1}\right), \quad A=\left(X^{\mathrm{T}} X+\Lambda_{*}^{-1}\right), \quad \Lambda_{*}=\tau^{2} \operatorname{diag}\left(\lambda_{1}^{2}, \ldots, \lambda_{p}^{2}\right)$
这种形式，所以就是用Gibbs算法对 $\beta$ 抽样。

现在介绍在p>n情况下的算法。

来源于Bhattacharya. et.al. (2016)，

假设我们要从 $N(\mu,\Sigma)$ 中抽样，有
$\Sigma=\left(\Phi^{\mathrm{T}} \Phi+D^{-1}\right)^{-1}, \quad \mu=\Sigma \Phi^{\mathrm{T}} \alpha$
算法如下

抽 $u\sim N(0,D)$ , $\delta\sim N(0,I_n)$
令 $v=\Phi u+\delta$
解线性方程 $(\Phi D \Phi^T+I_n)w=(\alpha-v)$ 得到w
令 $\theta=u+D\Phi^T w$

算法的推导过程是：

我们目标是从 $N(\mu,\Sigma)$ 中进行抽样，由于
$\Sigma=(\Phi^T\Phi+D^{-1})^{-1}=D-D\Phi(\Phi D\Phi^T+I_n^{-1})^{-1}\Phi D$
那么第一块D，则可以通过抽 $u\sim N(0,D)$ 直接得到，也就是算法第4步中的 $$u$$ . 问题就是怎么得到第二块。

第二块等于
$D\Phi(\Phi D\Phi^T+I_n^{-1})^{-1}\Phi D$
$=D\Phi^T(\Phi D\Phi^T+I_n)^{-1}(\Phi D \Phi ^T+I_n)(\Phi D\Phi^T+I_n)^{-1}\Phi D$
因为
$\mu=\Sigma\Phi^T\alpha=(\Phi^T\Phi+D^{-1})^{-1}\Phi^T\alpha=(D\Phi^T-D\Phi^T(\Phi D\Phi^T+I_n)^{-1}\Phi D\Phi^T)\alpha$
$=D\Phi^T (\Phi D\Phi^T+I_n)^{-1}\alpha$
所以相当于只需要对 $N(0,(\Phi D\Phi^T+I_n))$ 加上 $\alpha$ 再做线性变换 $D\Phi^T\cdot (\Phi D\Phi^{T}+I_n)^{-1}$ 就能得到第二块的分布，也就对应算法过程3和4。

而为了得到 $N(0,(\Phi D\Phi^T+I_n))$ ，则只需要对 $u\sim N(0,D)$ 乘以 $\Phi$ 加 $\delta \sim N(0,I_n)$ ，对应算法过程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分解进行抽样，故不再赘述。

下一块则是抽取 $\lambda_j$ 的过程，这块算法使用了 Polson. et. al. (2014)的supplement material中，使用Damien et.al. (1999)算法构造的Slice sampler。

首先，我们可以写出 $\lambda_j$ 的后验分布：
$p(\lambda_j|y,\tau,\sigma)\propto \frac{2}{\pi(1+\lambda_j^2)}\cdot \frac{1}{\lambda_j}\cdot exp(-\frac{1}{\lambda_j^2}\frac{\beta_j^2}{2\tau^2\sigma^2} )$
做变量代换: $\eta_j=1/\lambda_j^2$ , $\mu_j=\beta_j/(\tau\sigma)$ . 那么 $\lambda_j$ 的full-conditional distribution转变为 $\eta_j$ 的 full conditional distribution,由变量代换公式:
$p(\eta_j|y,\tau,\sigma)\propto \frac{1}{\eta_j+1}\cdot exp(-\frac{\mu_j^2}{2}\eta_j)$
然后通过Damien et. al. (1999)的Slice sampler:

如果要生成密度函数为 $$f$$ 的随机数， $$f$$ 可以写成如下形式:
$f(x) \propto \pi(x) \prod_{i=1}^{N} l_{i}(x)$
其中 $\pi(x)$ 是已知形式的密度函数， $$l_i$$ 是非负可逆函数，也就是如果有 $$l_i(x)>u$$ ，那么可以反推出 $A^i_u=\{x:l_i(x)>u \}$ 。那么就可以通过添加辅助变量的方法生成服从 $$f$$ 的随机数。

令辅助变量 $u_i|x \sim unif (0,l_i(x))$ , 那么 $$U$$ 和 $$x$$ 的联合分布就是
$f\left(x, u_{1}, \ldots, u_{N}\right) \propto \pi(x) \prod_{i=1}^{N} I\left\{u_{i}<l_{i}(x)\right\}$
这样，关于X的边际分布就是 $$f(x)$$ 。我们构造Gibbs sampler，其中 $$u_i$$ 的full conditional distribution就是够早的uniform, 而X的full conditional distribution，则是分布 $\pi$ , 但是局限在区域 $A_{u}=\left\{x: l_{i}(x)>u_{i}, i=1, \ldots, N\right\}$ .

因为 $\pi$ 是已知好抽的分布， $$l_i$$ 是可逆函数，所以就把问题变成求区域和从好抽的分布中抽的问题。

应用在如上例子中，就是 $l(\eta_j)=\frac{1}{\eta_j+1}$ , $\pi(\eta_j)=exp(-\frac{\mu_j^2}{2}\eta_j)$ 是指数分布。

那算法：

从 $(0,l(\eta_j))$ 中抽 $$u$$ .
求区域 $A_u = \{\eta_j,l(\eta_j)>u\}$ , 也就是 $(0,\frac{1}{u}-1)$
在限制于 $(0,\frac{1}{u}-1)$ 的exp分布上抽 $\eta_j$ 。

这样就能得到 $\eta_j$ 。然后再逆变化回 $\lambda_j$ 就能得到 $\lambda_j$ 的样本，又因为 $\lambda_j$ 之间是独立的，所以这一步能向量化操作同时对所有的p个 $\lambda_j$ 抽样。

其中第三步的算法实现是：

由于局限在了 $(0,\frac{1}{u}-1)$ 上，而指数分布是单调的。就可以用一般的分布抽样的方法，从定义在 $(0, F(\frac{1}{u}-1))$ 的uniform抽，F是累积分布函数, 在这里就是 $1-exp(-\mu_j\cdot (\frac{1}{u}-1))$ , 然后再逆回去，得到目标随机数。

算法：

从 $(0,1-exp(-\mu_j\cdot (\frac{1}{u}-1)))$ 中抽 $$a$$
$\eta_j= -log(1-a)/\mu_j$

对应代码：

## 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);

剩下的更新 $\tau$ 的部分类似。但是由于 $\tau$ 是全局的参数，所以单个式子的形式就变成了叠乘的形式，所以就从exp分布变味了gamma分布。

The density of $\tau$

Gamma distribution $f(x)=\frac{\beta^{\alpha}}{\Gamma(\alpha)} x^{\alpha-1} e^{-\beta x}$

参考文献：

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.