r2d2marg: Gibbs Sampler for Marginal R2-D2

In yandorazhang/R2D2: Bayesian Linear Regression Using the R2-D2 Shrinkage Prior

Description

r2d2marg adopts the Gibbs sampling algorithm, aiming to obtain a sequence of posterior samples, which are approximately from the Marginal R2-D2 prior.

Usage

r2d2marg(x, y, hyper, mcmc.n = 10000, eps = 1e-07, thin = 1, print = TRUE)

Arguments

x	An n-by-p input matrix, each row of which is an observation vector.
y	An n-by-one vector, representing the response variable.
hyper	The values of the hyperparameters in the prior, i.e., the values of (a_pi, b, a1, b1). a_pi is the concentration parameter of the Dirichlet distribution, which controls the local shrinkage of phi. Smaller values of a_pi lead to most phi close to zero; while larger values of a_pi lead to a more uniform phi. b is the second shape parameter of beta prior for the R-squared. a1 and b1 are shape and scale parameters of Inverse Gamma distribution on sigma^2. hyper is set to be (1/(p^(b/2)n^(b/2)logn), 0.5, 0.001, 0.001) by default.
mcmc.n	Number of MCMC samples, by default, 10000.
eps	Tolerance of convergence, by default, 1e-7.
thin	Thinning parameter of the chain. The default is 1, representing no thinning.
print	Boolean variable determining whether to print the progress of the MCMC iterations, i.e., which iteration the function is currently on. The default is TRUE.

Value

A list containing the following components:

• beta: Matrix (mcmc.n * p) of posterior samples for beta.

• sigma2: Vector (mcmc.n) of posterior samples for sigma^2.

• psi: Matrix (mcmc.n * p) of posterior samples for psi.

• w: Vector (mcmc.n) of posterior samples for the total prior probability w.

• xi: Vector (mcmc.n) of posterior samples for xi.

Examples

rho <- 0.5
# Number of predictors
p <- 25
# Number of observations
n <- 60
# Construct beta
n_nonzero <- 5
beta <- rep(0, p)
set.seed(1)
beta[11:(10 + n_nonzero)] <- stats::rt(n_nonzero, df = 3) * sqrt(0.5/(3 * n_nonzero/2))
# Construct x
sigma <- 1
times <- 1:p
H <- abs(outer(times, times, "-"))
V <- sigma * rho^H
x <- mvtnorm::rmvnorm(n, rep(0, p), V)
x <- scale(x, center = TRUE, scale = FALSE)
# Construct y
y <- x %*% beta + stats::rnorm(n)

# Gibbs sampling
mcmc.n <- 10000
fit.new <- r2d2marg(x = x, y = y, mcmc.n = mcmc.n, print = FALSE)

# Discard the early samples
burnIn <- 5000
beta.new <- fit.new$beta[burnIn:10000, ]

yandorazhang/R2D2 documentation built on Nov. 23, 2020, 7:05 a.m.