horseshoe_regression: Gibbs Sampler for Horseshoe Regression

In dcbdan/s525: Implementation of common Bayesian models

Description

The horseshoe regression model is given by

[y_i|x_i, β, σ] \sim N(x_i^tβ, σ^2), i = 1,...,n ,

[β_j|σ, λ_j] \sim N(0, σ^2λ_j^2),

[λ_j|A] \sim C^{+}(0, A),

A \sim Uniform(0, 10),

p(σ^2) \propto \frac{1}{σ^2}.

The half-Cauchy parameter expansion is used; given by

[η_j|γ_j] \sim Gamma(\frac{1}{2}, γ_j),

[γ_j] \sim Gamma(\frac{1}{2}, \frac{1}{A^2}).

Let η_j = λ_j^{-2} , τ_A = A^{-2} , τ = \frac{1}{σ^2} and Λ = diag(η_1, ..., η_p). The full conditionals are given by:

[β|Y, X, η, τ] \sim \mathcal{N}( (X'X+Λ)^{-1}X'Y, τ^{-1}(X'X+Λ)^{-1}),

[η_j|β_j, γ_j, τ] \sim \mathrm{exp}(\frac{τβ_j^2}{2} + γ_j),

[γ_j|η_j, τ_A] \sim \mathrm{exp}(η_j + τ_A),

[τ_A|γ] \sim \mathrm{Gamma}( \frac{p - 1}{2}, ∑ γ_i)\mathrm{I}_{(\frac{1}{100}, ∞)},

[τ|Y, X, β, η] \sim \mathrm{Gamma}( \frac{n + p}{2}, \frac{(y-Xβ)'(y-Xβ) + β'Λβ}{2}).

Usage

horseshoe_regression <- function(
  Y,
  X,
  niter = 10000)

Arguments

Y	n by 1
X	n by p predictor matrix
niter	number of gibbs sampling iterations

Details

This function returns the generated parameters from the gibbs sampling markov chain.

Value

beta	An niter x p matrix
lambda	An niter x p matrix
sigma	An niter x 1 matrix

Examples

# Load the data
 prostate.data = "https://web.stanford.edu/~hastie/ElemStatLearn/datasets/prostate.data"
prostate = read.table(file = prostate.data, sep="", header = TRUE)
# Training data:
prostate_train = prostate[which(prostate$train),-10]
# Testing data:
prostate_test = prostate[which(!prostate$train),-10]
# Response:
y = prostate_train$lpsa
# Center and scale the data:
y = scale(y)
# And the predictors
X = scale(prostate_train[,names(prostate_train) != "lpsa"])

gibbs_hs <- horseshoe_regression(y, X, niter=10000)
shrinkage_regression_plot(gibbs_hs$beta, y, X)

dcbdan/s525 documentation built on May 19, 2019, 10:48 p.m.