• 统计学
  • 使用Horseshoe 先验的Bayes回归及代码解析

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

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

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

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

在回归问题下的模型假设:
y=Xβ+ϵ,ϵN(0,σ2)y=X\beta+\epsilon, \epsilon\sim N(0,\sigma^2)
βjN(0,σ2λj2τ2)\beta_j \sim N(0,\sigma^2 \lambda_j^2 \tau^2)
λjHalfCauchy(0,1)\lambda_j \sim Half-Cauchy(0,1)
可以额外加上以下超参数模型
τHalfCauchy(0,1)\tau \sim Half-Cauchy (0,1)
σ21/σ2\sigma^2 \sim 1/\sigma^2

其中λj\lambda_j 叫做local shrinkage parameter,局部压缩参数,对于不同的θj\theta_j可以有不同的压缩系数,τ\tau 叫做global shrinkage parameter

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

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

这篇文章给出了一个用于判断各种压缩先验性质的方法:
E(θiy)=01(1κi)yip(κiy)dκi={1E(κiy)}yiE\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(κiy)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 (1/σ2\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的全条件后验分布基本上可以写成:
βy,λ,τ,σN(A1XTy,σ2A1),A=(XTX+Λ1),Λ=τ2diag(λ12,,λp2)\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(μ,Σ)N(\mu,\Sigma)中抽样,有
Σ=(ΦTΦ+D1)1,μ=ΣΦTα\Sigma=\left(\Phi^{\mathrm{T}} \Phi+D^{-1}\right)^{-1}, \quad \mu=\Sigma \Phi^{\mathrm{T}} \alpha
算法如下

  1. uN(0,D)u\sim N(0,D) , δN(0,In)\delta\sim N(0,I_n)
  2. v=Φu+δv=\Phi u+\delta
  3. 解线性方程(ΦDΦT+In)w=(αv)(\Phi D \Phi^T+I_n)w=(\alpha-v) 得到w
  4. θ=u+DΦTw\theta=u+D\Phi^T w

算法的推导过程是:

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

第二块等于
DΦ(ΦDΦT+In1)1ΦDD\Phi(\Phi D\Phi^T+I_n^{-1})^{-1}\Phi D
=DΦT(ΦDΦT+In)1(ΦDΦT+In)(ΦDΦT+In)1Φ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
因为
μ=ΣΦTα=(ΦTΦ+D1)1ΦTα=(DΦTDΦT(ΦDΦT+In)1ΦDΦT)α\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ΦT(ΦDΦT+In)1α=D\Phi^T (\Phi D\Phi^T+I_n)^{-1}\alpha
所以相当于只需要对N(0,(ΦDΦT+In))N(0,(\Phi D\Phi^T+I_n)) 加上α\alpha再做线性变换 DΦT(ΦDΦT+In)1D\Phi^T\cdot (\Phi D\Phi^{T}+I_n)^{-1} 就能得到第二块的分布,也就对应算法过程3和4。

而为了得到N(0,(ΦDΦT+In))N(0,(\Phi D\Phi^T+I_n)),则只需要对uN(0,D)u\sim N(0,D) 乘以Φ\PhiδN(0,In)\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分解进行抽样,故不再赘述。

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

首先,我们可以写出λj\lambda_j 的后验分布:
p(λjy,τ,σ)2π(1+λj2)1λjexp(1λj2βj22τ2σ2)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} )
做变量代换: ηj=1/λj2\eta_j=1/\lambda_j^2, μj=βj/(τσ)\mu_j=\beta_j/(\tau\sigma). 那么λj\lambda_j 的full-conditional distribution转变为ηj\eta_j的 full conditional distribution,由变量代换公式:
p(ηjy,τ,σ)1ηj+1exp(μj22ηj)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:

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

令辅助变量 uixunif(0,li(x))u_i|x \sim unif (0,l_i(x)), 那么UUxx 的联合分布就是
f(x,u1,,uN)π(x)i=1NI{ui<li(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)f(x)。我们构造Gibbs sampler,其中uiu_i 的full conditional distribution就是够早的uniform, 而X的full conditional distribution,则是分布π\pi, 但是局限在区域Au={x:li(x)>ui,i=1,,N}A_{u}=\left\{x: l_{i}(x)>u_{i}, i=1, \ldots, N\right\}.

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

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

那算法:

  1. (0,l(ηj))(0,l(\eta_j))中抽uu .

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

  3. 在限制于(0,1u1)(0,\frac{1}{u}-1) 的exp分布上抽ηj\eta_j

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

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

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

算法:

  1. (0,1exp(μj(1u1)))(0,1-exp(-\mu_j\cdot (\frac{1}{u}-1))) 中抽 aa
  2. ηj=log(1a)/μj\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 distributionf(x)=βαΓ(α)xα1eβxf(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.

    帖子的公式环境貌似不太对劲,改了半天,如何用正则来快速实现这个过程,最重要的是把 $公式$ 改成 \(公式\)。。。查了一下对应好像可以用s/(\s)\$(.+?)\$/\1<code class="katex-escape">\(\2\)</code>/sg来换,我试了一下在R里:

    library(dplyr)
    a=readtext::readtext("./horseshoe-bayes.Rmd")
    a%>% gsub("(\s)\$(.+?)\$","\1<code class="katex-escape">\(\2\)</code>")

    但是提示

    Error: '\s' is an unrecognized escape in character string starting ""(\s"

    感觉根本上的用法就不对,不太懂正则该怎么用...

      wglaive

      R里的\还得转义

      不过这个涉及多个$的情况下我估计正则很难准确判断两个$如何匹配,大概率会出问题,还是手动改吧

        写得很好,如果能首先上个图直观比较一下各个先验 shrinkage 的不同,最后跑个模拟/实际数据,比较一下几个方法的性能就更好了。如果能再加上 stan 实现更佳,毕竟换个先验都要手推一遍 Gibbs 抽样还是略痛苦了一些……

          nan.xiao 谢谢建议!我试着按这个思路改改....虽然不知道啥时候弄得出来。手推gibbs这个确实是...虽然这就是目前我在做的主要工作.....不过不太懂stan,得去研究一下....由于拖延症(情绪管理?)失败了好几次了.....头大
          Cloud2016 立字据(还是算了...)orz,感觉会拖很久的感觉哈哈哈哈哈哈哈哈哈哈哈

            wglaive 这样的话 建议可以写成系列文章 比如第一篇讲理论 第二篇讲计算推导 第三篇讲实现 每篇专注一件事 这样把大的任务分解 设立明确的milestone 比较容易做成

            3 年 后

            @wglaive
            想問你
            以下code是\tau \sim Half-Cauchy (0,1)
            果我想將\tau \sim Half-Cauchy (0,1)改成\tau \sim Half-Cauchy (0,0.9) code要怎麼修改

            if(method.tau == "truncatedCauchy"){
              tempt = sum((Beta/lambda)^2)/(2*sigma_sq)
              et = 1/tau^2
              utau = stats::runif(1,0,1/(1+et))
              ubt_1=1
              ubt_2 = min((1-utau)/utau,p^2)
              Fubt_1 = stats::pgamma(ubt_1,(p+1)/2,scale=1/tempt)
              Fubt_2 = stats::pgamma(ubt_2,(p+1)/2,scale=1/tempt)
              #Fubt = max(Fubt,1e-8) # for numerical stability
              ut = stats::runif(1,Fubt_1,Fubt_2)
              et = stats::qgamma(ut,(p+1)/2,scale=1/tempt)
              tau = 1/sqrt(et)
            }