二维正交旋转

# 二维正交旋转

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]])
#输出旋转前的载荷矩阵、共同度和Fj的方差贡献
print('--------旋转前的载荷矩阵--------')
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('总贡献',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()
#定义载荷间的总方差函数
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
#绘制载荷间的总方差函数随θ的变化图像并找到最大值对应的θ
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('旋转角(弧度制)');plt.ylabel('载荷间的总方差')
#计算使载荷间的总方差函数达到最大的旋转矩阵Q
Q=np.mat([[np.cos(θ),-np.sin(θ)],[np.sin(θ),np.cos(θ)]])
#计算旋转后的载荷矩阵
B=A*Q
#输出旋转后的载荷矩阵、共同度和Fj的方差贡献
print('--------旋转后的载荷矩阵--------')
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('总贡献',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)
#输出旋转角度θ和载荷间的总方差最大值
print('旋转角度θ='+str('%.0f'%(θ*180/np.pi))+'°')
print('载荷间的总方差='+str('%.3f'%ff(θ)))
#输出使载荷间的总方差函数达到最大的旋转矩阵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()