rafisher2cda.m

Fisher线性判别分类器应用于IRIS数据集的例子

if pp == 1;  %groups'observed probability (unknown)
   pr = n/sum(n);  
else pp == 2;  %groups'expected probability (known)
   pr = input('Give me the vector of the known prior probabilities:');
   if sum(pr)~= 1
      error('The sum of the known prior probabilities must to be one. Check it or adjust it.');
      return;
   end;
end;

[n,p] = size(Y);  %total number of data (n) and number of variables (p)

%Deviations of data from the grand mean (total).
dT = [];
for j = 1:p
   eval(['dT  = [dT,(Y(:,j) - mean(Y(:,j)))];']);
end;

%Deviations of samples from their own mean (within).
dW = [];
for k = 1:g
    q = find(G==k);
    mG = mean(Y(min(q):max(q),:));
    w = [];
    for l=1:p
       w = [w,Y(min(q):max(q),l)-mG(l)];
    end;
    dW = [dW;w];
end;

%Re-scaling procedure.
S = std(dW);
if any(S < n*max(S)*eps)
   error(sprintf(['Column %d in feature matrix X is constant within' ...
         ' groups.'], min(find(S < n*max(S)*eps))))
end;
S = diag(S);

pe = input('Do you want an unbiased (1) or maximum-likelihood parameter estimation? (2):');
if pe == 1;
   es = 0;
else (pe == 2);
   es = 1;
end;

[u s v] = svd(dW*S/sqrt(n-g*(1-es)),0);
r = sum(diag(s) > n*s(1)*eps);
if (r < p)
   warning(sprintf(['Nullity of within-groups covariance matrix is' ...
         ' %d.'], p-r))
   v = v(:,1:r);
   s = s(1:r,1:r);
