python实现reof

Alexander

花了几天时间研究了下python版本的reof，抛砖引玉，时间仓促还没用真实气象数据验证，鉴于本人水平有限可能许多地方理解有出入，重在理解。

1. import numpy as np
   
2. from tqdm import tqdm
   
3. 
   
4. #生成一个玩具数据集，840个站点1200个时次的数据
   
5. x  = np.random.randint(1,20,size=(840,1200))
   
6. def standardization(data):
   
7.     """
   
8.     标准化资料阵
   
9.     """
   
10.     mu = np.mean(data, axis=1,keepdims = True)
    
11.     #print(mu.shape)
    
12.     m = np.tile(mu,(1,data.shape[1]))
    
13.     #print(m.shape)
    
14.     sigma = np.std(data, axis=1,keepdims = True)
    
15.     s = np.tile(sigma,(1,data.shape[1]))
    
16.     fengzi = (data - m)
    
17.     fengmu = s
    
18.     return fengzi / fengmu
    
19. 
    
20. #对资料阵标准化
    
21. data = standardization(x)
    
22. #计算协方差阵,也时相关阵
    
23. covx = 1.0/x.shape[1]*np.matmul(x,x.T)
    
24. #covx = np.cov(x,rowvar=1)
    
25. print(covx.shape)
    
26. #计算特征值和特征向量
    
27. eigenvalues,eigenvecter =  np.linalg.eig(covx)
    
28. #print(np.matmul(eigenvecter.T,eigenvecter))
    
29. print(eigenvecter[:,0].var())
    
30. #根据特征值（方差贡献排序）
    
31. idx = np.argsort(eigenvalues)[::-1]
    
32. eigenvecter = eigenvecter [:,idx]
    
33. -----------------------------------------------------------------
    
34. #根据累计方差贡献大于85%来确定待旋转的特征向量
    
35. totalvar = np.sum(eigenvalues)
    
36. cumsum = 0
    
37. #eigvalue_sort  = sorted(eigenvalues, reverse=True)
    
38. eigvalue_sort = eigenvalues[idx]
    
39. #print(eigvalue_sort)
    
40. for number,i in enumerate(eigvalue_sort):
    
41.     cumsum += i/totalvar
    
42.     if cumsum > 0.85:
    
43.         break
    
44. #print(number)
    
45. -------------------------------------------------
    
46. rotate_eigvector = eigenvecter[:,:number+1]
    
47. eigvalue_sort = eigvalue_sort[:number+1]
    
48. 
    
49. #极大方差旋转，刚性旋转，轴是正交的
    
50. #1.使用wiki里的迭代解法，但无论用哪个版本都需要将特征向量转变成载荷矩阵
    
51. loadings = np.dot(rotate_eigvector, np.diag(np.sqrt(eigvalue_sort)))
    
52. #loadings = eVec_corr3*np.sqrt(eVal_corr3)
    
53. n_rows, n_cols = loadings.shape
    
54. #对loadings标准化，这里是横向给每个元素除共同度的开根也就是h
    
55. normalized_mtx = np.apply_along_axis(lambda x: np.sqrt(np.sum(x**2)), 1, loadings.copy())
    
56. normlize_loadings = (loadings.T / normalized_mtx).T
    
57. #print(normlize_loadings[:,0].mean())
    
58. #print((loadings[:,0]*loadings[:,1]).mean())
    
59. #assert loadings == n]ormlize_loadings
    
60. #初始化旋转矩阵，为1的对角矩阵，意味着第一次旋转相当于原地不动
    
61. rotation_mtx = np.eye(n_cols)
    
62. d = 0
    
63. for i in range(1000):
    
64.   print('process-->',i)
    
65.   old_d = d
    
66.   #basis应该是旋转后的载荷矩阵
    
67.   basis = np.dot(loadings, rotation_mtx)
    
