%Funcao Discriminante Canonica de Fisher
%Amostras
clc;
Dados=[1 -4 10
       1  0  6
       1 -2  2
       2  0 12
       2  4  8
       2  2  4
       3  2 -4
       3  0  0
       3 -2 -8];
X=Dados(:,2:end);

%medias entre grupos (ou médias dos grupos)
uw=[mean(X(1:3,:));mean(X(4:6,:));mean(X(7:end,:))];%media do grupo 1, média do grupo 2, média do grupo 3
%p: numero de medidas da amostra, ou variaveis ou campos
%n: tamanho de uma amostra de p dimensoes
[n,p]=size(X);
%media geral
u=mean(uw);
%numero de classes ou grupos
g=3;

%Imprime os dados nao transformados
subplot(1,2,1);
hold on;
xlim([min(X(:,1))-1 max(X(:,1))+1]);
ylim([min(X(:,2))-1 max(X(:,2))+1]);
for i=1:g,%para cada grupo
   
      k=(i-1)*3+1:(i-1)*3+3;  
      if i==1
          %vermelho
          plot(X(k,1),X(k,2),'ro','LineWidth',2);
          
      elseif i==2 
          %verde
          plot(X(k,1),X(k,2),'go','LineWidth',2);
         
      else
          %azul
          plot(X(k,1),X(k,2),'bo','LineWidth',2);
         
      end  
end
grid on;
title('Dados brutos');   
%soma dos quadrados dentro (Within) dos grupos (soma dos quadrados das médias dos grupos)
%pega cada grupo e calcula sua variacao (covariancia sem dividir por n), depois soma
W=zeros(p,p);
for i=1:g,
  %k=1..3,4..6,7..9  
  k=(i-1)*3+1:(i-1)*3+3;
  %(x-uw)*(x-uw), 2 colunas
  %             3x2'*3x2 = 2x3*3x2 => 2x2
  W=W+(X(k,:)-ones(3,1)*uw(i,:)).'*(X(k,:)-ones(3,1)*uw(i,:));
end;  
disp('Soma dos quadrados dentro dos grupos');
W
%soma dos quadrados entre (Between)
%soma da variacao entre a media de cada grupo e a media geral
B=3*(uw - ones(3,1)*u).'*(uw-ones(3,1)*u);
disp('Soma dos quadrados entre grupos (media(i) - mediageral)');
B
%(W^-1)*B = W\B
W_1B=W\B;
disp('Calculando W^(-1)*B:');
W_1B
%calcula autovetores V e autovalores D (de diagonal)
[V,D]=eig(W_1B);
%ordena tal que d1>d2>d3
[sD,index]=sort(diag(D),'descend');
%ordena os autovetores de acordo com seu autovalor correspondente
e=V(:,index);
%reescala os autovetores para ai'Wai=1, 
%se ei'Wei=k, entao ai=ei/(k^0.5) e faz ai'Wai que sera igual a 1
a=zeros(p,p);
a_=a;
for i=1:p,
    a(:,i)=e(:,i);
    k=e(:,i)'*W*e(:,i);
    if k~=1
        a(:,i)=e(:,i)/(k^0.5);
    end  
    a_(:,i)=a(:,i)/norm(a(:,1));
end;    
disp('Autovetores escalonados: ');
a
%calculando as funcoes discriminantes no vetor de medias amostrais
fprintf('\nFuncao discriminante das medias amostrais:\n');
uY=uw*a
%calculando
fprintf('\nEncontrando a classe para x:\n');
x=[4 6]
fprintf('\nFuncao discriminante para x:\n');
y=x*a
fprintf('\nDistancias quadradas de y a uY:\n');
d2=(ones(3,1)*y(1)-uY(:,1)).^2+(ones(3,1)*y(2)-uY(:,2)).^2
fprintf('\nMenor distancia e classe:\n');
[menorDistancia,classe]=min(d2)
Y=X*a;
subplot(1,2,2);
xlim([min(Y(:,1))-0.1 max(Y(:,1))+0.1]);
ylim([min(Y(:,2))-0.1 max(Y(:,2))+0.1]);
hold on;
for i=1:g,
      k=(i-1)*3+1:(i-1)*3+3;  
      if i==1
          %vermelho
          plot(Y(k,1),Y(k,2),'ro','LineWidth',2);
          plot(uY(i,1),uY(i,2),'rs','LineWidth',3);
          text(uY(i,1),uY(i,2), ' media C1');
      elseif i==2 
          %verde
          plot(Y(k,1),Y(k,2),'go','LineWidth',2);
          plot(uY(i,1),uY(i,2),'gs','LineWidth',3);
          text(uY(i,1),uY(i,2), ' media C2');
      else
          %azul
          plot(Y(k,1),Y(k,2),'bo','LineWidth',2);
          plot(uY(i,1),uY(i,2),'bs','LineWidth',3);
          text(uY(i,1),uY(i,2), ' media C3');
      end  
end
title('Dados transformados'); 
plot(y(1),y(2),'kp','LineWidth',3);  
text(y(1),y(2), '  ponto a classificar');
grid on;