CovEst.soft {CovTools}

Covariance Estimation via Soft Thresholding

Description

Soft Thresholding method for covariance estimation takes off-diagonal elements z of sample covariance matrix and applies

h_{\tau}(z) = \textrm{sgn}(z)(|z|-\tau)_{+}

where \textrm{sgn}(z) is a sign of the value z, and (x)_+ = \textrm{max}(x,0). If thr is rather a vector of regularization parameters, it applies cross-validation scheme to select an optimal value.

Usage

CovEst.soft(X, thr = 0.5, nCV = 10, parallel = FALSE)

Arguments

X

an (n\times p) matrix where each row is an observation.

thr

user-defined threshold value. If it is a vector of regularization values, it automatically selects one that minimizes cross validation risk.

nCV

the number of repetitions for 2-fold random cross validations for each threshold value.

parallel

a logical; TRUE to use half of available cores, FALSE to do every computation sequentially.

Value

a named list containing:

S

a (p\times p) covariance matrix estimate.

CV

a dataframe containing vector of tested threshold values(thr) and corresponding cross validation scores(CVscore).

References

Antoniadis A, Fan J (2001). “Regularization of Wavelet Approximations.” Journal of the American Statistical Association, 96(455), 939–967. ISSN 0162-1459, 1537-274X.

Donoho DL, Johnstone IM, Kerkyacharian G, Picard D (1995). “Wavelet Shrinkage: Asymptopia?” Journal of the Royal Statistical Society. Series B (Methodological), 57(2), 301–369. ISSN 00359246.

Examples

## generate data from multivariate normal with Identity covariance.
pdim <- 5
data <- matrix(rnorm(10*pdim), ncol=pdim)

## apply 4 different schemes
#  mthr is a vector of regularization parameters to be tested
mthr <- exp(seq(from=log(0.1),to=log(10),length.out=10))

out1 <- CovEst.soft(data, thr=0.1)  # threshold value 0.1
out2 <- CovEst.soft(data, thr=1)    # threshold value 1
out3 <- CovEst.soft(data, thr=10)   # threshold value 10
out4 <- CovEst.soft(data, thr=mthr) # automatic threshold checking

## visualize 4 estimated matrices
gcol <- gray((0:100)/100)
opar <- par(no.readonly=TRUE)
par(mfrow=c(2,2), pty="s")
image(out1$S[,pdim:1], col=gcol, main="thr=0.1")
image(out2$S[,pdim:1], col=gcol, main="thr=1")
image(out3$S[,pdim:1], col=gcol, main="thr=10")
image(out4$S[,pdim:1], col=gcol, main="automatic")
par(opar)