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$由信號$\beta$與噪聲$?$組成：
$$y=\beta +?$$

其中噪聲服從正態分布$N(0,{\sigma }^2{I}_n)$，$\beta$的先驗是Horseshoe prior，
$$\begin{gathered} {\beta }_i～N(0,{\lambda }_i^2{\tau }^2) \\ {\lambda }_i^2～C^{+}(0,1),{\tau }^2～C^{+}(0,1) \end{gathered}$$

$C^{+}(0,1)$是標準half-Cauchy分布，大家可以把它理解成是標準Cauchy分布的第二象限部分沿$y$軸對折到第一象限得到的分布。假設$\beta$非常稀疏，即$\{i:{\beta }_i=0\}$包含的元素個數遠小于$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$需要是一列向量，表示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$的取值，在理論研究的模擬實驗中通常假設$\sigma =1$，所以第五個輸入的默認值為1；如果選擇Jeffrey先驗，也就是$\pi (\sigma )\propto 1/\sigma$，那么模型中的每一個參數都有先驗，并且先驗不依賴于超參，這樣的模型叫做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$。

${\lambda }_i$的初值為$1/y_i^2$，這是常規的初值設置方式；$\tau$的初值就是函數的第三步輸入，實際上感覺這個code關于$\tau$的初值有一些設計，

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

但感覺執行之后Tau還是等于函數第三個輸入。

第一步：給定$\lambda ,\tau$，因為$\beta$作為一個normal mean，先驗又是正態分布，這正好就是共軛的情況，所以后驗也是正態分布：
$${\beta }_i|y,{\lambda }_i,\tau ～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)$$

這四行代碼就是在從這個分布中采樣

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

第二步：給定$\beta ,\lambda$，從$\pi (\tau |\beta ,\lambda ,y)$中采樣，
$$\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})}}$$

第三步：給定$\beta ,\tau$，從$\pi (\lambda |\beta ,\tau )$中采樣，
$${\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}})$$

---

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 = 1Beta = yTau = tauif (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)^2Tau = taufor (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)) * yBeta = 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)^2b = 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^2utau = stats::runif(1, 0, 1/(1 + et))ubt_1 = 1ubt_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/ThetaV2 = 1/stats::rgamma(n, 1, rate = (Lambda^2 + 1)/2)num1 = V2 * Theta^2den = 1 + num1s = sqrt(V2/den)m = (num1/den) * ZLambda = stats::rnorm(n, m, s)if (t > burn) {BetaSave[t - burn, ] = BetaTauSave[t - burn] = TauSigma2Save[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) }

總結

以上是生活随笔為你收集整理的Horseshoe prior的R package介绍：HS.normal.mean函数的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。