主成分分析（PCA）
概述
- 主成分分析的基本原理
- 主成分分析的算法实现
- PCA 和 Kernel PCA的比较
- 基于上证50指数的主成分分析
1、主成分分析的基本原理
问题的提出：为了全面系统的分析和研究问题，必须考虑许多指标，这些指标能从不同的侧面反映所研究的对象的特征，但指标过多，会增加分析的复杂性，也有很多是相关的，原始变量能不能减少为有代表性的少数几个综合变量，用它们来代表原来的指标？
主成分分析（principal component analysis，PCA）就是寻找用较少的（降维）综合变量代替原来较多的旧变量，而且使综合变量尽可能多地保留原来较多信息的方法。
在选取综合变量时，最简单的形式就是取原来变量的线性组合，适当调整组合系数，使新的变量之间相互独立且代表性最好。
综合变量或者说主成分应该怎么选取呢？
我们希望第一个线性组合即第一个综合变量记为F1，尽可能多的反映原来变量的信息。最经典的方法就是用方差来表达，即var(F1)越大，表示F1包含的信息越多。因此在所有的线性组合中所选取的F1应该是反差最大的，故称之为第一主成分。如果第一主成分不足以代表原来旧变量的信息，再考虑选取F2即第二个线性组合，成为第二主成分。F2是与F1不相关的旧变量的所有线性组合中方差最大的。F3、F4…..依次类推。
以二维空间为例：
可以看出这n个样本点无论是沿着X1轴方向或X2轴方向都具有较大的离散性，其离散程度可以分别用观测变量X1的方差和X2的方差定量地表示。显然，如果只考虑X1和X2中的任何一个，那么包含在原始数据中的信息将会有较大的损失。如果我们将X1轴和X2轴先平移，再同时旋转，得到新坐标轴Z1和Z2。Z1和Z2是两个新变量，也是变量X1和X2的线性组合。
二维平面上的各点的方差大部分都归结在Z1轴上，而Z2轴上的方差很小，这里我们可以只考虑Z1就能将原始数据的大部分信息保存。Z1，Z2除了可以对包含在X1，X2中的信息起着浓缩作用以外，还具有不相关的性质，这就使得在研究复杂的问题时避免了信息重叠所带来的虚假性。
主成分分析的计算步骤
1. 计算相关系数矩阵R；
2. 计算特征值与特征向量；
解特征方程|λE-R|=0,求出特征值λ(i)（从大到小的顺序排列）和特征向量e(i)。
（最大方差问题转化为求特征值问题）
3. 计算主成分贡献率及累计贡献率；
贡献率：
累计贡献率：
4. 计算主成分载荷，得出主成分模型。
主成分模型为：
，
2、主成分分析的算法实现
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
# step 1 随机生成6*4的整数矩阵X
np.random.seed(1234)
X = np.random.randint(0,10,size=(6,4))
# step 2 计算X相关系数矩阵R
R = np.corrcoef(X.T)
# step 3 计算特征值Lambda和特征向量e
Lambda, e = np.linalg.eig(R)
# step 4 计算主成分贡献率及累计贡献率
k = 3
eigen_pairs = [(Lambda[i], e[:, i]) for i in range(len(Lambda))]
eigen_pairs = sorted(eigen_pairs, reverse=True, key=lambda k: k[0])
# 将特征值从大到小排列，reverse倒序，k[0]表示按第一个元素的大小进行排序，也可以写成x：x[0],可以任意更改
get_we = lambda x: x / x.sum()
weights = get_we(np.stack(eigen_pairs[i][0] for i in range(len(Lambda)))) # 主成分贡献率
acc_we = weights[:k].sum() # 累计贡献率
E = np.column_stack((eigen_pairs[i][1] for i in range(k))) # 提取前k个特征值对应的特征向量
# step 4 主成分模型
X_pca = X.dot(E)3、PCA 和 Kernel PCA的比较
PCA不进行分类的动作，而只做做数据预处理，将样本变换到一个容易分类（向最大化方差的方向，principal component axes，投影）的更低维的新的特征空间中。Kernel PCA比PCA多了一步，也即先升维（RBF包括多项式核均是升高到无穷维）再进行投影的动作，因为有些非线性可分的数据集只有在升维的视角下才线性可分。
下面我们以一组非线性数据集为例。
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import make_moons
x2, y2 = make_moons(n_samples=100, random_state=123)
sc = StandardScaler()
x2_std = sc.fit_transform(x2) # 标准正态化
plt.scatter(x2_std[y2==0, 0], x2_std[y2==0, 1], color='red', marker='^', alpha=0.5)
plt.scatter(x2_std[y2==1, 0], x2_std[y2==1, 1], color='blue', marker='o', alpha=0.5)
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.show()
我们先用PCA进行处理：
from sklearn.decomposition import PCA
pca = PCA(n_components=2)
x_spca = pca.fit_transform(x2_std)
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(14,6))
ax[0].scatter(x_spca[y2==0, 0], x_spca[y2==0, 1], color='red', marker='^', alpha=0.5)
ax[0].scatter(x_spca[y2==1, 0], x_spca[y2==1, 1], color='blue', marker='o', alpha=0.5)
ax[1].scatter(x_spca[y2==0, 0], np.zeros((50,1))+0.02, color='red', marker='^', alpha=0.5)
ax[1].scatter(x_spca[y2==1, 0], np.zeros((50,1))+0.02, color='blue', marker='o', alpha=0.5)
ax[0].set_xlabel('PC1')
ax[0].set_ylabel('PC2')
ax[1].set_ylim([-1, 1])
ax[1].set_yticks([])
ax[1].set_xlabel('PC1')
plt.show()
从图中可以看出，经过主成分分析之后，数据仍旧是非线性不可分的。接下来，我们用基于核函数的主成分分析来试下看。
from sklearn.decomposition import KernelPCA
kpca = KernelPCA(n_components=2, kernel='rbf', gamma=15)
x_kpca = kpca.fit_transform(x2_std)
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(14,6))
ax[0].scatter(x_kpca[y2==0, 0], x_kpca[y2==0, 1], color='red', marker='^', alpha=0.5)
ax[0].scatter(x_kpca[y2==1, 0], x_kpca[y2==1, 1], color='blue', marker='o', alpha=0.5)
ax[1].scatter(x_kpca[y2==0, 0], np.zeros((50,1))+0.02, color='red', marker='^', alpha=0.5)
ax[1].scatter(x_kpca[y2==1, 0], np.zeros((50,1))+0.02, color='blue', marker='o', alpha=0.5)
ax[0].set_xlabel('PC1')
ax[0].set_ylabel('PC2')
ax[1].set_ylim([-1, 1])
ax[1].set_yticks([])
ax[1].set_xlabel('PC1')
plt.show()
由图可知，只需要把数据集投影到经变换后的特征PC1上就可以实现线性划分了，这个时候只需要一个特征PC1就够了。
4、基于上证50指数的主成分分析
参考我在优矿上发表的文章：基于上证50指数的主成分分析（PCA）
参考文献：
从PCA到Kernel PCA（Python）
主成分分析（PCA）和基于核函数的主成分分析（KPCA）入门
Yves Hilpisch. Python for Finance