经验正交函数 (EOF) / 主成分 (PCA) 分解及Python实现

前言

EOF 经验正交函数 分解原理及用法:
https://wenku.baidu.com/view/93bc8a9b680203d8ce2f2468.html

方法实现过程

EOF

经验正交函数 (EOF) / 主成分 (PCA) 分解及Python实现

PCA

经验正交函数 (EOF) / 主成分 (PCA) 分解及Python实现

Python实现

import numpy as np
from scipy import linalg

# 生成数据并打乱
np.random.seed(9)
data = np.random.randint(low=1, high=50, size=600).reshape(20, 30) * 10

# 时间距平
mean = np.repeat(data.mean(axis=1), 30).reshape(20, 30)
data = data - mean

# 交叉积
C = np.dot(data, data.T) / 30
print(C)

# 求特征值和特征向量 (空间模态) eigh求实对称矩阵的特征值 向量
la, v = linalg.eigh(C) # 空间模态矩阵 行代表着某个站点, 列代表第几空间模态

# 根据特征值对特征向量升序排序
index = np.argsort(la)[::-1]
v = v[:, index]

# 计算时间系数矩阵 第i行代表i-th时间模态 第k列代表k时刻时间系数  
# 理论上这里就是pca, np.dot(v.T[0:2, :], data) 取前三个主成分
# 如果用sklearn 的PCA API 要先将 data 转成特征成列分布的,而不是成行分布,结果差不多,
# 可能符号有区别 此外,sklearn的PCA考虑了计算成本是使用SVD分解实现的
# from sklearn.decomposition import PCA
# pca_sk = PCA(n_components=2).fit_transform(data.T)

PC = np.dot(v.T, data)

# 证明PC具有正交性
pcpcpc = np.around(np.dot(PC, PC.T) / 30, 5)
vvv = np.around(np.dot(v, v.T), 5)
vvv_ = np.around(np.dot(v.T, v), 5)

# 还原数据
data_reverse = np.dot(v, PC)
上一篇:PCA方法校正群体结构,GWAS该用多少个主成分?


下一篇:day2021_9_26