1.题目:假定某地固定资产投资率x1、通货膨胀率x2、失业率x3、相关系数矩阵为
,试用主成分分析法求因子分析模型
(《Python数学实验与建模》(司守奎)例11.9)
2.代码如下
注意:
(1)要先明白因子分析中的基本概念,如贡献、贡献率、共同度、载荷矩阵
(2)要先安装factor_analyzer库,在Anaconda Prompt里pip install factor_analyzer,这个库的作用就是根据相关系数矩阵,旋转因子得到理想的载荷矩阵
版本声明:
python 3.7.1
#特征值要降序,对应的特征向量也要排序
#我参考了
import numpy as np
from sklearn import decomposition as dc
from scipy.linalg import eig
from factor_analyzer import FactorAnalyzer#这个包需要自己安装,用于因子旋转以得到更好结果,见p345
#1.应用原理计算方法
r=np.array([[1,1/5,-1/5],
[1/5,1,-2/5],
[-1/5,-2/5,1]])
val,vec=np.linalg.eig(r)#求相关系数矩阵的特征值和特征向量,但是这样求出来的特征值是没有排序的
index=np.argsort(val)[::-1]
print('特征值降序的下标:',index)
A1=np.tile(np.sqrt(val),(3,1))*vec#利用同维数矩阵逐个元素相乘求载荷矩阵
A2=vec*np.sqrt(val)#向量放第二位自动转置,好像也可以解释成广播机制,不过我更倾向于第一种理解
A3=A2[:,index]#排序后的载荷矩阵
print('排序前特征值:',val.round(4))
print('排序后特征值:',np.sort(val)[::-1].round(4))
print('排序前特征向量:\n',vec.round(4))
print('排序后特征向量:\n',vec[:,index].round(4))
print('排序前载荷矩阵:\n',A1.round(4))
print('排序后载荷矩阵:\n',A3.round(4))
num=int(input(('请输入选择公共因子的个数:')))
A=A3[:,:num]#提出num个因子的载荷矩阵
#print('\n',A)
Ac=np.sum(A**2,axis=0)#逐列元素求和
Ar=np.sum(A**2,axis=1)#逐行元素求和
print('对x的贡献为:',Ac)
print('对x的贡献率为:',Ac/3)
print('对x的累计贡献率为:',np.sum(Ac/3))
print('共同度为:',Ar)
#根据贡献率来看,可得知因子f1是最重要的
#共同度是指某一原变量在所有公因子上载荷平方和反应公公元素对原变量的方差的解释程度,大于0.8即可得到很好的结果
#2.应用factor_analyzer求解
f=FactorAnalyzer(rotation='varimax', n_factors=num, method='principal')
'''
(1)rotation:旋转的方式,包括None:不旋转,'varimax':最大方差法,'promax':最优斜交旋转;
(2)n_factors:公因子的数量;
(3)method:因子分析的方法,包括'minres':最小残差因子法,'principal':主成分分析法;
'''
f.fit(r)
print('\n因子旋转后的载荷矩阵:\n',f.loadings_)
Ac=np.sum(f.loadings_**2,axis=0)#逐列元素求和
Ar=np.sum(f.loadings_**2,axis=1)#逐行元素求和
print('对x的贡献为:',Ac)
print('对x的贡献率为:',Ac/3)
print('对x的累计贡献率为:',np.sum(Ac/3))
print('共同度为:',Ar)
输入选择几个公共因子,当num=2时,结果为:
特征值降序的下标: [1 0 2]
排序前特征值: [0.8536 1.5464 0.6 ]
排序后特征值: [1.5464 0.8536 0.6 ]
排序前特征向量:
[[-0.8881 0.4597 0. ]
[ 0.3251 0.628 0.7071]
[-0.3251 -0.628 0.7071]]
排序后特征向量:
[[ 0.4597 -0.8881 0. ]
[ 0.628 0.3251 0.7071]
[-0.628 -0.3251 0.7071]]
排序前载荷矩阵:
[[-0.8205 0.5717 0. ]
[ 0.3003 0.7809 0.5477]
[-0.3003 -0.7809 0.5477]]
排序后载荷矩阵:
[[ 0.5717 -0.8205 0. ]
[ 0.7809 0.3003 0.5477]
[-0.7809 -0.3003 0.5477]]
请输入选择公共因子的个数:2
对x的贡献为: [1.54641016 0.85358984]
对x的贡献率为: [0.51547005 0.28452995]
对x的累计贡献率为: 0.7999999999999998
共同度为: [1. 0.7 0.7]
因子旋转后的载荷矩阵:
[[ 0.17265918 0.98498163]
[ 0.99698659 0.07757406]
[-0.85195502 -0.52361497]]
对x的贡献为: [1.74962082 1.25037918]
对x的贡献率为: [0.58320694 0.41679306]
对x的累计贡献率为: 1.0
共同度为: [1. 1. 1.]
对比可发现旋转因子后,共同度变大了,反应公共因素对原变量的方差的解释效果很好
factor_analyzer学习参考