【机器学习】数据降维—主成分分析（PCA）

本文代码推荐使用Jupyter notebook跑，这样得到的结果更为直观。

主成分分析（PCA）

特征抽取通常用于提高计算效率，降低维度灾难。

主成分分析（Principe component analysis，PCA）：

       是一种广泛应用于不同领域的无监督线性数据转换技术，作用是降维。

 

常用领域：

股票交易市场数据的探索性分析和信号去噪、生物信息学领域的基因组和基因表达水平数据分析

 

PCA可以基于特征之间的关系识别出数据内在模式

PCA的目标：在高维数据中找到最大方差的方向，并将数据映射到一个维度不大于原始数据的新的子空间上。

 

PCA图：

新特征的坐标是相互正交为约束条件，子空间上的正交的坐标轴（PC）为方差最大方向。

x1和x2为原始特征坐标轴，pc1和pc2为主成分。

构建一个d x k维的转换矩阵W，将一个样本向量x映射到一个新的k维特征子空间上，此空间维度小于原始的d维特征空间。

完成从原始的d维数据到新的k维子空间转换后，第一主成分的方差应该最大，由于各主成分是正交的，后续主成分也可能具备更大方差。

 

主成分方向对数据值的范围高度敏感，如果特征值不同维度应该先对特征标准化处理，让各特征具有相同的重要性。

PCA算法流程：

 1、       对原始d维数据集做标准化

 2、       构造样本的协方差矩阵

 3、       计算协方差矩阵的特征值和相应的特征向量

 4、       选择与前k个最大特征值对应的特征向量，其中k为新特征空间维度(k<=d)

 5、       通过前k个特征向量构建映射矩阵W

 6、       通过映射矩阵W将d维的输入数据集X转换到新的k维特征子空间

第一步，加载数据集，标准化数据集。

# 加载葡萄酒数据集
import pandas as pd
df_wine = pd.read_csv('https://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data', header=None)
# 将数据分成70%的培训和30%的测试子集。
from sklearn.cross_validation import train_test_split
X, y = df_wine.iloc[:, 1:].values, df_wine.iloc[:, 0].values
X_train, X_test, y_train, y_test = \
        train_test_split(X, y, test_size=0.3, random_state=0)
# 使用单位方差标准化
from sklearn.preprocessing import StandardScaler
sc = StandardScaler()
X_train_std = sc.fit_transform(X_train)
X_test_std = sc.fit_transform(X_test)

第二步，构造协方差矩阵，dxd维协方差矩阵是沿对角线对称，d为数据集的维度，矩阵存储了不同特征之间的协方差。

协方差公式：

μj和μk分别为特征j和k的均值。

标准化后均值为0。

两个特征之间的协方差为正，则两个特征同时递减。

协方差为负，则两个特征反向移动

协方差矩阵的特征向量代表主成分，对应的特征值大小决定特征向量的重要性。

# 协方差矩阵的特征分解，计算数据集协方差矩阵的特征对。
import numpy as np
cov_mat = np.cov(X_train_std.T)
eigen_vals, eigen_vecs = np.linalg.eig(cov_mat)
print('\nEigenvalues \n%s' % eigen_vals)

cov函数得到标准化处理的训练集协方差矩阵

eig函数进行特征分解，得到特征向量及其对应的特征值。

只选择包含最多信息的特征向量组成子集。

特征值决定了特征向量的重要性，需要将特征值按降序排列，取排序在前k的特征值对应的特征向量。

 

绘制特征值的方差贡献率图像，某个特征值与所有特征值和的比较

# 使用NumPy的cumsum函数，计算累计方差。
tot = sum(eigen_vals)
var_exp = [(i / tot) for i in sorted(eigen_vals, reverse=True)]
cum_var_exp = np.cumsum(var_exp)
# 通过Plo的step绘制
import matplotlib.pyplot as plt
%matplotlib inline
plt.bar(range(1, 14), var_exp, alpha=0.5, align='center',
        label='individual explained variance')
plt.step(range(1, 14), cum_var_exp, where='mid',
         label='cumulative explained variance')
plt.ylabel('Explained variance ratio')
plt.xlabel('Principal components')
plt.legend(loc='best')
plt.tight_layout()
# plt.savefig('./figures/pca1.png', dpi=300)
plt.show()

第一主成分占方差总和40%，前两个主成分占比近60%

