Horseshoe prior的R package介绍：HS.normal.mean函数

最近做的一些事情需要和Horseshoe prior对比，所以一直在看Horseshoe的一些资料。上周做了一点simulation发现Horseshoe在normal mean model上的表现还挺不错的，所以打算扒一下horseshoe这个包里面的HS.normal.means这个函数看看那几位搞Horseshoe的大牛是怎么写的。这里贴一个那个包的说明文本：https://cran.r-project.org/web/packages/horseshoe/horseshoe.pdf

本文末尾附了HS.normal.means的源代码，懒得自己去找的可以拉到最后看看。

---

首先Normal mean model非常简单，假设样本 y y y由信号 β \beta β与噪声 ϵ \epsilon ϵ组成：
y = β + ϵ y = \beta + \epsilon y=β+ϵ

其中噪声服从正态分布 N ( 0 , σ 2 I n ) N(0,\sigma^2I_n) N(0,σ2In​)， β \beta β的先验是Horseshoe prior，
β i ∼ N ( 0 , λ i 2 τ 2 ) λ i 2 ∼ C + ( 0 , 1 ) , τ 2 ∼ C + ( 0 , 1 ) \beta_i \sim N(0,\lambda_i^2\tau^2) \\ \lambda_i^2 \sim C^+(0,1),\tau^2 \sim C^+(0,1) βi​∼N(0,λi2​τ2)λi2​∼C+(0,1),τ2∼C+(0,1)

C + ( 0 , 1 ) C^+(0,1) C+(0,1)是标准half-Cauchy分布，大家可以把它理解成是标准Cauchy分布的第二象限部分沿 y y y轴对折到第一象限得到的分布。假设 β \beta β非常稀疏，即 { i : β i ≠ 0 } \{i:\beta_i \ne 0\} {i:βi​​=0}包含的元素个数远小于 n n n，这个模型可以实现比较简单的稀疏信号复原。

---

函数定义

function (y, method.tau = c("fixed", "truncatedCauchy", "halfCauchy"), 
  tau = 1, method.sigma = c("fixed", "Jeffreys"), Sigma2 = 1, 
  burn = 1000, nmc = 5000, alpha = 0.05) 

这是HS.normal.means函数定义的部分：
第一个输入 y y y需要是一列向量，表示observation；
第二个输入method.tau表示 τ \tau τ的先验，这里可以选择 τ \tau τ作为固定常数、truncated cauchy或者halfcauchy；因为 τ \tau τ在模型中是global variance component，它的作用是在对每一个signal做belief update的时候可以从其他signal那里borrow information，所以如果是选择fixed，一般也是先用MLE（在horseshoe这个包中可以用HS.MMLE这个函数）估计出了 τ \tau τ的值，然后再作为固定值在这里输入（这种操作叫empirical Bayesian）；
第三个输入 τ \tau τ默认值为1，如果第二个输入选择的是fixed，那么第三个输入就是输入 τ \tau τ的固定取值，如果第二个输入选择的是那两种prior之一，那么这个输入不会被采用；
第四个输入是噪声的标准差 σ \sigma σ的取值方法，如果选择fixed，那么就需要第五个输入确认 σ \sigma σ的取值，在理论研究的模拟实验中通常假设 σ = 1 \sigma=1 σ=1，所以第五个输入的默认值为1；如果选择Jeffrey先验，也就是 π ( σ ) ∝ 1 / σ \pi(\sigma) \propto 1/\sigma π(σ)∝1/σ，那么模型中的每一个参数都有先验，并且先验不依赖于超参，这样的模型叫做full Bayesian；
第六个输入设置burn-in sample的数量，第七个输入设置chain的长度，因为这个函数大框架是Metropolis-Hastings，所以也是需要这两个参数的；
因为这个函数输出包含 β \beta β的置信区间，所以第八个输入alpha的作用是控制置信水平（1-alpha是置信水平）。

---

输出值定义

result <- list(BetaHat = BetaHat, LeftCI = left.points, 
    RightCI = right.points, BetaMedian = BetaMedian, Sigma2Hat = Sigma2Hat, 
    TauHat = TauHat, BetaSamples = BetaSave, TauSamples = TauSave, 
    Sigma2Samples = Sigma2Save)
  return(result)

