function [RAFisher2cda] = RAFisher2cda(X,pp,c,alpha)
% RAFISHER2CDA Canonical Discriminant Analysis.
%
%[As well as the previous RAFisher1, this file was named RAFisher2 in honor to
%Sir Ronald Aylmer Fisher (17 Feb.,1890-29 Jul.,1962), one of the most insightful
%and influential statisticians of all times. His pioneering work on discriminant
%functions was made in 1936 (a pdf file of this paper is attached in this zip),
%presenting an analysis on data collected by Edgar Anderson on three species of
%iris flowers (a classical set of multivariate data and also attached in this zip).
%Also, for all those curious people we attach the jpg-images of the three Iris plant
%species.]
%
%-Again, we encourage to all the users to initiate a chain of petitions to make
%possible to carry to the big or small screen the life of this notable genius of our
%time. A possible title of the film would be 'Fisher's geometric mind'.-His work shaped
%the world of statistics, he made it as we know it today.-
%
%While RAFisher1 is a procedure that produces very different functions for classification
%that are also called linear discriminant analysis, RAFisher2cda is a dimension-reduction
%technique related to principal component analysis and canonical correlation called
%canonical discriminant analysis. It derives the canonical coefficients parallels that
%of one-way MANOVA and it finds linear combinations of the quantitative variables that
%provide maximal separation between the classes or groups in much the same way that
%principal components summarize total variation. The output produced are the canonical
%coefficients and the scored canonical variables. The canonical coefficients are rotated.
%Also, it proceeds with a Bartlett's approximate chi-squared statistic for testing the
%canonical correlation coefficients.
%
%In summary, the canonical discriminant analysis:
% - Transform the variables so that the pooled within-group covariance matrix is
% an identity matrix.
% - Compute group means on the transformed variables.
% - Perform a principal component analysis on the means, weighting each mean by
% the number of observations in the group. The eigenvalues are equal to the ratio
% of between-group variation to the within-group variation in the direction of
% each principal component. Here, the principal component analysis is runned by
% the singular value decomposition.
% - Back-transform the principal components into the space of the original variables,
% obtaining the canonical variables.
%
%File gives you the option to get an unbiased or maximum-likelihood parameter estimation.
%
% Syntax: function [RAFisher2cda] = RAFisher2cda(X,pp,c,alpha)
%
% Inputs:
% X - multivariate data matrix.
% pp - vector of prior probabilities (unknown,pp = 1;
% known [default], pp = 2 [you must to give it]).
% c - asking for ellipse confidence bounds (yes:c = 1;
% [default], no:c = 2).
% alpha - significance level (default = 0.05). Used for the
% Chi-square tests and for the ellipse confidence bounds.
%
% Output:
% - Canonical discriminant functions. It includes:
% constant, variates, eigenvalues, percentages and
% cumulative percentages. Also the Bartlett's approximate
% chi-squared statistic for testing the canonical correlation
% coefficients. Optionally, it also can show you the pair-wise
% canonical scores with its ellipses confidence bounds.
%
% Example: From the Table 11.5 (Iris data) of Johnson and Wichern (1992, p. 562),
% with 150 observations (n = 150), four variables (p = 4) [all = cm] and
% three groups (g = 3; sizes: 50,50,50). We are interested to apply a
% Canonical Discriminant Analysis with an unbiased parameter estimation.
% It is considering equal prior probabilities (pp = 1) and the groups
% ellipse bound with a probability of 0.95 (alpha value of 0.05).
%
% ------------------------------------------------------------------------------
% 1 2 3
% ------------------------------------------------------------------------------
% x1 x2 x3 x4 x1 x2 x3 x4 x1 x2 x3 x4
% ------------------------------------------------------------------------------
% 5.1 3.5 1.4 0.2 7.0 3.2 4.7 1.4 6.3 3.3 6.0 2.5
% 4.9 3.0 1.4 0.2 6.4 3.2 4.5 1.5 5.8 2.7 5.1 1.9
% 4.7 3.2 1.3 0.2 6.9 3.1 4.9 1.5 7.1 3.0 5.9 2.1
% 4.6 3.1 1.5 0.2 5.5 2.3 4.0 1.3 6.3 2.9 5.6 1.8
% 5.0 3.6 1.4 0.2 6.5 2.8 4.6 1.5 6.5 3.0 5.8 2.2
% 5.4 3.9 1.7 0.4 5.7 2.8 4.5 1.3 7.6 3.0 6.6 2.1
% . . . . . . . . . . . .
% . . . . . . . . . . . .
% . . . . . . . . . . . .
% 5.1 3.8 1.6 0.2 5.7 2.9 4.2 1.3 6.3 2.5 5.0 1.9
% 4.6 3.2 1.4 0.2 6.2 2.9 4.3 1.3 6.5 3.0 5.2 2.0
% 5.3 3.7 1.5 0.2 5.1 2.5 3.0 1.1 6.2 3.4 5.4 2.3
% 5.0 3.3 1.4 0.2 5.7 2.8 4.1 1.3 5.9 3.0 5.1 1.8
% ------------------------------------------------------------------------------
%
% Group 1 = Iris setosa (n1 = 50); Group 2 = Iris versicolor (n2 = 50);
% Group 3 = Iris virginica (n3 = 50).
% Var1 (x1) = sepal length; Var2 (x2) = sepal with;
% Var3 (x3) = petal length; Var4 (x4) = petal with.
%
% Total data matrix must be:
% You can get the X-matrix by calling to iris data file provided in
% the zip as
% load path-drive:irisdata
%
% Calling on Matlab the function:
% RAFisher2cda(X,1)
%
% Answer is:
%
% It was asking for the ellipses confidence bounds.
% Do you want an unbiased (1) or maximum-likelihood parameter estimation? (2):1
%
% Canonical Discriminant Functions.
% --------------------------------------------------------------------------------
% Constant =
% -2.1051 6.6615
%
% Variates =
% 0.8294 0.0241
% 1.5345 2.1645
% -2.2012 -0.9319
% -2.8105 2.8392
%
% Eigenvalue =
% 32.1919 0.2854
%
% Percentage =
% 99.1213 0.8787
%
% CumulativePercentage =
% 99.1213 100.0000
% --------------------------------------------------------------------------------
% Functions = columns. On variates, Variate1 = first row and so forth to 4
%
% Chi-square Tests with Successive Roots Removed.
% -------------------------------------------------------------------------
% Removed Eigenvalue CanCor LW Chi-sqr. df P
% -------------------------------------------------------------------------
% 32.1919 0.9848 0.0234 546.1153 8 0.0000
% 1 0.2854 0.4712 0.7780 36.7807 3 0.0000
% -------------------------------------------------------------------------
% With a given significance of: 0.05
% If P-value >= alpha, results not significative. Else, it is significative.
%
% The pair-wise plots you can get are: 1
% --------------
% plots =
%
% 1 2
%
% --------------
% Give me the pair of factors to plot [a b]: [1 2]
% Do you need another plot? (y/n): n
%
%
% Created by A. Trujillo-Ortiz, R. Hernandez-Walls and S. Perez-Osuna
% Facultad de Ciencias Marinas
% Universidad Autonoma de Baja California
% Apdo. Postal 453
% Ensenada, Baja California
% Mexico.
% atrujo@uabc.mx
%
% Copyright (C) April 15, 2004.
%
% To cite this file, this would be an appropriate format:
% Trujillo-Ortiz, A., R. Hernandez-Walls and S. Perez-Osuna. (2004).
% RAFisher2cda:Canonical Discriminant Analysis. A MATLAB file.
% [WWW document]. URL http://www.mathworks.com/matlabcentral/fileexchange/
% loadFile.do?objectId=4836
%
% References:
%
% Johnson, R. A. and Wichern, D. W. (1992), Applied Multivariate Statistical Analysis.
% 3rd. ed. New-Jersey:Prentice Hall. Chapter 11, pp. 493-572.
% Fisher, R. A. (1936), The use of multiple measurements in taxonomic problems. Annals
% of Eugenics, 7: 179-188.
%
if nargin < 4,
alpha = 0.05; %(default)
end;
if (alpha <= 0 | alpha >= 1)
fprintf('Warning: significance level must be between 0 and 1\n');
return;
end;
if nargin < 3,
c = 1; %(default)
end;
disp(' ')
if c == 1;
fprintf('It was asking for the ellipses confidence bounds.');
else c == 2;
fprintf('It was not asking for the ellipses confidence bounds.');
end;
disp(' ')
if nargin < 2,
pp = 1; %(default)
end;
if nargin < 1,
error('Requires at least one input arguments.');
return;
end;
Y = X(:,2:end);
G = X(:,1);
g = max(G); %number of groups
%Vector of sample sizes per group.
n = [];
indice = X(:,1);
for i = 1:g
Xe = find(indice==i);
eval(['X' num2str(i) '= X(Xe,2);']);
eval(['n' num2str(i) '= length(X' num2str(i) ') ;'])
eval(['xn = n' num2str(i) ';'])
n = [n,xn];
end;
r = 1;
r1 = n(1);
%Partition of the group mean and covariance matrices.
M = [];
for k = 1:g
eval(['M' num2str(k) '= mean(Y(r:r1,:));']);
eval(['m = M' num2str(k) ';'])
if k < g
r = r+n(k);
r1 = r1+n(k+1);
end;
M = [M;m];
end;
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.')
disp('--------------------------------------------------------------------------------\');
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);
disp('--------------\');
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;
