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

2021-04-08 11:29:55  阅读：6  来源： 互联网

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

y = β + ϵ y = \beta + \epsilon y=β+ϵ

β 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)


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


λ 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


β 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)


τ = ∣ 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​)∣​

λ 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)
}