# 二維正交旋轉(zhuǎn)

import numpy as np
import matplotlib.pyplot as plt
import matplotlib;matplotlib.rc("font",family='Microsoft YaHei')
#載荷矩陣的初始解
A=np.mat([[0.796,0.424],[0.731,0.610],[0.964,-0.235],[0.953,-0.285],
? ? ? ? ?[0.940,-0.323],[0.919,-0.379],[0.793,0.284],[0.881,0.160]])
#輸出旋轉(zhuǎn)前的載荷矩陣、共同度和Fj的方差貢獻(xiàn)
print('--------旋轉(zhuǎn)前的載荷矩陣--------')
print('分量|組件 ?組件1 ? 組件2 ? 共同度')
for i in range(A.shape[0]):
? ?print('分量'+str(i+1),end='\t')
? ?for j in range(A.shape[1]):print('%.3f'%A[i,j],end='\t')
? ?print('%.3f'%sum(A[i,j]**2 for j in range(A.shape[1])))
print('總貢獻(xiàn)',end='\t')
for j in range(A.shape[1]):print('%.3f'%sum(A[i,j]**2 for i in range(A.shape[0])),end='\t')
print('%.3f'%sum(sum(A[i,j]**2 for i in range(A.shape[0])) for j in range(A.shape[1])),end='\t')
print()
#定義載荷間的總方差函數(shù)
def ff(θ):
? ?Q=np.mat([[np.cos(θ),-np.sin(θ)],
? ? ? ? ? ? ?[np.sin(θ),np.cos(θ)]])
? ?h=[sum(A[i,j]**2 for j in range(A.shape[1])) for i in range(A.shape[0])]
? ?B=A*Q
? ?junzhi=[sum(B[i,j]**2/h[i] for i in range(A.shape[0]))/A.shape[0] for j in range(A.shape[1])]
? ?ff=sum(sum((B[i,j]**2/h[i]-junzhi[j])**2 for i in range(A.shape[0])) for j in range(A.shape[1]))
? ?return ff
#繪制載荷間的總方差函數(shù)隨θ的變化圖像并找到最大值對應(yīng)的θ
x=np.arange(0,np.pi/2,0.01);y=[ff(i) for i in x]
θ=x[y.index(max(y))];plt.plot(x,y,'r-')
plt.plot([θ,θ],[min(y),max(y)],'b*--')
plt.plot([0,θ],[max(y),max(y)],'b*--')
plt.xlabel('旋轉(zhuǎn)角(弧度制)');plt.ylabel('載荷間的總方差')
#計算使載荷間的總方差函數(shù)達(dá)到最大的旋轉(zhuǎn)矩陣Q
Q=np.mat([[np.cos(θ),-np.sin(θ)],[np.sin(θ),np.cos(θ)]])
#計算旋轉(zhuǎn)后的載荷矩陣
B=A*Q
#輸出旋轉(zhuǎn)后的載荷矩陣、共同度和Fj的方差貢獻(xiàn)
print('--------旋轉(zhuǎn)后的載荷矩陣--------')
print('分量|組件 ?組件1 ? 組件2 ? 共同度')
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,j]**2 for j in range(B.shape[1])))
print('總貢獻(xiàn)',end='\t')
for j in range(B.shape[1]):print('%.3f'%sum(B[i,j]**2 for i in range(B.shape[0])),end='\t')
print('%.3f'%sum(sum(B[i,j]**2 for i in range(B.shape[0])) for j in range(B.shape[1])))
print('-'*29)
#輸出旋轉(zhuǎn)角度θ和載荷間的總方差最大值
print('旋轉(zhuǎn)角度θ='+str('%.0f'%(θ*180/np.pi))+'°')
print('載荷間的總方差='+str('%.3f'%ff(θ)))
#輸出使載荷間的總方差函數(shù)達(dá)到最大的旋轉(zhuǎn)矩陣Q
print('-----------變換矩陣-----------')
for i in range(Q.shape[0]):
? ?for j in range(Q.shape[1]):print('%.3f'%Q[i,j],end='\t')
? ?print()
print('-'*29);plt.show()