mcorran2.m

这是一个matlab版的基础响应分析功能文件

%  Are you interested to plot the origin reference of the variable levels? (y/n): n
%  Are you interested to plot the origin quadrature? (y/n): y
%  
%  Do you need another plot? (y/n):n
%
%  Created by A. Trujillo-Ortiz, R. Hernandez-Walls and K. Barba-Rojo
%             Facultad de Ciencias Marinas
%             Universidad Autonoma de Baja California
%             Apdo. Postal 453
%             Ensenada, Baja California
%             Mexico.
%             [e-mail]
%  Copyright. December 24, 2008.
%
%  References:
%  Greenacre, M.J. (2006), From simple to multiple correspondence analysis,
%       Chapter 2, in Greenacre, M.J. and Blasius, J. (eds.), Multiple 
%       Correspondence Analysis and Related Methods, Chapman-Hall/CRC:London,
%       pp.41-76.
%  Rencher, A. C. (2000), Methods of Multivariate Analysis. 2nd ed. 
%       John Wiley & Sons. p. 514-525
%

g = input('Give me the vector of categories of each variable. Please, use [ ]: ');
p = length(g); %number of variables
b1 = size(X,2); %number of categories
b2 = sum(g); %input number of categories
if b1 ~= b2;
    disp('Warning:the number of categories in the input vector doesn''t correspond to the given in the Burt matrix. Please, check.');
    return
end
q = sum(g)-p;
B = X;
r = sum(B,2); %row total
Dr = diag(r);
Z = (sqrt(inv(Dr))*B*sqrt(inv(Dr)));
[U,S,V] = svd(Z); %S are eigenvalues (= singular values)
I = diag(S).^2; %principal inertias
I = I(2:q+1); %principal inertias equal to one and zero are deleted
PI = I./sum(I)*100; %percent of inertia
PA = cumsum(PI); %cumulative percent of inertia
N = sum(sum(B)); %grand total
P = B/N; %correspondence matrix
r = sum(P); %row totals (masses equal to column masses)

%Adjustment of Inertias
SV = sqrt(I);
adj = SV>=(1/p); %adjustment of principal inertias
SVad = SV(adj);
SVadj =(p/(p-1))^2*(SVad - (1/p)).^2;
EVm = (p/(p-1))*(sum(I)-(q/(p^2)));
PIadj = SVadj/EVm*100; %percent of explained adjusted inertia
PAadj = cumsum(PIadj); %cumulative percent of adjusted inertia

a = []; %i-th row (or column) standard coordinates for the s-th dimension
for i = 1:length(r)
    ai = V(:,i)./sqrt(r');
    a = [a ai];
end

D = diag(S);
f = []; %principal coordinates for the s-th dimension
for i = 1:length(r)
    fi = a(:,i).*D(i);
    f = [f fi];
end
a = a(:,2:q+1);
D = D(2:q+1);
disp(' ')
disp('Inertias of the Multiple Correspondence Analysis')
disp('for the Burt matrix given.')
fprintf('----------------------------------------\n');  
disp(' Eigenvalue     Percent     Cummulative');           
disp('               (Inertia)      Percent');
fprintf('----------------------------------------\n');
fprintf('%9.4f   %10.2f    %10.2f\n',[I,PI,PA].');
fprintf('----------------------------------------\n');
fprintf('  Total    %11.2f\n', sum(PI));
disp(['Variable categories = ',num2str(g)]);
fprintf('----------------------------------------\n');
disp(' ');
disp('Adjusted inertias of the Multiple Correspondence')
disp('Analysis for the Burt matrix given.')
fprintf('----------------------------------------\n');  
disp(' Eigenvalue     Percent     Cummulative');           
disp('               (Inertia)      Percent');
fprintf('----------------------------------------\n');
fprintf('%9.5f   %10.2f    %10.2f\n',[SVadj,PIadj,PAadj].');
fprintf('----------------------------------------\n');
fprintf('  Total    %11.2f\n', sum(PIadj));
disp(['Variable categories = ',num2str(g)]);
fprintf('----------------------------------------\n');
disp(' ');

Dm = [];
for i = 1:q
    for s = i+1:q
        if s ~= i
            d = [i s];
            Dm = [Dm;d];
        end
    end
end

co = (q*(q - 1))/2;
gf = input('Are you interested to get the dimensions plots? (y/n): ','s');
if gf == 'y'
    disp('--------------');
    fprintf('The pair-wise plots you can get are: %.i\n', co);
    disp('--------------');
    Dm
    fprintf('--------------\n');
    
    d = 'y';
    while d == 'y';
        p = input('Give me the interested dimensions to plot. Please, use [a b]:');
        figure;
        Y = a;
        pol = 0;
        hold on
        lg = [];
        for c = 1:length(g),
            plot(Y(1+pol:g(c)+pol,p(1)),Y(1+pol:g(c)+pol,p(2)),'*','Color',paletc(c));
            lg = [lg,['''Levels of categorical variable ' num2str(c) ''',']];
            A1 = Y(1+pol:g(c)+pol,p(1));
            A2 = Y(1+pol:g(c)+pol,p(2));
            for Kl = 1:length(A1),
                text(A1(Kl)+.02,A2(Kl)+.02,num2str(Kl));
            end
            pol = g(c)+pol;
        end
        lg(end) = ' ';
        eval(['legend(' lg ')']);     
        title({'Plot of variable levels of the performed Multiple Correspondence Analysis';'from the Burt matrix.'});      
        xlabel(['Dimension ' num2str(p(1))';]);
        ylabel(['Dimension ' num2str(p(2))';]);
        
        disp(' ')
        vl = input('Are you interested to plot the origin reference of the variable levels? (y/n): ','s');
        if vl == 'y'     
            disp(' ')
            disp('Note: At the end of the program execution. If you are interested to fit the variable labels')
            disp('on the generated figures you can turn-on active button ''Edit Plot'', do click on the')
            disp('selected label and drag to fix it on the desired position. Then turn-off active ''Edit Plot''.')
            disp(' ')
            pol = 0;
            for c = 1:length(g),
                hold on
                A1 = Y(1+pol:g(c)+pol,p(1));
                A2 = Y(1+pol:g(c)+pol,p(2));
                plot2org(A1,A2,'--k'), %this m-file was taken from Jos's plot2org
                pol = g(c)+pol;
            end
        else
        end
        o = input('Are you interested to plot the origin quadrature? (y/n): ','s');
        if o == 'y'     
            hline(0,'k'); %this m-file was taken from Brandon Kuczenski's hline and vline zip
            vline(0,'k'); %this m-file was taken from Brandon Kuczenski's hline and vline zip
            disp(' ')
            %d=input('Do you need another plot? (y/n):','s');
        else            
        end
        d=input('Do you need another plot? (y/n):','s');
    end
else
end

function res = paletc(n)
% - utility program for selecting plot colors
while (n>10)
   n = n - 10;
end
% build cell-array of colors to choose from:
colors{1}  = [0 0 1];          % blue
colors{2}  = [1 0 0];          % red
colors{3}  = [0 1 0];          % green
colors{4}  = [0.5 0 0.5];      % purple
colors{5}  = [0 0.5 0];        % dark green
colors{6}  = [0.28 0.73 0.94]; % sky
colors{7}  = [0 0 0];          % black
colors{8}  = [0 0 0.5];        % navy
colors{9}  = [1 0.5 0];        % orange
colors{10} = [1 0 1];          % magenta
colors{11} = [1 1 0];          % yellow
res = colors{n};

return,