Mhorseshoe: Approximate Algorithm for Horseshoe Prior

Description

Mhorseshoe is a package for a high-dimensional Bayesian linear modeling algorithm using a horseshoe prior. This package provides two different algorithm functions : exact_horseshoe, approx_horseshoe. approx_horseshoe is version that can lower the computational cost than the existing horseshoe estimator through the approximate MCMC algorithm in the case of p >> N for p predictors and N observations. You can see examples of the use of the two algorithms through the vignette, browseVignettes("Mhorseshoe").

Author(s)

Maintainer: Kang Mingi leehuimin115@g.skku.edu

Authors:

• Lee Kyoungjae leekjstat@gmail.com

Run approximate MCMC algorithm for horseshoe prior

Description

The approximate MCMC algorithm for the horseshoe prior

Usage

approx_horseshoe(
  y,
  X,
  burn = 1000,
  iter = 5000,
  auto.threshold = TRUE,
  tau = 1,
  s = 0.8,
  sigma2 = 1,
  alpha = 0.05,
  ...
)

Arguments

Details

This function implements the approximate algorithm introduced in Section 2.2 of Johndrow et al. (2020) and the method proposed in this package, which improves computation speed when p >> N. The approximate algorithm introduces a threshold and uses only a portion of the total p columns for matrix multiplication, reducing the computational cost compared to the existing MCMC algorithms for the horseshoe prior. The "auto.threshold" argument determines whether the threshold used in the algorithm will be selected by the adaptive method proposed in this package. For more information, browseVignettes("Mhorseshoe").

Value

References

Johndrow, J., Orenstein, P., & Bhattacharya, A. (2020). Scalable Approximate MCMC Algorithms for the Horseshoe Prior. In Journal of Machine Learning Research, 21, 1-61.

Examples

# Making simulation data.
set.seed(123)
N <- 200
p <- 100
true_beta <- c(rep(1, 10), rep(0, 90))

X <- matrix(1, nrow = N, ncol = p) # Design matrix X.
for (i in 1:p) {
  X[, i] <- stats::rnorm(N, mean = 0, sd = 1)
}

y <- vector(mode = "numeric", length = N) # Response variable y.
e <- rnorm(N, mean = 0, sd = 2) # error term e.
for (i in 1:10) {
  y <- y + true_beta[i] * X[, i]
}
y <- y + e

# Run with auto.threshold set to TRUE
result1 <- approx_horseshoe(y, X, burn = 0, iter = 100,
                            auto.threshold = TRUE)

# Run with fixed custom threshold
result2 <- approx_horseshoe(y, X, burn = 0, iter = 100,
                            auto.threshold = FALSE, threshold = 1/(5 * p))

# posterior mean
betahat <- result1$BetaHat

# Lower bound of the 95% credible interval
leftCI <- result1$LeftCI

# Upper bound of the 95% credible interval
RightCI <- result1$RightCI

Run exact MCMC algorithm for horseshoe prior

Description

The exact MCMC algorithm for the horseshoe prior introduced in section 2.1 of Johndrow et al. (2020).

Usage

exact_horseshoe(
  y,
  X,
  burn = 1000,
  iter = 5000,
  tau = 1,
  s = 0.8,
  sigma2 = 1,
  alpha = 0.05,
  ...
)

Arguments

Details

The exact MCMC algorithm introduced in Section 2.1 of Johndrow et al. (2020) is implemented in this function. This algorithm uses a blocked-Gibbs sampler for (\tau, \beta, \sigma^2), where the global shrinkage parameter \tau is updated by an Metropolis-Hastings algorithm. The local shrinkage parameter \lambda_{j},\ j = 1,2,...,p is updated by the rejection sampler.

Value

References

Johndrow, J., Orenstein, P., & Bhattacharya, A. (2020). Scalable Approximate MCMC Algorithms for the Horseshoe Prior. In Journal of Machine Learning Research, 21, 1-61.

Examples

# Making simulation data.
set.seed(123)
N <- 50
p <- 100
true_beta <- c(rep(1, 10), rep(0, 90))

X <- matrix(1, nrow = N, ncol = p) # Design matrix X.
for (i in 1:p) {
  X[, i] <- stats::rnorm(N, mean = 0, sd = 1)
}

y <- vector(mode = "numeric", length = N) # Response variable y.
e <- rnorm(N, mean = 0, sd = 2) # error term e.
for (i in 1:10) {
  y <- y + true_beta[i] * X[, i]
}
y <- y + e

# Run exact_horseshoe
result <- exact_horseshoe(y, X, burn = 0, iter = 100)

# posterior mean
betahat <- result$BetaHat

# Lower bound of the 95% credible interval
leftCI <- result$LeftCI

# Upper bound of the 95% credible interval
RightCI <- result$RightCI

rejection sampler to update local shrinkage parameters lambda.

Description

rejection sampler to update local shrinkage parameters lambda.

Usage

rejection_sampler(eps, a, b)