%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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)