#此代碼是我手動搭建的

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

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

#如需將此代碼放在自己個人社交平臺，請給個引用"B站 耿大哥講算法" 作者：耿大哥講算法?

import numpy as np
import matplotlib.pyplot as plt
from sklearn import preprocessing
from sklearn.preprocessing import normalize
import matplotlib;matplotlib.rc("font", family='Microsoft YaHei')
from factor_analyzer import calculate_kmo,calculate_bartlett_sphericity
#讀取數(shù)據(jù)&校準化
# ? ? ? ? ? ? ? 指標1 指標2 指標3 指標4指標5 指標6 指標7 指標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標準化
X=preprocessing.scale(data,axis=0)
#KMO和bartlett檢驗
bartlett,KMO=calculate_bartlett_sphericity(data),calculate_kmo(data)[1]
print('-'*20,'KMO和bartlett檢驗','-'*19)
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)計量='+str('%.3f'%bartlett[0]))
print('bartlett 自由度='+str(int(X.shape[1]*(X.shape[1]-1)/2)))
print('bartlett P值='+str('%.3f'%bartlett[1]))
X=preprocessing.scale(X,axis=0)
#求樣本協(xié)方差矩陣
C=np.dot(X.T,X)/(X.shape[0])
#主因子迭代25次(SPSS默認25次)
m,D,S0=2,np.eye(C.shape[0]),0
list1,list2,F,A=[[1+i/1000] for i in range(X.shape[1])],[1],0,0
for k in range(25):
? ?λ,Q0=np.linalg.eig(C-D);Q0=np.mat(Q0)
? ?Q0=normalize(np.array(Q0),axis=0)
? ?# 特征根和特征向量排序
? ?tezhenggen1=[i for i in λ];tezhenggen1.sort(reverse=True)
? ?Q=np.zeros(shape=(Q0.shape[1],m))
? ?for i in range(Q0.shape[1]):
? ? ? ?for j in range(m):Q[i,j]=Q0[i,list(λ).index(tezhenggen1[j])]
? ?#計算載荷矩陣和D
? ?diag=np.diag([tezhenggen1[i]**0.5 for i in range(m)])
? ?A=Q*np.mat(diag)
? ?D=np.diag([(C-A*A.T)[i,i] for i in range(Q0.shape[1])])
? ?list2.append(np.linalg.det(D))
? ?for i in range(X.shape[1]):list1[i].append(D[i,i])
? ?#print('第'+str(k+1)+'次迭代完成！')
? ?S0=((A.T*np.mat(D).I*A).I*A.T*np.mat(D).I).T
#繪制迭代過程中D的行列式和D的主對角線元素的變換曲線
plt.subplot(1,2,1)
for i in range(X.shape[1]):plt.plot(range(len(list1[i])),list1[i])
plt.xlabel('迭代次數(shù)');plt.ylabel('各分量的方差')
plt.subplot(1,2,2)
plt.plot(range(len(list2)),list2,'r-')
plt.xlabel('迭代次數(shù)');plt.ylabel('D的行列式')
plt.show()
# 對載荷矩陣A施以二維正交旋轉(zhuǎn)(后面的代碼與主成分法后面的代碼相同)
# 定義載荷間的總方差函數(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
#計算并修正載荷矩陣的旋轉(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ù)達到最大的旋轉(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 標準化最大方差法')
# 輸出載荷矩陣的旋轉(zhuǎn)解
print('-'*54);print('載荷矩陣的旋轉(zhuǎn)解')
print('指標|組件',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('指標'+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('總貢獻',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('指標|組件',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('指標'+str(i+1),end='\t')
? ?for j in range(S.shape[1]):print('%.3f'%S[i,j],end='\t')
? ?print()
#計算各樣本的綜合得分并排序
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,'指標'+str(i+1))
plt.show()