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非常简单，假设样本yyy由信号β\betaβ与噪声ϵ\epsilonϵ组成：
y=β+ϵy = \beta + \epsilony=β+ϵ

其中噪声服从正态分布N(0,σ2In)N(0,\sigma^2I_n)N(0,σ2In)，β\betaβ的先验是Horseshoe prior，
βi∼N(0,λi2τ2)λi2∼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分布的第二象限部分沿yyy轴对折到第一象限得到的分布。假设β\betaβ非常稀疏，即{i:βi≠0}\{i:\beta_i \ne 0\}{i:βi=0}包含的元素个数远小于nnn，这个模型可以实现比较简单的稀疏信号复原。

函数定义

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函数定义的部分：
第一个输入yyy需要是一列向量，表示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/yi21/y_i^21/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)^2Tau = tau

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

第一步：给定λ,τ\lambda,\tauλ,τ，因为β\betaβ作为一个normal mean，先验又是正态分布，这正好就是共轭的情况，所以后验也是正态分布：
βi∣y,λi,τ∼N(λi2τ21+λi2τ2y,λi2τ21+λi2τ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τ2y,1+λi2τ2λi2τ2σ2)

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

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

第二步：给定β,λ\beta,\lambdaβ,λ，从π(τ∣β,λ,y)\pi(\tau|\beta,\lambda,y)π(τ∣β,λ,y)中采样，
τ=∣N(∑i=1nβiyi1+∑i=1nβi2,11+∑i=1nβi2)∣Gamma(n+12,1+∑i=1nβi22λi2)\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=1n2λi2βi2)∣N(1+∑i=1nβi2∑i=1nβiyi,1+∑i=1nβi21)∣

第三步：给定β,τ\beta,\tauβ,τ，从π(λ∣β,τ)\pi(\lambda|\beta,\tau)π(λ∣β,τ)中采样，
λi=N(βiλiyiGamma(1,1+λi22)+βi2λi2,1Gamma(1,1+λi22)+βi2λi2)\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βiyi,Gamma(1,21+λi2)+λi2βi21)

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函数相关推荐

Error in install.packages : cannot remove prior installation of package
Error in install.packages : cannot remove prior installation of package 目录 Error in install.packages ...
零基础自学R语言 1 R语言介绍 1.5 RStudio软件
零基础自学R语言文章目录零基础自学R语言 1 R语言介绍 1.5 RStudio软件 1.5.1 介绍 1.5.2 项目 1.5.3 帮助 1.5.4 使用技巧 1.5.4.1 使用历史命令 1. ...
R语言学习笔记——入门篇：第一章-R语言介绍
R语言 R语言学习笔记--入门篇:第一章-R语言介绍文章目录 R语言一.R语言简介 1.1.R语言的应用方向 1.2.R语言的特点二.R软件的安装 2.1.Windows/Mac 2.2.Lin ...
R语言实战-读书笔记（第1 章 R语言介绍）
*R语言实战所有学习笔记,如涉及侵权,请联系撤稿.* **标题号与书中标题号对应** R语言实战第1章 R语言介绍 1.2 R的获取与安装 R可以在CRAN(Comprehe ...
零基础自学R语言 1 R语言介绍 1.3 R扩展软件包的安装与管理
零基础自学R语言文章目录零基础自学R语言 1 R语言介绍 1.3 R扩展软件包的安装与管理 1.3.1 扩展包使用 1.3.2 安装 1.3.3 Github和BioConductor的扩展包 1 ...
初识R语言介绍以及常见的问题
R语言是用于统计分析,图形表示和报告的编程语言和软件环境. R语言由Ross Ihaka和Robert Gentleman在新西兰奥克兰大学创建,目前由R语言开发核心团队开发. R语言的核心是解释计算 ...
R语言从入门到精通Day1之【R语言介绍】
R语言开篇–R语言介绍开篇不再介绍R语言是如何下载和R语言的代码,如果您想真正的了解R,学习R,利用R做一些实际性的应用,不妨花点时间先了解一下当前数据科学的进展.了解R语言的历史和发展进程,R语言 ...
Making fast, good decisions with the FFTrees R package
"-[W]e are suspicious of rapid cognition. We live in a world that assumes that the quality of a ...
零基础自学R语言 1 R语言介绍 1.4 基本R软件的用法
零基础自学R语言文章目录零基础自学R语言 1 R语言介绍 1.4 基本R软件的用法 1.4.1 基本运行 1.4.2 项目目录 1 R语言介绍 1.4 基本R软件的用法 1.4.1 基本运行在M ...

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

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

Horseshoe prior的R package介绍：HS.normal.mean函数相关推荐

最新文章

热门文章