% 3-20-2007
% Script file: irisbayes.m
% Training Quadratic Classifier from datairis.txt (75 patterns)
%
d=4; N=150; n1=25; n2=25; n3=25;  n=n1+n2+n3;
fin=fopen('datairis.txt');
fgetl(fin); fgetl(fin); fgetl(fin);   % skip three header lines
A=fscanf(fin,'%f',[d+1  N]);
X=A(1:d,:);
X1=X(:, 1:25);  X2=X(:,51:75);   X3=X(:,101:125);
Y1=X(:,25:50);  Y2=X(:,76:100);  Y3=X(:,126:150);
%
% Compute prior probabilities: p1,p2,p3\n
%         mean vectors u1\n,  u2\n,  u3\n
%         class-conditional covariance matrices S1, S2, S3 for each category
%         withn-class scatter matrix W
%
p1=1/3;  p2=1/3;  p3=1/3;
u1=mean(X1,2);  u2=mean(X2,2);  u3=mean(X3,2);  
S1=cov(X1');    S2=cov(X2');    S3=cov(X3');   % using an unbiased estimator
W=((n1-1)*S1+(n2-1)*S2+(n3-1)*S3)/3;           % within-class scatter matrix

fout=fopen('irisbayes.txt','w');
fprintf(fout,'probabilities, mean vectors, ');
fprintf(fout, 'within-class and toal scatter matrices\n');
fprintf(fout,'%8.4f\t%8.4f\t%8.4f\n',p1,p2,p3);
for i=1:d-1
  fprintf(fout,'%8.4f\t',u1(i));
end
fprintf(fout,'%8.4f\n',u1(d));
for i=1:d-1
  fprintf(fout,'%8.4f\t',u2(i));
end
fprintf(fout,'%8.4f\n',u2(d));
for i=1:d-1
  fprintf(fout,'%8.4f\t',u3(i));
end
fprintf(fout,'%8.4f\n',u3(d));
%
% Print class-conditional covariance matrices: S1, S2, S3
%
for i=1:d
  for j=1:d
    fprintf(fout,'%8.4f',S1(i,j));
  end
  fprintf(fout,'\n');
end

% Line 49: End-of-Page
for i=1:d
  for j=1:d
    fprintf(fout,'%8.4f',S2(i,j));
  end
  fprintf(fout,'\n');
end
for i=1:d
  for j=1:d
    fprintf(fout,'%8.4f',S3(i,j));
  end
  fprintf(fout,'\n');
end
%
% Print within-class scatter matrix W
%
for i=1:d
  for j=1:d
    fprintf(fout,'%8.4f',W(i,j));
  end
  fprintf(fout,'\n');
end
%
% For Testing 
%
err1=0;  err2=0;  err3=0;
  r12=log(p1*p1*det(S2)/p2*p2*det(S1));
  r13=log(p1*p1*det(S3)/p3*p3*det(S1));
  r23=log(p2*p2*det(S3)/p3*p3*det(S2));
  if (i<=n1)
    x=Y1(:,i);
    R1=mhdist(x,u1,S1); R2=mhdist(x,u2,S2); R3=mhdist(x,u3,S3);
    if ( (R1-R2<r12) || (R1-R3)<r13 ) 
      err1=err1+1;
    end
  end
  if (i<=n2)
    x=Y2(:,i);
    R1=mhdist(x,u1,S1); R2=mhdist(x,u2,S2); R3=mhdist(x,u3,S3);
    if ( (R2-R1<-r12) || (R2-R3)<r23 ) 
      err2=err2+1;
    end
  end
  if (i<=n3)
    x=Y3(:,i);
    R1=mhdist(x,u1,S1); R2=mhdist(x,u2,S2); R3=mhdist(x,u3,S3);
    if ( (R3-R1<-r13) || (R3-R2)<-r23 ) 
      err3=err3+1;
    end
  end
  
[err1,err2,err3,n]
fprintf(fout,'----------------------------------\n');
fprintf(fout,'(err1,err2,err3,n)=(%2d,%2d,%2d,%3d)\n',err1,err2,err3,n);
fclose(fin);
fclose(fout);