%  从相关系数矩阵出发进行因子分析

%***************************定义相关系数矩阵PHO*****************************

PHO = [1     0.79    0.36    0.76    0.25    0.51

0.79  1       0.31    0.55    0.17    0.35

0.36  0.31    1       0.35    0.64    0.58

0.76  0.55    0.35    1       0.16    0.38

0.25  0.17    0.64    0.16    1       0.63

0.51  0.35    0.58    0.38    0.63    1

PHO=xlsread('fa1.xls',2)

%******************调用factoran函数根据相关系数矩阵作因子分析*****************

% 从相关系数矩阵出发，进行因子分析(不进行因子旋转)，公共因子个数是2，delta设定特殊方差psi的下界

% ‘covariance’表示协方差或相关系数矩阵，'rotate'='none'不旋转

% 输出的lambda是因子载荷L矩阵，psi是估计的个体方差，T为旋转矩阵

[lambda,psi,T] = factoran(PHO,2,'xtype','covariance','delta',0,'rotate','none')

% 定义元胞数组，以元胞数组形式显示结果

head = {'变量', '因子f1', '因子f2'};

varname = {'身高','坐高','胸围','手臂长','肋围','腰围','',''}';

Contribut = 100*sum(lambda.^2)/6;

CumCont = cumsum(Contribut);

result1 = num2cell([lambda; Contribut; CumCont]);

result1 = [head; varname, result1]

biplot(lambda,'LineWidth',2,'MarkerSize',20)

scatter(lambda(:,1),lambda(:,2))

%从相关系数矩阵出发，进行因子分析(进行因子旋转)，默认是使用最大方差旋转法

[lambda,psi,T] =factoran(PHO,2,'xtype','covariance','delta',0)

Contribut = 100*sum(lambda.^2)/6

CumCont = cumsum(Contribut)

biplot(lambda,'LineWidth',2,'MarkerSize',20)

scatter(lambda(:,1),lambda(:,2))

% 从相关系数矩阵出发，进行因子分析，公共因子数为3(进行因子旋转)

[lambda,psi,T] = factoran(PHO,3,'xtype','covariance','delta',0)

Contribut = 100*sum(lambda.^2)/6

CumCont = cumsum(Contribut)

biplot(lambda,'LineWidth',2,'MarkerSize',20)

% 不建议做m>=4，报错

[lambda,psi,T] = factoran(PHO,4,'xtype','covariance','delta',0)

%*********************************读取数据*********************************

[X,textdata] = xlsread('fa2.xls');

X = X(:,3:end);

varname = textdata(4,3:end);

obsname = textdata(5:end,2);

%******************调用factoran函数根据原始观测数据作因子分析*****************

% 从原始数据(实质还是相关系数矩阵)出发，进行因子分析，公共因子数为4

% 进行因子旋转(最大方差旋转法)，stats返回一些拟合的检验统计量

[lambda,psi,T,stats] = factoran(X,4)

Contribut = 100*sum(lambda.^2)/8

CumCont = cumsum(Contribut)

% 从原始数据(实质还是相关系数矩阵)出发，进行因子分析，公共因子数为2

% 进行因子旋转(最大方差旋转法)，F返回因子得分

[lambda,psi,T,stats,F] = factoran(X, 2)

Contribut = 100*sum(lambda.^2)/8

CumCont = cumsum(Contribut)

[varname', num2cell(lambda)]

biplot(lambda,'LineWidth',2,'MarkerSize',20)

plot(lambda(:,1),lambda(:,2),'k*')

%**************将因子得分F分别按耐力因子得分和速度因子得分进行排序*************

obsF = [obsname, num2cell(F)];

F1 = sortrows(obsF, 2);    % 按耐力因子得分排序，默认按列排序

F2 = sortrows(obsF, 3);    % 按速度因子得分排序

head = {'国家/地区','耐力因子','速度因子'};

result1 = [head; F1];

result2 = [head; F2];

%*************************绘制因子得分负值的散点图***************************

plot(-F(:,1),-F(:,2),'k.');

xlabel('耐力因子得分(负值)');

ylabel('速度因子得分(负值)');