end;
S = S*inv(triu(qr(s*v')));

Ms = diag(sqrt(n*pr/(g-es)))*(M-repmat(pr*M,g,1))*S;
[u s v] = svd(Ms,0);
r = sum(diag(s) > n*eps*s(1));
cdf = S*v(:, 1:r);  %canonical discriminant functions (cdf)
ncdf = min(p,g-1);  %number of cdf to retain
rcdfr = cdf(:,1:ncdf);   %retained canonical discriminant functions
                      %variates of the canonical discriminant analysis
Z = Y;

tmpsc = Z*rcdfr; %scores of the canonical discriminant analysis
ct = mean(tmpsc);  %constants of the canonical discriminant analysis
Constant = ct;
Variates = rcdfr;

dB = dT - dW;  %deviations between samples

W = dW'*dW;  %within sum of squares
B = dB'*dB;  %between sum of squares

L = eig(inv(W)*B);
L = max([L'; zeros(size(L'))]);  %ignore negative eigenvalues
L = fliplr(sort((L)));
L = L(1,1:ncdf);  %retain subset
Eigenvalue = L;
R2 = L./(1+L);
R = sqrt(R2);  %canonical coorrelation coefficients

pvar = 100*L/sum(L);  %percentage to total variance
Percentage = pvar;
cpvar = cumsum(pvar);  %cumulative percentage of variance
CumulativePercentage = cpvar;

disp(' ')
disp('Canonical Discriminant Functions.')
fprintf('--------------------------------------------------------------------------------\');
Constant
Variates
Eigenvalue
Percentage
CumulativePercentage
fprintf('--------------------------------------------------------------------------------\n');
fprintf('Functions = columns. On variates, Variate1 = first row and so forth to %.i\n', p);

%Bartlett's approximate chi-squared statistic for testing
%the canonical correlation coefficients
d = ncdf;
k = 0:(d-1);
df = (p-k).*(g-k-1);  %Chi-square statistic degrees of freedom
LL = 1./(1+L);
LW = fliplr(cumprod(fliplr(LL)));  %Wilk's lambda vector
v = n-1;  %total degrees of freedom
X2 = -(v - .5*((p-k)+(g-k-1)+1)).* log(LW);  %Chi-square statistic
P = 1 - chi2cdf(X2,df);  %p-value associated to the Chi-square statistic

disp(' ')
disp('Chi-square Tests with Successive Roots Removed.')
fprintf('-------------------------------------------------------------------------\n');
disp('Removed    Eigenvalue   CanCor      LW       Chi-sqr.     df       P')
fprintf('-------------------------------------------------------------------------\n');
fprintf('%4.i%15.4f%11.4f%10.4f%13.4f%7.i%10.4f\n',[k',L',R',LW',X2',df',P'].');
fprintf('-------------------------------------------------------------------------\n');
fprintf('With a given significance of: %.2f\n', alpha);
disp('If P-value >= alpha, it is not significative. Else, it results significative.')

if c == 1;
   
   disp(' ');
   m = size(tmpsc,2);
   CO = [];
   for i = 1:m
      for s = i+1:m
         if s ~= i
            co = [i s];
            CO = [CO;co];
         end;
      end;
   end;
   
   plots = CO;
   ct = (m*(m - 1))/2;
   fprintf('The pair-wise plots you can get are: %.i\n', ct);
   fprintf('--------------\');
   plots
   fprintf('--------------\n');
   
   d='y';
   while d=='y';
      sd=input('Give me the pair of factors to plot [a b]: ');
      sc = [tmpsc(:,sd(1)) tmpsc(:,sd(2))];
      
      figure;
      hold on;
      
      sc=[G sc];
      k = size(sc,2)-1;
      
      lg = [];
      for k = 1:g
         gr = find(G==k); % get row indices for each group
         plot(sc(gr,2),sc(gr,3),dcsymb0(k),'MarkerFaceColor',dcrgb0(k),'MarkerEdgeColor',dcrgb0(k));
         lg = [lg,['''Group ' num2str(k) ''',']];
      end;
      lg(end)=' ';
      eval(['legend(' lg ')']);
      
      [r c] = size(sc);
      x = sc(:,2:c);
      uG  = unique(G);  %unique groups
      centroid(g,size(x,2)) = 0;  %preallocate results array
      
      for k = 1:g
         subsr = find(G==uG(k));  %get indices of rows to extract
         subsx = x(subsr,:);  %extraction of group separately
         centroid(k,:) = mean(subsx,1);  %compute group centroid
      end;
      
      for k = 1:g
         h = text(centroid(k,1),centroid(k,2),num2str(k));
         set(h,'HorizontalAlignment','center','FontWeight','bold','FontSize',10,'Color',[0 0 0]);
      end;
      
      %Plotting of ellipse
      title('Canonical  Discriminant  Analysis');
      xlabel(['Canonical  Function ' num2str(sd(1)) ;]);
      ylabel(['Canonical  Function ' num2str(sd(2)) ;]);
      indice = sc(:,1);
      Dx = [];Dy = [];aa=[];
      for k = 1:g
         Xe = find(indice==k);
         eval(['x' num2str(k) '= sc(Xe,2);']);
         eval(['y' num2str(k) '= sc(Xe,3);']);
         eval(['n' num2str(k) '= length(x' num2str(k) ') ;']); 
         %Group means vector for the canonical scores.
         eval(['mx' num2str(k) '= mean(x' num2str(k) ') ;']);
         eval(['my' num2str(k) '= mean(y' num2str(k) ') ;']);
         eval(['X' num2str(k) '= [x' num2str(k) ',y' num2str(k) '] ;']);
         %Group covariances for the canonical scores.
         eval(['S' num2str(k) '= cov(X' num2str(k) ') ;']);
         %Group eigenstructure of the covariances.
         eval(['[V' num2str(k) ', L' num2str(k) '] = eig(S' num2str(k) ') ;']);
         %Group eigenvalues.
         eval(['L' num2str(k) '= diag(L' num2str(k) ') ;']);
         %Group angle rotation in radians for a 2-D ellipse.
         eval(['a' num2str(k) '= acos(V' num2str(k) '(1,1)) ;']);
         %Critical value for confidence bounds for the group ellipses.
         eval(['c' num2str(k) '= 2*finv(1-alpha,2,n' num2str(k) '-2) ;']);
         %Group length for the major axis.
         eval(['A' num2str(k) '= sqrt(c' num2str(k) '*max(L' num2str(k) ')) ;']);
         %Group length for the minor axis.
         eval(['B' num2str(k) '= sqrt(c' num2str(k) '*min(L' num2str(k) ')) ;']);
         %Parametrize of the group ellipses by angles around unit circle.
         eval(['t' num2str(k) '= linspace(0,2*pi,n' num2str(k) ') ;']);
         eval(['Xdata' num2str(k) '= A' num2str(k) '*cos(t' num2str(k) ') ;']);
         eval(['Ydata' num2str(k) '= B' num2str(k) '*sin(t' num2str(k) ') ;']);
         eval(['Ydata' num2str(k) '(end)= 0 ;']);
         %Form of the group contour points.
         eval(['x' num2str(k) '= cos(a' num2str(k) ')*Xdata' num2str(k) '-sin(a' num2str(k) ')*Ydata' num2str(k) '+mx' num2str(k) ' ;']);
         eval(['y' num2str(k) '= sin(a' num2str(k) ')*Xdata' num2str(k) '+cos(a' num2str(k) ')*Ydata' num2str(k) '+my' num2str(k) ' ;']);
         eval(['a = a' num2str(k) ';']);
         eval(['dx = x' num2str(k) ';']);
         eval(['dy = y' num2str(k) ';']);
         eval(['mx = mx' num2str(k) ';']);
         eval(['my = my' num2str(k) ';']);
         eval(['ss= S' num2str(k) ';']);
         Dx = [Dx,dx];Dy = [Dy,dy];aa=[aa;a];
      end;
      ed = [G Dx' Dy'];
      
      for k = 1:g
         gr = find(G==k); % get row indices for each group
         plot(ed(gr,2),ed(gr,3),dcsymb1(k),'MarkerFaceColor',dcrgb1(k),'MarkerEdgeColor',dcrgb1(k));  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      end;
      
      d=input('Do you need another plot? (y/n): ','s');
   end; 
else
end;