%%%%%%%   correspondence analysis(ANACOR)  %%%%%%%%%%%%%%%%%%%%%
% X為原始數(shù)據(jù) p X n維 即有n個樣本，每個樣本由p個變量來描述。返回F為R型因子分析后的結(jié)果，Q為最后的結(jié)果,A=ZZ'。
%function [G,F,A]=ANACOR(X)
function [G,F]=ANACOR(X)
%為了確保表示的向量大小與表示樣本的向量大小有同樣的意義，對原始數(shù)據(jù)進(jìn)行處理
T=sum(sum(X));
Xi=sum(X');  %% 對行求和，即：變量的邊緣概率
Xj=sum(X);   %% 對列求和，即：樣本的邊緣概率
%%% 計算矩陣Z=（zij）

Z=(X-(Xi'*Xj/T))./sqrt(Xi'*Xj);
A=Z*Z';
%%%% 求特征值和特征向量,A為輸入方陣，D為由A的全部特征值構(gòu)成的對角矩陣，V的各列為對應(yīng)于特征值的特征向量構(gòu)成的矩陣
%[V,D]=eig(A);
D=eig(A,'nobalance');
% 特征值降序排列，并求累積量>=85%
b=sort(D);
D=b(end:-1:1);
d=sum(D);
D1=D/d;
i=1;
D2=D1(1);
while D2<=0.95
    i=i+1;
    D2=D2+D1(i);
end
%求前i項特征值對應(yīng)的特征向量
option = struct('disp',0);
[eigvector, eigvalue] = eigs(A,i,'la',option);  %'LA' or 'SA' - Largest or Smallest Algebraic

%%% R型因子分析

U=eigvector;  %%修正后的特征值對應(yīng)的特征向量，p X i維， i為使累積量>=85%的特征值數(shù)量
[m,n]=size(U);
U1=zeros(m,n);
E=D(1:i);     %%修正后的特征值，  在此為行向量
for x=1:i
    U1(:,x)=U(:,x)*sqrt(E(x));
end
%%% 求載荷矩陣 F
F=U1;
%[m,n]=size(F);
%% plot(F,'r*');
%% 畫圖
%for i=1:m
%    plot(F(i,1),F(i,2),'r*');
%    hold on;
%end

%%%求Q型因子分析

%B=Z'*Z;
V=Z'*U;
[m,n]=size(V);
V1=zeros(m,n);
E=D(1:i);     %%修正后的特征值，  在此為行向量
for x=1:i
    V1(:,x)=V(:,x)*sqrt(E(x));
end
G=V1;
%% 畫圖
[m,n]=size(F);
for i=1:m
    plot(F(i,1),F(i,2),'r*');
    hold on;
end

hold on

[x,y]=size(G);
for i=1:x
    plot(G(i,1),G(i,2),'bo');
    hold on;
end