MHLS {EAinference}

Metropolis-Hastings lasso sampler under a fixed active set.

Description

Metropolis-Hastings sampler to simulate from the sampling distribution of lasso given a fixed active set.

Usage

MHLS(X, PE, sig2, lbd, weights = rep(1, ncol(X)), B0, S0, A = which(B0 !=
  0), tau = rep(1, ncol(X)), niter = 2000, burnin = 0, PEtype = "coeff",
  updateS.itv = 1, verbose = FALSE, ...)

Arguments

X

predictor matrix.

PE, sig2, lbd

parameters of target distribution. (point estimate of beta or E(y) depends on PEtype, variance estimate of error and lambda).

weights

weight vector with length p(the number of covariates). Default is weights = rep(1, p).

B0

numeric vector with length p. Initial value of lasso estimator.

S0

numeric vector with length p. Initial value of subgradients. If not given, this will be generated in a default way.

A

numeric vector. Active coefficient index. Every active coefficient index in B0 must be included. Default is A = which(B0 != 0).

tau

numeric vector with length p. Standard deviation of proposal distribution for each coefficient.

niter

integer. The number of iterations. Default is niter = 2000

burnin

integer. The length of burin-in periods. Default is burnin = 0

PEtype

Type of PE which is needed to characterize the target distribution. Users can choose either "coeff" or "mu".

updateS.itv

integer. Update subgradients every updateS.itv iterations. Set this value larger than niter if one wants to skip updating subgradients.

verbose

logical. If true, print out the progress step.

...

complementary arguments.

• FlipSA : optional parameter. This has to be a subset of active set, A. If the index is not listed in FlipSA, the sign of coefficients which correspond to the listed index will remain fixed. The default is FlipSA=A

• SFindex : optional parameter. subgradient index for the free coordinate.

• randomSFindex : logical. If true, resample SFindex every updateSF.itv iterations.

• updateSF.itv : integer. In every updateSF.itv iterations, randomize SFindex.

Details

Given appropriate initial value, provides Metropolis-Hastings samples under the fixed active set.
From the initial values, B0 and S0, MHLS draws beta and subgrad samples. In every iteration, given t-th iteration values, t-th beta and t-th subgrad, a new set of proposed beta and subgradient is sampled. We either accept the proposed sample and use that as (t+1)-th iteration values or reuse t-th iteration values.
See Zhou(2014) for more details.

Value

MHLS returns an object of class "MHLS". The functions summary.MHLS and plot.MHLS provide a brief summary and generate plots.

beta

lasso samples.

subgrad

subgradient samples.

acceptHistory

numbers of acceptance and proposal.

niter, burnin, PE, type

same as function arguments.

References

Zhou, Q. (2014), "Monte Carlo simulation for Lasso-type problems by estimator augmentation," Journal of the American Statistical Association, 109, 1495-1516.

Examples

#-------------------------
# Low dim
#-------------------------
set.seed(123)
n <- 10
p <- 5
X <- matrix(rnorm(n * p), n)
Y <- X %*% rep(1, p) + rnorm(n)
sigma2 <- 1
lbd <- .37
weights <- rep(1, p)
LassoResult <- lassoFit(X = X, Y = Y, lbd = lbd, type = "lasso", weights = weights)
B0 <- LassoResult$B0
S0 <- LassoResult$S0
MHLS(X = X, PE = rep(0, p), sig2 = 1, lbd = 1,
     weights = weights, B0 = B0, S0 = S0, niter = 50, burnin = 0,
     PEtype = "coeff")
MHLS(X = X, PE = rep(0, n), sig2 = 1, lbd = 1,
     weights = weights, B0 = B0, S0 = S0, niter = 50, burnin = 0,
     PEtype = "mu")

#-------------------------
# High dim
#-------------------------
set.seed(123)
n <- 5
p <- 10
X <- matrix(rnorm(n*p),n)
Y <- X %*% rep(1,p) + rnorm(n)
weights <- rep(1,p)
LassoResult <- lassoFit(X = X,Y = Y,lbd = lbd, type = "lasso", weights = weights)
B0 <- LassoResult$B0
S0 <- LassoResult$S0
MHLS(X = X, PE = rep(0, p), sig2 = 1, lbd = 1,
     weights = weights, B0 = B0, S0 = S0, niter = 50, burnin = 0,
     PEtype = "coeff")
MHLS(X = X, PE = rep(0, n), sig2 = 1, lbd = 1,
     weights = weights, B0 = B0, S0 = S0, niter = 50, burnin = 0,
     PEtype = "mu")