dl: Gibbs Sampler for Dirichlet-Laplace Prior

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

Description

dl aims to generate the posterior samples for Dirichlet-Laplace prior.

Usage

dl(x, y, hyper, a.prior, mcmc.n = 10000, thin = 1, eps = 1e-07, 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 hyperparameters (a1,b1) in the prior of sigma^2.
a.prior	The value of hyperparameter in the prior, which can be a value between (1/max(n,p), 1/2). By default, it is tuned through the MCMC process.
mcmc.n	Number of MCMC samples, by default, 10000.
thin	Thinning parameter of the chain. The default is 1, representing no thinning.
eps	Tolerance of convergence, by default, 1e-7.
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.

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

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

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

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

Examples

rho <- 0.5
# Number of predictors
p <- 25
# Number of observations
n <- 60
# Construct beta
n_nonzero <- 5
beta <- rep(0, p)
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
a.prior0 <- p/n
c.prior0 <- a.prior0/p
a.dl.prior0 <- 2 * c.prior0
mcmc.n <- 10000
fit.dl <- dl(x = x, y = y, hyper = list(a1 = 0.001, b1 = 0.001),
                  a.prior = a.dl.prior0, mcmc.n = mcmc.n, print = FALSE)

# Discard the early samples and get statistics of beta
burnIn <- 5000
beta.dl = fit.dl$beta[burnIn:mcmc.n, ]
mean.beta = apply(beta.dl, 2, mean)

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