Nothing
# Original name : Donoho95
#' Covariance Estimation via Soft Thresholding
#'
#' Soft Thresholding method for covariance estimation takes off-diagonal elements \eqn{z} of sample covariance matrix and applies
#' \deqn{h_{\tau}(z) = \textrm{sgn}(z)(|z|-\tau)_{+}}
#' where \eqn{\textrm{sgn}(z)} is a sign of the value \eqn{z}, and \eqn{(x)_+ = \textrm{max}(x,0)}. If \code{thr} is rather a vector of regularization parameters, it applies
#' cross-validation scheme to select an optimal value.
#'
#' @param X an \eqn{(n\times p)} matrix where each row is an observation.
#' @param thr user-defined threshold value. If it is a vector of regularization values, it automatically selects one that minimizes cross validation risk.
#' @param nCV the number of repetitions for 2-fold random cross validations for each threshold value.
#' @param parallel a logical; \code{TRUE} to use half of available cores, \code{FALSE} to do every computation sequentially.
#'
#' @return a named list containing: \describe{
#' \item{S}{a \eqn{(p\times p)} covariance matrix estimate.}
#' \item{CV}{a dataframe containing vector of tested threshold values(\code{thr}) and corresponding cross validation scores(\code{CVscore}).}
#' }
#'
#' @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)
#'
#' @references
#' \insertRef{antoniadis_regularization_2001}{CovTools}
#'
#' \insertRef{donoho_wavelet_1995}{CovTools}
#'
#' @rdname CovEst.soft
#' @export
CovEst.soft <- function(X, thr=0.5, nCV=10, parallel=FALSE){
#-----------------------------------------------------
## PREPROCESSING
fname = "CovEst.soft"
pnameTHR = "'thr'"
pnamenCV = "'nCV"
pnamenthrs = "'nthrs'"
# 1. data matrix
checker1 = invisible_datamatrix(X, fname)
# 2. thr
if (length(as.vector(thr))==1){
checker3 = invisible_PosReal(thr, fname, pnameTHR)
CV = FALSE
} else { # vector threshold value case
nthrs = length(thr)
for (i in 1:nthrs){
checker3 = invisible_PosReal(thr[i], fname, pnameTHR)
}
CV = TRUE
}
# 4. nCV
if (CV==TRUE){
nCV = as.integer(nCV)
checker4 = invisible_PosIntMM(nCV, fname, pnamenCV, 1, nrow(X))
}
# 5. parallel
if (!parallel){
nCore = 1
} else {
nCore = max(round(detectCores()/2),1)
}
#-----------------------------------------------------
## MAIN COMPUTATION
if (CV==FALSE){
output = thr1.once(X,thr,func_soft)
} else {
output = thr1.multiple(X,nCV,nCore,func_soft_givenS,thr)
}
#-----------------------------------------------------
## RETURN OUTPUT
return(output)
}
# Auxiliary Functions for Soft Thresholding -------------------------------
#' @keywords internal
#' @noRd
func_soft <- function(X, thr){
S = cov(X)
diagS = diag(S)
output = array(0,dim(S))
idxPos = which(S>abs(thr))
idxNeg = which(S<(-abs(thr)))
output[idxPos] = S[idxPos]-thr
output[idxNeg] = S[idxNeg]+thr
diag(output) = diagS
return(output)
}
#' @keywords internal
#' @noRd
func_soft_givenS <- function(S, thr){
diagS = diag(S)
output = array(0,dim(S))
idxPos = which(S>abs(thr))
idxNeg = which(S<(-abs(thr)))
output[idxPos] = S[idxPos]-thr
output[idxNeg] = S[idxNeg]+thr
diag(output) = diagS
return(output)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.