% Primary Component Analysis
clear
px=linspace(-2,2,50);
py=linspace(-2,2,50);
px=px+randn(1, 50)*0.2+1; py=py+randn(1,50)*0.5+1;
plot(px,py,'bx');
Ap=cov(px,py);
hold on;
xlabel('feature s1');
ylabel('feature s2');
G_px=mean(px); G_py=mean(py);
plot(G_px,G_py,'ro');
axis([-5,5,-5,5]);
grid on;
ramda1=(Ap(1,1)+Ap(2,2)+sqrt((Ap(1,1)+Ap(1,2))^2-4*(Ap(1,1)*Ap(2,2)-Ap(1,2)^2)))/2;
ramda2=(Ap(1,1)+Ap(2,2)-sqrt((Ap(1,1)+Ap(1,2))^2-4*(Ap(1,1)*Ap(2,2)-Ap(1,2)^2)))/2;
if abs(ramda1)>abs(ramda2)
    ramda=ramda1;
else
    ramda=ramda2;
end;
a=Ap(1,2)/sqrt(Ap(1,2)^2+(ramda-Ap(1,1))^2);
b=(ramda-Ap(1,2))/sqrt(Ap(1,2)^2+(ramda-Ap(1,1))^2);
for n=1:50
    D(n)=a*(px(n)-G_px)+b*(py(n)-G_py);
end;

for x=-5:0.1:5
     y=a/b*x;
     plot(x,y,'k.');
end

hold off;
