// Programa: mahalanobis.sce
// Autor: Hemerson Pistori (http://www.ec.ucdb.br/~pistori) 
// Última Alteração: 23/09/2005
//
// Descrição:
// 
// Script para scilab que ilustra a utilização da distância
// de mahalanobis em um espaço bidimensional.
//
// Instruções para utilização:
// 1. Entre no scilab e verifique em que diretório você se encontra:
// pwd();
//
// 2. Mude, se necessário, para o diretório que contém o script
// chdir('diretorio_contendo_o_script'); 
//
// 3. Carregue o script
// exec('mahalanobis.sce');
//
// Leia instruções abaixo sobre como utilizar cada uma
// das funções. Para um exemplo, digite, após carregar o
// script: exemplo().
//



//////////////////////////////////////////////////////////////////////////
// Cálcula a distância Euclidiana entre um vetor e uma distribuição
// estatística gaussiana.
//
// Parâmetros:
//    vetor = Um vetor
//    media = Vetor média da distribuição gaussiana
//    sigma = Matriz de variância e covariância
//
// Saída:
//    f = distância Euclidiana
//
// Exemplo de utilização
//
// --> distE([10 12]',[20 20]',[12 0:0 15]');
//
function f=distE(vetor,media,sigma)
  f=norm(vetor-media);
endfunction


//////////////////////////////////////////////////////////////////////////
// Cálcula a distância de Mahalanobis entre um vetor e uma distribuição
// estatística Gaussiana.
//
// Parâmetros:
//    vetor = Um vetor
//    media = Vetor média da distribuição gaussiana
//    sigma = Matriz de variância e sigmaiância
//
// Saída:
//    f = distância de Mahalanobis
//
// Exemplo de utilização
//
// --> distM([10 12]',[20 20]',[12 0:0 15]');
//
function f=distM(vetor,media,sigma)
  f=sqrt((vetor-media)'*inv(sigma)*(vetor-media));
endfunction



//////////////////////////////////////////////////////////////////////////
// Desenha duas distribuições gaussianas e mostra a distância de
// de Mahalanobis de um ponto em relação as Gaussianas.
//
// Exemplo de utilização
//
// --> exemplo()
//
function exemplo()

   media1=[10 10]';
   sigma1=[6 0; 0 25];
   media2=[20 20]';
   sigma2=[5 10;10 25];
   intervalo=[0:30];
   ponto=[15 20]';


   z1 = zgauss(intervalo,intervalo,media1,sigma1);
   z2 = zgauss(intervalo,intervalo,media2,sigma2);
   z = max(z1,z2);  // Agrupando as duas distribuições Gaussianas
   z(ponto(1),ponto(2))=180;  // Criando uma saliência para destacar o lugar onde está o ponto

   [xx,yy,zz] = genfac3d(intervalo,intervalo,z); // Para mostrar com cores diferentes 
   plot3d1(xx,yy,zz,alpha=3,theta=50);

   printf("Média distribuição 1: \n");
   disp(media1);
   printf("Média distribuição 2: \n");
   disp(media2);
   printf("Ponto:\n");
   disp(ponto);
   printf("Distância Euclidiana para Distribuição 1: \n");
   disp(distE(ponto,media1,sigma1));
   printf("Distância Euclidiana para Distribuição 2: \n");
   disp(distE(ponto,media2,sigma2));
   printf("Distância de Mahalanobis para Distribuição 1: \n");
   disp(distM(ponto,media1,sigma1));
   printf("Distância de Mahalanobis para Distribuição 2: \n");
   disp(distM(ponto,media2,sigma2));
   printf("\n");
 

endfunction



//////////////////////////////////////////////////////////////////////
// Função auxiliar para apresentar graficamente os modelos inferidos.
// Retorna a matriz pronta para ser utilizada por plo3d
// (mais eficiente pois não recalcula, a cada iteração, as "constantes"
//  da função gaussiana)
// x - pontos no eixo X
// y - pontos no eixo Y
// x e y devem possui o mesmo tamanho
// mu - vetor média
// sigma - matriz de variância e sigmaiância
function z=zgauss(x,y,mu,sigma)

  c1 = (1/(  (2*%pi)^(length(mu)/2)    *sqrt(det(sigma))));
  c1 = c1*10000; // Aumenta escala para melhorar visualização plot3d
  isigma = inv(sigma);
  n = length(x);
  z = zeros(n,n);
  for i = 1:n
    for j = 1:n
       v = [x(i) y(j)]';
       z(i,j) = c1 * exp( -0.5 * (v-mu)'*isigma*(v-mu))
    end;
  end;


endfunction