function pca_eval()

close all;

Ntrain = 100;
Ndim = 2;
K = 1;
epsilon = 0.1;

x = 2*(rand(Ndim, Ntrain)-0.5);
x = 0.3*x + repmat(linspace(-0.7,0.7, Ntrain), Ndim, 1);
plot(x(1,:), x(2,:), '+'); axis([-1 1 -1 1]); grid; title('Raw data'); xlabel('x_1'); ylabel('x_2');

sigma = x * x' / size(x, 2)
[U,S,V] = svd(sigma)
xrot = U' * x;
figure;
plot(xrot(1,:), xrot(2,:), '+'); axis([-1 1 -1 1]); grid; title('Rotated data'); xlabel('x_{rot,1}'); ylabel('x_{rot,2}');

xpca = diag(sqrt(1./(diag(S) + epsilon))) * xrot;
figure;
plot(xpca(1,:), xpca(2,:), '+'); axis([-1 1 -1 1]); grid; title('Whitened data dim=2'); xlabel('x_{pca,1}'); ylabel('x_{pca,2}');


xtilde = zeros(size(xrot));
xtilde(1:K, :) = xrot(1:K, :);

lambda = diag(S);
retained_variance = sum(lambda(1:K))./sum(lambda)

figure;
plot(xtilde(1,:), xtilde(2,:), '+'); axis([-1 1 -1 1]); grid; title('Rotated data reduced to dim = 1'); xlabel('x_{tilde,1}'); ylabel('x_{tilde,2}');

size (U)
size(xrot)
xhat = U(:, 1:K) * xrot(1:K, :);

figure;
plot(xhat(1,:), xhat(2,:), '+'); axis([-1 1 -1 1]); grid; title('Data reduced to dim = 1'); xlabel('x_{hat,1}'); ylabel('x_{hat,2}');

xpca1 = diag(sqrt(1./(diag(S)))) * xhat;
figure;
plot(xpca1(1,:), xpca1(2,:), '+'); axis([-1 1 -1 1]); grid; title('Whitened data dim=1'); xlabel('x_{pca,1}'); ylabel('x_{pca,2}');