这些输出分别是 β \beta β的后验均值、1-alpha置信区间下界、1-alpha置信区间上界、 β \beta β的MAP估计、 σ \sigma σ的后验均值、 τ \tau τ的后验均值、 β \beta β的chain、 τ \tau τ的chain、 σ \sigma σ的chain

---

核心算法

我主要是想看 σ \sigma σ固定、 τ \tau τ用halfcauchy的情况，所以就以此为例看看这个函数的源代码。首先，核心算法的大框架是blocked Sampling，第一步采样 β ∣ λ , τ \beta|\lambda,\tau β∣λ,τ；第二步采样 τ ∣ λ , β \tau|\lambda,\beta τ∣λ,β；第三步采样 λ ∣ τ , β \lambda|\tau,\beta λ∣τ,β。

λ i \lambda_i λi​的初值为 1 / y i 2 1/y_i^2 1/yi2​，这是常规的初值设置方式； τ \tau τ的初值就是函数的第三步输入，实际上感觉这个code关于 τ \tau τ的初值有一些设计，

  if (method.tau != "fixed") {
    Tau = max(sum(abs(y) >= sqrt(2 * Sigma2 * log(n)))/n, 
      1/n)
  }
  Lambda = 1/abs(y)^2
  Tau = tau

但感觉执行之后Tau还是等于函数第三个输入。

第一步：给定 λ , τ \lambda,\tau λ,τ，因为 β \beta β作为一个normal mean，先验又是正态分布，这正好就是共轭的情况，所以后验也是正态分布：
β i ∣ y , λ i , τ ∼ N ( λ i 2 τ 2 1 + λ i 2 τ 2 y , λ i 2 τ 2 1 + λ i 2 τ 2 σ 2 ) \beta_i|y,\lambda_i,\tau \sim N(\frac{\lambda_i^2\tau^2}{1+\lambda_i^2\tau^2}y,\frac{\lambda_i^2\tau^2}{1+\lambda_i^2\tau^2}\sigma^2) βi​∣y,λi​,τ∼N(1+λi2​τ2λi2​τ2​y,1+λi2​τ2λi2​τ2​σ2)

这四行代码就是在从这个分布中采样

    a = (Tau^2) * (Lambda^2)
    s = sqrt(Sigma2 * a/(1 + a))
    m = (a/(1 + a)) * y
    Beta = stats::rnorm(n, m, s)

第二步：给定 β , λ \beta,\lambda β,λ，从 π ( τ ∣ β , λ , y ) \pi(\tau|\beta,\lambda,y) π(τ∣β,λ,y)中采样，
τ = ∣ N ( ∑ i = 1 n β i y i 1 + ∑ i = 1 n β i 2 , 1 1 + ∑ i = 1 n β i 2 ) ∣ G a m m a ( n + 1 2 , 1 + ∑ i = 1 n β i 2 2 λ i 2 ) \tau = \frac{|N(\frac{\sum_{i=1}^n \beta_i y_i}{1+\sum_{i=1}^n \beta_i^2},\frac{1}{1+\sum_{i=1}^n \beta_i^2})|}{\sqrt{Gamma(\frac{n+1}{2},1+\sum_{i=1}^n \frac{\beta_i^2}{2\lambda_i^2})}} τ=Gamma(2n+1​,1+∑i=1n​2λi2​βi2​​) ​∣N(1+∑i=1n​βi2​∑i=1n​βi​yi​​,1+∑i=1n​βi2​1​)∣​

第三步：给定 β , τ \beta,\tau β,τ，从 π ( λ ∣ β , τ ) \pi(\lambda|\beta,\tau) π(λ∣β,τ)中采样，
λ i = N ( β i λ i y i G a m m a ( 1 , 1 + λ i 2 2 ) + β i 2 λ i 2 , 1 G a m m a ( 1 , 1 + λ i 2 2 ) + β i 2 λ i 2 ) \lambda_i = N(\frac{\frac{\beta_i}{\lambda_i}y_i}{Gamma(1,\frac{1+\lambda_i^2}{2})+\frac{\beta_i^2}{\lambda_i^2}},\frac{1}{Gamma(1,\frac{1+\lambda_i^2}{2})+\frac{\beta_i^2}{\lambda_i^2}}) λi​=N(Gamma(1,21+λi2​​)+λi2​βi2​​λi​βi​​yi​​,Gamma(1,21+λi2​​)+λi2​βi2​​1​)