PCA是一种无监督方法，可以忽略类标信息

随机森林通过类标信息计算节点的不纯度，方差度量的是特征值在轴线是的分布。

# 按特征值的降序排列特征对
# 列出(特征值，特征向量)元组。
eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vecs[:,i]) for i in range(len(eigen_vals))]

# 从高到低排序(特征值，特征向量)元组。
eigen_pairs.sort(reverse=True)
# 本案例只选择前60%的两个特征向量
w = np.hstack((eigen_pairs[0][1][:, np.newaxis],
               eigen_pairs[1][1][:, np.newaxis]))
print('Matrix W:\n', w)

# 通过计算矩阵点积，将整个训练集转换到包含两个主成分的子空间上。
X_train_pca = X_train_std.dot(w)

# 可视化
colors = ['r', 'b', 'g']
markers = ['s', 'x', 'o']

for l, c, m in zip(np.unique(y_train), colors, markers):
    plt.scatter(X_train_pca[y_train==l, 0], 
                X_train_pca[y_train==l, 1], 
                c=c, label=l, marker=m)

plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.legend(loc='lower left')
plt.tight_layout()
# plt.savefig('./figures/pca2.png', dpi=300)
plt.show()

# 使用Scikit-learn进行主成分分析
from sklearn.decomposition import PCA

pca = PCA()
X_train_pca = pca.fit_transform(X_train_std)
# pca.explained_variance_ratio_
plt.bar(range(1, 14), pca.explained_variance_ratio_, alpha=0.5, align='center')
plt.step(range(1, 14), np.cumsum(pca.explained_variance_ratio_), where='mid')
plt.ylabel('Explained variance ratio')
plt.xlabel('Principal components')
plt.show()

pca = PCA(n_components=2)
X_train_pca = pca.fit_transform(X_train_std)
X_test_pca = pca.transform(X_test_std)
plt.scatter(X_train_pca[:,0], X_train_pca[:,1])
plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.show()

# 使用plot_decision_regions函数进行可视化决策区域
from matplotlib.colors import ListedColormap

def plot_decision_regions(X, y, classifier, resolution=0.02):

    #设置标记生成器和颜色映射。
    markers = ('s', 'x', 'o', '^', 'v')
    colors = ('red', 'blue', 'lightgreen', 'gray', 'cyan')
    cmap = ListedColormap(colors[:len(np.unique(y))])

    # plot 决定表面
    x1_min, x1_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    x2_min, x2_max = X[:, 1].min() - 1, X[:, 1].max() + 1
    xx1, xx2 = np.meshgrid(np.arange(x1_min, x1_max, resolution),
                         np.arange(x2_min, x2_max, resolution))
    Z = classifier.predict(np.array([xx1.ravel(), xx2.ravel()]).T)
    Z = Z.reshape(xx1.shape)
    plt.contourf(xx1, xx2, Z, alpha=0.4, cmap=cmap)
    plt.xlim(xx1.min(), xx1.max())
    plt.ylim(xx2.min(), xx2.max())

    # plot 类样本
    for idx, cl in enumerate(np.unique(y)):
        plt.scatter(x=X[y == cl, 0], y=X[y == cl, 1],
                    alpha=0.8, c=cmap(idx),
                    marker=markers[idx], label=cl)

# 使用前两个主要组件训练逻辑回归分类器
from sklearn.linear_model import LogisticRegression

lr = LogisticRegression()
lr = lr.fit(X_train_pca, y_train)
plot_decision_regions(X_train_pca, y_train, classifier=lr)
plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.legend(loc='lower left')
plt.tight_layout()
# plt.savefig('./figures/pca3.png', dpi=300)
plt.show()

SKlearn实现的PCA和之前实现的PCA是经过y轴旋转后的。

特征分析方法：特征向量可以为正或者为负

有时候需要乘上-1在实现图像的镜像。

# 绘制逻辑回归在转换后的测试数据上得到的决策区域
plot_decision_regions(X_test_pca, y_test, classifier=lr)
plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.legend(loc='lower left')
plt.tight_layout()
# plt.savefig('./figures/pca4.png', dpi=300)
plt.show()

# 获取相应的方差贡献率
pca = PCA(n_components=None)
X_train_pca = pca.fit_transform(X_train_std)
pca.explained_variance_ratio_