68.   """
    
69.   下面这个地方最难理解并且是最核心的地方，因为这个公式实在太抽象了我试着猜一下大概是怎么回事，这里参考我下面贴的svd旋转矩阵的链接
    
70.   首先套公式argmax（x.T*y）,如果载荷矩阵的转置是x的话，那么很容易猜到这个后面一长串的公式就是这个y，根据帖子中x,y都是多维空间中的点
    
71.   群，x是原来的载荷矩阵，那么y就应该是更具最大方差原则处理后的载荷矩阵，并且我猜这里没有正交的约束，因为svd分解在后面，它的跳出循环的
    
72.   条件应该就是不满足正交矩阵了。这个basis3次方困惑了我很久，根据
    
73.   黄嘉佑老师书中的公式，可以确定的最大化方差的计算中每一个因子（特征向量）（列）中每个值都是平方的，在求方差应该
    
74.   是有四次方项的，但这个y肯定包含了方差最大化的条件，所以我从公式的右半侧开始看np.diag(np.diag(np.dot(basis.T, basis)这个
    
75.   最好理解它的作用是先求每个因子的方差然后对角化在除于因子（特征向量）的长度，相当于因子方差平均到每个feature（站点）上的大小
    
76.   接下来就纯属于我个人理解（猜测），因子旋转的目的是让因子中的某几个feature（站点或格点）具有高载荷，它不同于pca（或eof）只需要让
    
77.   方差集中在前几个因子上就行了，那么旋转的目的就是让大的更大小的更小，所以根据这个思路这个basis**3应该拆开来看，这里的三次方是逐元素
    
78.   的幂，将其拆成basis*（basis**2）这里的*是哈达玛逐元素乘积，那么这个公式应该可以合并成basis*（basis**2 - 每一个因子元素的平均方差）
    
79.   那么根据这个思路（如果正确的话）就比较好理解了，大于平均方差的项朝着更大的方向转，小于的则更小这也符合旋转的原则，最后在这个方向上根
    
80.   据svd求出正交的部分（旋转矩阵）直到按照这个方向旋转求不出正交部分为止，这根家园fortran版本里的解法应该是相反的（我只是简单的看了
    
81.   下），fortran的版本是先固定正交的条件，然后两两排列组合旋转，所以也非常好理解。
    
82.   
    
83.   """
    
84.   transformed = np.dot(loadings.T, basis**3 - (1.0 / n_rows) * np.dot(basis, np.diag(np.diag(np.dot(basis.T, basis)))))
    
85.   U, S, V = np.linalg.svd(transformed)
    
86.   #下面参考svd求旋转矩阵的连接wiki里也有简要的推导
    
87.   rotation_mtx = np.dot(U, V)
    
88.   d = np.sum(S)
    
89.   #这里的意思应该是保证对角线上奇异值的和保持不变从而保证旋转的矩阵是正交的
    
90.   if old_d != 0 and d / old_d < 1 + 1e-05 :
    
91.       break
    
92. new_loadings = np.dot(loadings, rotation_mtx)
    
93. print(new_loadings)
    
94. #通过旋转后的载荷矩阵计算新的方差贡献，载荷矩阵对应的因子(模态)的总方差跟未旋转的总方差相同
    
95. new_variance = np.apply_along_axis(lambda x: np.sum(x**2), 0, new_loadings)
    
96. old_variance = np.apply_along_axis(lambda x: np.sum(x**2), 0, loadings)
    
97. new_variance = new_variance/new_variance.sum()
    
98. old_variance = old_variance/old_variance.sum()
    
99. 
    
100. 
     
101. 
     
102. #参考链接
     
103. #wiki源码链接：https://en.wikipedia.org/wiki/Talk%3aVarimax_rotation
     
104. #因子分析：https://blog.csdn.net/mengjizhiyou/article/details/83241469
     
105. #varimax ：https://github.com/odeviron/varimax/blob/master/varimax_lib.py
     
106. #https://github.com/trhgu/Varimax_Rotation_and_Thereafter_module/blob/master/Varimax_Rotation_and_Thereafter.ipynb
     
107. #https://github.com/EducationalTestingService/factor_analyzer/blob/master/factor_analyzer/factor_analyzer.py
     
108. #https://github.com/EducationalTestingService/factor_analyzer/blob/master/factor_analyzer/rotator.py
     
109. #svd求解旋转矩阵
     
110. #https://blog.csdn.net/dfdfdsfdfdfdf/article/details/53213240

复制代码

Alexander

希望能跟大家交流一下，我给程序里写了详细的注释，希望有高手指点一下我对transformed计算的那一大段的理解是否正确{:eb302:}{:eb302:}{:eb302:}{:eb302:}

lancerv8

感谢分享，可以尝试用python进行一次reof分析了。

heiheihie

多谢楼主分享呀

灭火器

之前好像在关于因子分析的心理学的论文里看到过推导，似乎是直接当最优化问题来求，实在是看不懂，只能调包来用了。

Alexander

灭火器 发表于 2020-8-31 10:11
之前好像在关于因子分析的心理学的论文里看到过推导，似乎是直接当最优化问题来求，实在是看不懂，只能调包 ...

fortran版本的reof也是解最优化，包括各种气象统计书上的算法都是每次只旋转两个，但维基上的算法是直接对载荷矩阵操作，就很迷