Python中的主成分分析(PCA)

80

我有一个 (26424 x 144) 的数组,想使用Python对其执行PCA。然而,网络上没有特定的位置来解释如何完成此任务(有一些网站只是按照自己的方式执行PCA - 没有通用的方法可供查找)。任何能提供任何形式帮助的人都将是伟大的。


你的数组是否稀疏(大部分为0)?你关心前2-3个组件捕获多少方差 - 50%,90%吗? - denis
不,它不是稀疏的,我已经过滤掉了错误值。是的,我对找出需要解释> 75%和> 90%方差的主成分数量很感兴趣...但不确定如何做。有任何想法吗? - khan
查看Doug答案中排序后的evals--如果您愿意,可以在此处或新问题中发布前几个和总和。并参见维基百科PCA累积能量 - denis
可以在这里链接找到使用仅numpy和/或scipy的基本PCA方法的比较,附有timeit结果。 - djvg
11个回答

0

这段示例代码加载了日本收益率曲线,并创建了PCA组件。 然后使用PCA估计给定日期的移动,并将其与实际移动进行比较。

%matplotlib inline

import numpy as np
import scipy as sc
from scipy import stats
from IPython.display import display, HTML
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
import datetime
from datetime import timedelta

import quandl as ql

start = "2016-10-04"
end = "2019-10-04"

ql_data = ql.get("MOFJ/INTEREST_RATE_JAPAN", start_date = start, end_date = end).sort_index(ascending= False)

eigVal_, eigVec_ = np.linalg.eig(((ql_data[:300]).diff(-1)*100).cov()) # take latest 300 data-rows and normalize to bp
print('number of PCA are', len(eigVal_))

loc_ = 10
plt.plot(eigVec_[:,0], label = 'PCA1')
plt.plot(eigVec_[:,1], label = 'PCA2')
plt.plot(eigVec_[:,2], label = 'PCA3')
plt.xticks(range(len(eigVec_[:,0])), ql_data.columns)
plt.legend()
plt.show()

x = ql_data.diff(-1).iloc[loc_].values * 100 # set the differences
x_ = x[:,np.newaxis]
a1, _, _, _ = np.linalg.lstsq(eigVec_[:,0][:, np.newaxis], x_) # linear regression without intercept
a2, _, _, _ = np.linalg.lstsq(eigVec_[:,1][:, np.newaxis], x_)
a3, _, _, _ = np.linalg.lstsq(eigVec_[:,2][:, np.newaxis], x_)

pca_mv = m1 * eigVec_[:,0] + m2 * eigVec_[:,1] + m3 * eigVec_[:,2] + c1 + c2 + c3
pca_MV = a1[0][0] * eigVec_[:,0] + a2[0][0] * eigVec_[:,1] + a3[0][0] * eigVec_[:,2]
pca_mV = b1 * eigVec_[:,0] + b2 * eigVec_[:,1] + b3 * eigVec_[:,2]

display(pd.DataFrame([eigVec_[:,0], eigVec_[:,1], eigVec_[:,2], x, pca_MV]))
print('PCA1 regression is', a1, a2, a3)


plt.plot(pca_MV)
plt.title('this is with regression and no intercept')
plt.plot(ql_data.diff(-1).iloc[loc_].values * 100, )
plt.title('this is with actual moves')
plt.show()

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接