#此代碼是我基于因子分析的專用工具包(factor_analyzer)搭建的

#此代碼能解決需要輸出多個(gè)(n_factors=1,2,3,...m)潛變量的因子分析問題

#如需將此代碼放在自己個(gè)人社交平臺(tái)，請給個(gè)引用"B站 耿大哥講算法"?

import numpy as np
from sklearn import preprocessing
from factor_analyzer import calculate_kmo
from factor_analyzer import FactorAnalyzer
from factor_analyzer import calculate_bartlett_sphericity
#讀取數(shù)據(jù)&校準(zhǔn)化
# ? ? ? ? ? ? ? 指標(biāo)1 指標(biāo)2 指標(biāo)3 指標(biāo)4指標(biāo)5 指標(biāo)6 指標(biāo)7 指標(biāo)8
data=np.array([[40.4,24.7,7.20,6.10,8.30,8.70,2.442,20.0], #樣本1
? ? ? ? ? ? ? [25.0,12.7,11.2,11.0,12.9,20.2,3.542,9.10], #樣本2
? ? ? ? ? ? ? [13.2,3.30,3.90,4.30,4.40,5.50,0.578,3.60], #樣本3
? ? ? ? ? ? ? [22.3,6.70,5.60,3.70,6.00,7.40,0.176,7.30], #樣本4
? ? ? ? ? ? ? [34.3,11.8,7.10,7.10,8.00,8.90,1.726,27.5], #樣本5
? ? ? ? ? ? ? [35.6,12.5,16.4,16.7,22.8,29.3,3.017,26.6], #樣本6
? ? ? ? ? ? ? [22.0,7.80,9.90,10.2,12.6,17.6,0.847,10.6], #樣本7
? ? ? ? ? ? ? [48.4,13.4,10.9,9.90,10.9,13.9,1.772,17.8], #樣本8
? ? ? ? ? ? ? [40.6,19.1,19.8,19.0,29.7,39.6,2.449,35.8], #樣本9
? ? ? ? ? ? ? [24.8,8.00,9.80,8.90,11.9,16.2,0.789,13.7], #樣本10
? ? ? ? ? ? ? [12.5,9.70,4.20,4.20,4.60,6.50,0.874,3.90], #樣本11
? ? ? ? ? ? ? [1.80,0.60,0.70,0.70,0.80,1.10,0.056,1.00], #樣本13
? ? ? ? ? ? ? [32.3,13.9,9.40,8.30,9.80,13.3,2.126,17.1], #樣本13
? ? ? ? ? ? ? [38.5,9.10,11.3,9.50,12.2,16.4,1.327,11.6]])#樣本14
# 對樣本Z-score標(biāo)準(zhǔn)化
X=preprocessing.scale(data,axis=0)
#KMO和barlett檢驗(yàn)
bartlett,KMO=calculate_bartlett_sphericity(data),calculate_kmo(data)[1]
print('-'*54);print('KMO和bartlett檢驗(yàn)')
print('KMO='+str('%.3f'%calculate_kmo(data)[1]),end='\t')
if ? KMO<=1 ?and KMO>=0.9:print('非常好')
elif KMO<0.9 and KMO>=0.8:print('好')
elif KMO<0.8 and KMO>=0.7:print('一般')
elif KMO<0.7 and KMO>=0.6:print('差')
elif KMO<0.6 and KMO>=0.5:print('很差')
else:print('不能接受')
print('bartlett 統(tǒng)計(jì)量='+str('%.3f'%bartlett[0]))
print('bartlett 自由度='+str(int(data.shape[1]*(data.shape[1]-1)/2)))
print('bartlett P值='+str('%.3f'%bartlett[1]))
#(常用)method='minres'迭代主因子,'principal'主成分
#(不常用)method='ml'極大似然,'mle'極大似然估計(jì),'uls'無條件最小二乘,
#(常用)rotation=varimax標(biāo)準(zhǔn)化最大方差法
#(不常用)rotation=promax,oblimin,oblimax,quartimin,quartimax,equamax
FA=FactorAnalyzer(n_factors=2,method='minres',rotation='varimax')
FA.fit(X)
#計(jì)算載荷矩陣
A=FA.loadings_
#計(jì)算樣本得分矩陣
F=FA.transform(X)
#計(jì)算得分系數(shù)矩陣
S=np.mat(X).I*F
#計(jì)算共同度
h=FA.get_communalities()
#計(jì)算方差貢獻(xiàn)
g=FA.get_factor_variance()
#計(jì)算旋轉(zhuǎn)矩陣
Q=FA.rotation_matrix_
# 輸出使載荷間的總方差函數(shù)達(dá)到最大的旋轉(zhuǎn)矩陣Q
print('-'*50)
print('旋轉(zhuǎn)矩陣Q')
for i in range(Q.shape[0]):
? ?for j in range(Q.shape[1]):print('%.3f'%Q[i,j],end='\t')
? ?print()
# 輸出載荷矩陣的旋轉(zhuǎn)解
print('-'*50);print('載荷矩陣的旋轉(zhuǎn)解')
print('指標(biāo)|組件',end='\t')
for i in range(A.shape[1]):print('潛變量f'+str(i+1),end='\t')
print('共同度')
for i in range(A.shape[0]):
? ?print('指標(biāo)'+str(i+1),end='\t')
? ?for j in range(A.shape[1]):print('%.3f'%A[i,j],end='\t')
? ?print('%.3f'%h[i],end='\t')
? ?print()
print('總貢獻(xiàn)',end='\t')
for i in g[0]:print('%.3f'%i,end='\t')
print('%.3f'%sum(g[0]))
print('貢獻(xiàn)率',end='\t')
for i in g[1]:print('%.1f'%(i*100)+'%',end='\t')
print('-----')
print('貢獻(xiàn)累積',end='\t')
for i in g[2]:print('%.1f'%(i*100)+'%',end='\t')
print('-----')
print('-'*50)
#輸出旋轉(zhuǎn)后的成分得分系數(shù)矩陣
print('成分得分系數(shù)矩陣(旋轉(zhuǎn)后)')
print('指標(biāo)|組件',end='\t')
for i in range(A.shape[1]):print('潛變量f'+str(i+1),end='\t')
print()
for i in range(S.shape[0]):
? ?print('指標(biāo)'+str(i+1),end='\t')
? ?for j in range(S.shape[1]):print('%.3f'%S[i,j],end='\t')
? ?print()
#計(jì)算各樣本的綜合得分并排序
print('-'*50)
print('各樣本的得分及排序');print('樣本|fi',end='\t')
for i in range(A.shape[1]):print('潛變量f'+str(i+1),end='\t')
print('綜合變量f',end='\t')
for i in range(A.shape[1]):print('按f'+str(i+1)+'排序',end='\t')
print('按f排序')
def xvhao(x):
? ?list2=[i for i in x]
? ?list2.sort(reverse=True)
? ?return [list2.index(i)+1 for i in x]
list3=[i/sum(g[1]) for i in g[1]]
list4=[sum(list3[j]*F[i,j] for j in range(F.shape[1])) for i in range(F.shape[0])]
for i in range(F.shape[0]):
? ?print('樣本'+str(i+1),end='\t')
? ?for j in range(F.shape[1]):print('%.3f'%F[i,j],end='\t')
? ?print('%.3f'%list4[i],' ?',end='\t')
? ?for j in range(F.shape[1]):print(xvhao(F[:,j])[i],' ?',end='\t')
? ?print(xvhao(list4)[i],' ?',end='\t')
? ?print()
print('-'*50)