使用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, Chakraborty, and Mallick (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情况下的算法。

来源于@Bhattacharya2016 ，

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

1. 抽 \(u\sim N(0,D)\) , \(\delta\sim N(0,I_n)\)
2. 令 \(v=\Phi u+\delta\)
3. 解线性方程\((\Phi D \Phi^T+I_n)w=(\alpha-v)\) 得到w
4. 令\(\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\). 问题就是怎么得到第二块。

第二块等于 \[ \begin{align*} &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 \end{align*} \] 因为 \[ \begin{align*} \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 \end{align*} \] 所以相当于只需要对\(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, Scott, and Windle (2014) 的supplement material中，使用 Damien, Wakefield, and Walker (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)\) 是指数分布。

那算法：

1. 从\((0,l(\eta_j))\)中抽\(u\) .

2. 求区域 \(A_u = \{\eta_j,l(\eta_j)>u\}\), 也就是\((0,\frac{1}{u}-1)\)

3. 在限制于\((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))\), 然后再逆回去，得到目标随机数。

算法：

1. 从\((0,1-exp(-\mu_j\cdot (\frac{1}{u}-1)))\) 中抽 \(a\)
2. \(\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}\)