function ex2(n, m)
% Illustration of using principal component analysis
% to reduce the data dimensionality from 3 to m

% bounds of the range of the input x
lx = -5;
ux = +5;

u = ones(3, 1);
if m > 0 & m <= 3
    u(1:3-m) = 0;
end

if n > 0
    YnoiseLevel = 0.4;
    ZnoiseLevel = 0.5;

    % define the target function
    f = inline('3*sin(x)', 'x');
    X = linspace(lx, ux, 100);
    Y = f(X);
    Z = zeros(1, 100);

    % plot the target function
    plot3(X, Y, Z, '-b');
    grid on
    hold on
    xlabel('x')
    ylabel('y')
    zlabel('z')
    
    % data
    x = linspace(lx, ux, n);

    % row vector of target values
    y = f(x) + YnoiseLevel * (rand(size(x)) - 0.5);

    % row of noise
    z = ZnoiseLevel * (rand(size(x)) - 0.5);
    
    % plot the data
    plot3(x, y, z, '+r');
    
    xlim([lx ux]); 
    ylim([lx ux]);   
    zlim([-ZnoiseLevel ZnoiseLevel])   
    pause;

    % compute covariance matrix and its eigen-decomposition
    C = cov([x; y; z]');
    [V,D] = eig(C);    
	
    % plot 1st principal component
    plotv(zeros(3, 1), sqrt(D(3,3))*V(:,3), 'x´');
	
    % 2nd principal component
    plotv(zeros(3, 1), sqrt(D(2,2))*V(:,2), 'y´');

    % 3rd principal component
    plotv(zeros(3, 1), sqrt(D(1,1))*V(:,1), 'z´');

    % plot new data
    Q(1,1) = V(3,2:3)*(V(1:2,2:3)\[lx; lx]);
    Q(1,2) = V(3,2:3)*(V(1:2,2:3)\[ux; lx]);
    Q(2,1) = V(3,2:3)*(V(1:2,2:3)\[lx; ux]);
    Q(2,2) = V(3,2:3)*(V(1:2,2:3)\[ux; ux]);
    surf([lx ux], [lx ux], Q, 'FaceAlpha', 0.2);
    
    xlim([lx ux]); 
    ylim([lx ux]);   
    zlim([-ZnoiseLevel ZnoiseLevel])   
    pause;

    % map data to new coordinates
    Q = V*diag(u)*V'*[X; Y; Z];
    plot3(Q(1,:), Q(2,:), Q(3,:), '--g');
    Q = V*diag(u)*V'*[x; y; z];
    plot3(Q(1,:), Q(2,:), Q(3,:), 'og');

    xlim([lx ux]); 
    ylim([lx ux]);   
    zlim([-ZnoiseLevel ZnoiseLevel])   
    hold off
end


function plotv(p0, p1, str)
% plot 3d vector

x0 = p0(1);
y0 = p0(2);
z0 = p0(3);
x1 = p1(1);
y1 = p1(2);
z1 = p1(3);
plot3([x0;x1], [y0;y1], [z0;z1], '-k');   % Draw a line between p0 and p1

p = p1 - p0;
alpha = 0.1;  % Size of arrow head relative to the length of the vector
beta = 0.1;  % Width of the base of the arrow head relative to the length

hu = [x1-alpha*(p(1)+beta*(p(2)+eps)); x1; x1-alpha*(p(1)-beta*(p(2)+eps))];
hv = [y1-alpha*(p(2)-beta*(p(1)+eps)); y1; y1-alpha*(p(2)+beta*(p(1)+eps))];
hw = [z1-alpha*p(3); z1; z1-alpha*p(3)];

plot3(hu(:), hv(:), hw(:), '-k')  % Plot arrow head

if nargin > 2
    % plot str
    text(x1*1.04, y1*1.04, z1*1.04, str);
end
