%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % prerequites: % % % % the m/data files need to be put in a path where matlab/octave % % can find them, e.g., the directory from which matlab is called. % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % shrinkage covariance estimator % compared with usual unbiased sample estimator %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Example script (to be run in MATLAB or OCTAVE) % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % small data data set % with 10 variables (columns) % and 6 repetitions (rows) load smalldata.txt % tab-delimited data X=smalldata; size(X) % 6 x 10 % estimate 10x10 covariance matrix s1 = cov(X); % usual unbiased sample estimate [s2, lamcor, lamvar] = covshrinkKPM(X, 1); % shrinkage estimate % (fast (!) function by Kevin P. Murphy) lamvar % shrinkage intensity variances: 0.60155 lamcor % shrinkage factor correlations: 0.73152 % compare ranks and conditions rank (s1) % 5 cond (s1) % 3.7370e+17 (Inf) rank (s2) % 10 cond (s2) % 4.5130 % compare positive definiteness all ( eig(s1) > 0) % not positive definite all ( eig(s2) > 0) % positive definite % which estimate can be inverted? inv(s1); % this cannot: matrix singular to machine precision inv(s2); % no problem at all (s2 has full rank and is positive definite) %%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%% % large data data set % with 100 variables (columns) % and 20 repetitions (rows) load largedata.txt % tab-delimited data X=largedata; size(X) % 20 x 100 % estimate 100x100 covariance matrix s1 = cov(X); % usual unbiased sample estimate [s2, lamcor, lamvar] = covshrinkKPM(X, 1); % shrinkage estimate % (function by Kevin P. Murphy) lamvar % shrinkage intensity variances: 0.77718 lamcor % shrinkage intensity correlations: 0.88511 % compare ranks and conditions rank (s1) % 19 cond (s1) % 1.3859e+18 (Inf) rank (s2) % 100 cond (s2) % 2.8671 % compare positive definiteness all ( eig(s1) > 0) % not positive definite all ( eig(s2) > 0) % positive definite % which estimate can be inverted? inv(s1); % this cannot: matrix singular to machine precision inv(s2); % no problem at all (s2 has full rank and is positive definite)