#此代碼是我手動(dòng)搭建的

#除KMO和bar檢驗(yàn)外未使用因子分析的專用工具包(factor_analyzer)

#此代碼僅能解決需要輸出兩個(gè)潛變量的因子分析問題

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

import numpy as np
import matplotlib.pyplot as plt
from sklearn import preprocessing
import matplotlib;matplotlib.rc("font",family='Microsoft YaHei')
from sklearn.preprocessing import normalize
from factor_analyzer import calculate_kmo,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
# 對(duì)樣本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]))
#求Z-score標(biāo)準(zhǔn)化后的樣本間的協(xié)方差矩陣C
C=np.dot(X.T,X)/(data.shape[0])
#求協(xié)方差矩陣C的特征值λ和特征向量V0
λ,V0=np.linalg.eig(C);V0=np.mat(V0)
#協(xié)方差矩陣C的特征向量施密特正交化和單位化
a=[V0[:,i] for i in range(V0.shape[1])]
def b(i):
? ?if i==0:return a[0]
? ?else:return a[i]-sum(float(a[i].T*b(k)/(b(k).T*b(k)))*b(k) for k in range(i))
V1=b(0)
for i in range(1,V0.shape[1]):V1=np.column_stack((V1,b(i)))
V1=normalize(np.array(V1),axis=0,norm='l2')
# 按照協(xié)方差矩陣C的特征值大小對(duì)特征向量排序
tezhenggen1=[i for i in λ];tezhenggen1.sort(reverse=True)
V=np.zeros(shape=(data.shape[1],data.shape[1]))
for i in range(data.shape[1]):
? ?for j in range(data.shape[1]):V[i,j]=V1[i,list(λ).index(tezhenggen1[j])]
# 計(jì)算載荷矩陣A(兩列)和得分系數(shù)矩陣S0(兩列)
Z=np.mat(np.zeros(shape=(data.shape[1],data.shape[1])))
for i in range(data.shape[1]):Z[i,i]=1/(tezhenggen1[i]**0.5)
A,S0=(V*Z.I)[:,0:2],(V*Z)[:,0:2]
# 對(duì)載荷矩陣A施以二維正交旋轉(zhuǎn)
#(1)定義載荷間的總方差函數(shù)ff(θ)
def ff(a):
? ?B=np.mat(A)*np.mat([[np.cos(a),-np.sin(a)],[np.sin(a),np.cos(a)]])
? ?zongfangcha=0
? ?for k in range(B.shape[1]):
? ? ? ?junzhi=sum(B[i,k]**2/sum(A[i,j]**2 for j in range(A.shape[1]))
? ? ? ? ? ? ? ? ? ? for i in range(B.shape[0]))/B.shape[0]
? ? ? ?zongfangcha+=sum((B[i,k]**2/sum(A[i,j]**2
? ? ? ? ? ? ? ? ? ?for j in range(A.shape[1]))-junzhi)**2
? ? ? ? ? ? ? ? ? ?for i in range(B.shape[0]))
? ?return zongfangcha
#(2)計(jì)算并修正載荷矩陣的旋轉(zhuǎn)解和得分系數(shù)矩陣
x=np.arange(0,np.pi/2,0.01);y=[ff(i) for i in x];θ=x[y.index(max(y))]
for i in range(4):
? ?θ+=i*np.pi/2
? ?Q=np.mat([[np.cos(θ),-np.sin(θ)],[np.sin(θ),np.cos(θ)]])
? ?B=A*Q
? ?list0=[abs(B[j,0]) for j in range(B.shape[0])]
? ?list1=[abs(B[j,1]) for j in range(B.shape[0])]
? ?a,b=B[list0.index(max(list0)),0],B[list1.index(max(list1)),1]
? ?if a>0 and b>0:break
? ?else:pass
Q=np.mat([[np.cos(θ),-np.sin(θ)],[np.sin(θ),np.cos(θ)]])
B=A*Q;S=S0*Q
list2=[(B.T*B)[j,j] for j in range(B.shape[1])]
if list2.index(max(list2))==0:pass
else:B=B*np.mat([[0,1],[1,0]]);S=S*np.mat([[0,1],[1,0]])
list2.sort(reverse=True)
#輸出旋轉(zhuǎn)角度θ和載荷間的總方差最大值max(ff)
print('-'*54)
print('旋轉(zhuǎn)角度θ='+str('%.2f'%(θ*180/np.pi))+'°',end='\t')
print('載荷間的總方差='+str('%.3f'%ff(θ)))
#輸出使載荷間的總方差函數(shù)達(dá)到最大的旋轉(zhuǎn)矩陣Q
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()
print('旋轉(zhuǎn)方法：Kaiser 標(biāo)準(zhǔn)化最大方差法')
#輸出載荷矩陣的旋轉(zhuǎn)解
print('-'*54);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(B.shape[0]):
? ?print('指標(biāo)'+str(i+1),end='\t')
? ?for j in range(B.shape[1]):print('%.3f'%B[i,j],end='\t')
? ?print('%.3f'%sum(B[i,k]**2 for k in range(B.shape[1])),end='\t')
? ?print()
print('總貢獻(xiàn)',end='\t')
for i in list2:print('%.3f'%i,end='\t')
print('%.3f'%sum(list2));print('-'*53)
#輸出旋轉(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('-'*54)
print('各樣本的得分及排序');print('樣本|fi',end='\t')
for i in range(A.shape[1]):print('潛變量f'+str(i+1),end='\t')
print('綜合變量f',end='\t');F=X*S
print('按f1排序',end='\t');print('按f2排序',end='\t');print('按f排序')
def xvhao(x):
? ?list3=[i for i in x]
? ?list3.sort(reverse=True)
? ?return [list3.index(i)+1 for i in x]
list4=[i/sum(list2) for i in list2]
list5=[sum(list4[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'%list5[i],' ?',end='\t');print(xvhao(F[:,0])[i],' ?',end='\t')
? ?print(xvhao(F[:,1])[i],' ?',end='\t');print(xvhao(list5)[i],' ?',end='\t')
? ?print()
print('-'*54)
#繪制載荷圖(旋轉(zhuǎn)后的空間組件圖)
plt.plot(B[:,0],B[:,1],'o');plt.xlabel('潛變量f1');plt.ylabel('潛變量f2')
plt.title('載荷圖(旋轉(zhuǎn)后的空間組件圖)')
for i in range(B.shape[0]):plt.text(B[i,0]+0.01,B[i,1]+0.01,'指標(biāo)'+str(i+1))
plt.show()