function a=Fisher_discriminants(X1,X2,X3)
%function Fisher_discriminants(X1,X2,X3)
% tre popolazioni
% Grafico 2D di discriminanza
% a = [a1, a2]  primo e secondo discriminante 


[n1,p1]=size(X1); [n2,p1]=size(X2); [n3,p1]=size(X3);

g=3;
n = n1+n2+n3;
m1 = mean(X1); m2 = mean(X2); m3 = mean(X3);
S1 = cov(X1); S2 = cov(X2); S3 = cov(X3);

W = (n1-1)*S1+ (n2-1)*S2+ (n3-1)*S3;
m = mean([X1;X2;X3]);           %m = mean([m1;m2;m3]);

B1 = (m1-m)'*(m1-m); B2 = (m2-m)'*(m2-m); B3 = (m3-m)'*(m3-m);
B = n1*B1+ n2*B2+ n3*B3;


[V,E]=eig(W\B);

[e,ie]=sort(diag(E),'descend'); V=V(:,ie);
e
Spooled=W/(n1+n2+n3-g);
a(:,1:2)=V(:,1:2);
a=a/diag(sqrt(diag( a'*Spooled*a)))


plot(X1*a(:,1),X1*a(:,2),'ro')
hold on
plot(X2*a(:,1),X2*a(:,2),'*')
plot(X3*a(:,1),X3*a(:,2),'+k')
legend('Fam. 1','Fam. 2','Fam. 3')
xlabel('y_1')
ylabel('y_2')
plot(mean(X1*a(:,1)),mean(X1*a(:,2)),'>r')
plot(mean(X2*a(:,1)),mean(X2*a(:,2)),'p')
plot(mean(X3*a(:,1)),mean(X3*a(:,2)),'kd')
hold off


return
