clear

%  Generate simple data set
dat = [1 -1; 2 -1; -5 3; 2 -1];

%  Check for zero mean:
mean(dat)

%  The mean is zero.  If we needed to remove the mean, we could do the
%  following:
dat = dat - (ones(4,1)*mean(dat));



%  Perform EOF analysis by SVD:
[u,s,v] = svd(dat)

%  Look at percent variance explained:
per = diag(s).^2/sum(diag(s).^2)



%  Redo this by calculating eigenvectors of the covariance matrix
c = dat'*dat/3;
[lds, lam] = eig(c)

%  Note that the 'eig' routine orders columns (EOF's) by increasing
%  eigenvalue (in contrast to the SVD routine above).  So, sort these
%  out:
[lam, ind] = sort(diag(lam));
ind = flipud(ind);
lam = lam(ind)
lds = lds(:,ind)



%  Compare the two methods:

%  Compare 'lds' and 'v'.
v
lds

%  Calculate PC's
pcs = dat*lds

%  Make sure we have consistency between lam, pcs, and s:
lam
var(pcs)
diag(s.^2/3)

per
lam/sum(lam)

%  Last, check out 'pcs' and 'u':
pcs
u*s