---

function (y, method.tau = c("fixed", "truncatedCauchy", "halfCauchy"), 
  tau = 1, method.sigma = c("fixed", "Jeffreys"), Sigma2 = 1, 
  burn = 1000, nmc = 5000, alpha = 0.05) 
{
  method.tau = match.arg(method.tau)
  method.sigma = match.arg(method.sigma)
  n <- length(y)
  BetaSave = matrix(0, nmc, n)
  TauSave = rep(0, nmc)
  Sigma2Save = rep(0, nmc)
  Lambda2 = rep(1, n)
  nu = rep(1, n)
  xi = 1
  Beta = y
  Tau = tau
  if (method.sigma == "Jeffreys") {
    Sigma2 = 0.95 * stats::var(y)
  }
  if (method.tau != "fixed") {
    Tau = max(sum(abs(y) >= sqrt(2 * Sigma2 * log(n)))/n, 
      1/n)
  }
  Lambda = 1/abs(y)^2
  Tau = tau
  for (t in 1:(nmc + burn)) {
    if (t%%1000 == 0) {
      print(t)
    }
    a = (Tau^2) * (Lambda^2)
    s = sqrt(Sigma2 * a/(1 + a))
    m = (a/(1 + a)) * y
    Beta = stats::rnorm(n, m, s)
    Theta = Beta/(Lambda)
    if (method.tau == "halfCauchy") {
      G = 1/sqrt(stats::rgamma(1, (n + 1)/2, rate = (1 + 
        sum(Theta^2))/2))
      Z = y/(Theta * Lambda)
      a = (Lambda * Theta)^2
      b = sum(a)
      s2 = 1/(1 + b)
      m = {
        s2
      } * sum(a * Z)
      Delta = stats::rnorm(1, m, sqrt(s2))
      Tau = abs(Delta) * G
    }
    if (method.tau == "truncatedCauchy") {
      tempt = sum(Beta^2/Lambda^2)/(2 * Sigma2)
      et = 1/Tau^2
      utau = stats::runif(1, 0, 1/(1 + et))
      ubt_1 = 1
      ubt_2 = min((1 - utau)/utau, n^2)
      Fubt_1 = stats::pgamma(ubt_1, (n + 1)/2, scale = 1/tempt)
      Fubt_2 = stats::pgamma(ubt_2, (n + 1)/2, scale = 1/tempt)
      ut = stats::runif(1, Fubt_1, Fubt_2)
      et = stats::qgamma(ut, (n + 1)/2, scale = 1/tempt)
      Tau = 1/sqrt(et)
    }
    if (method.sigma == "Jeffreys") {
      Sigma2 = 1/stats::rgamma(1, n, rate = sum((y - Beta)^2)/2 + 
        sum(Beta^2/(Tau^2 * Lambda^2))/2)
    }
    Z = y/Theta
    V2 = 1/stats::rgamma(n, 1, rate = (Lambda^2 + 1)/2)
    num1 = V2 * Theta^2
    den = 1 + num1
    s = sqrt(V2/den)
    m = (num1/den) * Z
    Lambda = stats::rnorm(n, m, s)
    if (t > burn) {
      BetaSave[t - burn, ] = Beta
      TauSave[t - burn] = Tau
      Sigma2Save[t - burn] = Sigma2
    }
  }
  BetaHat = colMeans(BetaSave)
  BetaMedian = apply(BetaSave, 2, stats::median)
  TauHat = mean(TauSave)
  Sigma2Hat = mean(Sigma2Save)
  left <- floor(alpha * nmc/2)
  right <- ceiling((1 - alpha/2) * nmc)
  BetaSort <- apply(BetaSave, 2, sort, decreasing = F)
  left.points <- BetaSort[left, ]
  right.points <- BetaSort[right, ]
  result <- list(BetaHat = BetaHat, LeftCI = left.points, 
    RightCI = right.points, BetaMedian = BetaMedian, Sigma2Hat = Sigma2Hat, 
    TauHat = TauHat, BetaSamples = BetaSave, TauSamples = TauSave, 
    Sigma2Samples = Sigma2Save)
  return(result)
}