# sicb2013 / regularization_functions.r

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65``` ```#' regularization_by_thresholding #' #' Regularizes a correlation matrix by thresholding #' #' @param cor_matrix A correlation matrix #' #' @param n The number of independent observations that were used to construct #' the #' #' @return The regularized correlation matrix #' #' @references #' Bickel, P. J. & Levina, E. Covariance regularization by thresholding. #' Ann. Statist. 36, 2577–2604 (2008). http://dx.doi.org/10.1214/08-AOS600 #' regularization_by_thresholding <- function ( cor_matrix, n ) { # Regularization by thresholding # Bickel, P. J. & Levina, E. Covariance regularization by thresholding. # Ann. Statist. 36, 2577–2604 (2008). http://dx.doi.org/10.1214/08-AOS600 p <- ncol( cor_matrix ) regularized <- cor_matrix * ( abs(cor_matrix) > sqrt(log(p)/n) ) return( regularized ) } #' regularization_by_convex_minimization #' #' Regularizes a correlation matrix by convex minimization #' #' @param cor_matrix A correlation matrix #' #' @param n The number of independent observations that were used to construct #' the #' #' @return The regularized correlation matrix #' #' @references #' Luo, X. High Dimensional Low Rank and Sparse Covariance Matrix Estimation via #' Convex Minimization. arXiv.org (2011). http://arxiv.org/abs/1111.1133 #' regularization_by_convex_minimization <- function ( cor_matrix, n ) { p <- ncol( cor_matrix ) re.lorec <- lorec( cor_matrix, diag(1, 100), diag(1, 100), sqrt(p/n), sqrt(log(p)/n) ) re.lorec.eig <- eigen( re.lorec\$L ) # threholding both sel <- re.lorec.eig\$values > sqrt( p/n ) if ( sum(sel) > 0 ) { V <- re.lorec.eig\$vectors V <- V * ( abs(V)>sqrt(1/p) ) V <- V[ ,1:sum(sel) ] Lhat <- V%*%diag( re.lorec.eig\$values[1:sum(sel)] )%*%t(V) } else { Lhat <- matrix(0, p, p) } Shat <- re.lorec\$S Shat <- Shat*(abs(Shat)>sqrt(1/p)) regularized <- Lhat + Shat return( regularized ) } ```