%	tsvd		- test properties of SVD
%
%	This program
%	- generates n 2-D points
%	- find the principal components via eigen decomposition (PCA)
%	- plot the principal components
%	- find the principal components via Singular Value Decomposition	
%	- plot the principal components
%
%	vary n from 3 to 100 see how well the 2 methods agree
%

clf;

n	= 3;		%	# data points
dev	= .9;		%	deviation


x	= 2*rand(2,n)-1;
x(2,:)	= 2*x(1,:)+(2*rand(1,n)-1)*dev;


%plotmd(x,'r.');
plot(x(1,:),x(2,:),'k.');

%[pc_data,pvar,paxis]	= pca(x);
[v,d]	= eig(cov(x'))
lambda	= diag(d);

% reorder in descending order of eigenvalues
if lambda(2)>lambda(1)
	v	= [v(:,2) v(:,1)];
	d	= diag([lambda(2) lambda(1)]);
end

%plotaxis(paxis,'b-','-plotover 1');

component=1;
hold on;
plot([v(1,component) 0],[v(2,component) 0],'r-');
component=2;
hold on;
plot([v(1,component) 0],[v(2,component) 0],'g--');

[u,s,v]	= svd(x);
lambda	= diag(s);

% reorder in descending order of singular values
if lambda(2)>lambda(1)
	u	= [u(:,2) u(:,1)];
	d	= diag([lambda(2) lambda(1)]);
end

component=1;
hold on;
plot([u(1,component) 0],[u(2,component) 0],'b-.');
component=2;
hold on;
plot([u(1,component) 0],[u(2,component) 0],'m:');

legend('data','PCA1','PCA2','SVD1','SVD2');
axis ([-1 1 -1 1]);
axis square