以下是使用PCA算法处理实际问题的例子，同样使用鸢尾花数据集，目的依旧是完成降维任务
基本的流程如下：
1.数据预处理，只有数值数据才可以进行PCA降维
2.计算样本数据的协方差方阵
3.求解协方差矩阵的特征值和特征向量
4.将特征值按照从大到小的顺序排列，选择其中较大的K个，然后将其对应的K个特征向量组成投影矩阵
5.将样本点投影计算，完成PCA降维任务
1、导入数据
import numpy as np import pandas as pd # 读取数据集 df = pd.read_csv('iris.data') # 原始数据没有给定列名的时候需要我们自己加上 df.columns=['sepal_len', 'sepal_wid', 'petal_len', 'petal_wid', 'class'] df.head()
2、展示数据特征
# 把数据分成特征和标签 X = df.iloc[:,0:4].values y = df.iloc[:,4].values from matplotlib import pyplot as plt # 展示我们标签用的 label_dict = {1: 'Iris-Setosa', 2: 'Iris-Versicolor', 3: 'Iris-Virgnica'} # 展示特征用的 feature_dict = {0: 'sepal length [cm]', 1: 'sepal width [cm]', 2: 'petal length [cm]', 3: 'petal width [cm]'} # 指定绘图区域大小 plt.figure(figsize=(8, 6)) for cnt in range(4): # 这里用子图来呈现4个特征 plt.subplot(2, 2, cnt+1) for lab in ('Iris-setosa', 'Iris-versicolor', 'Iris-virginica'): plt.hist(X[y==lab, cnt], label=lab, bins=10, alpha=0.3,) plt.xlabel(feature_dict[cnt]) plt.legend(loc='upper right', fancybox=True, fontsize=8) plt.tight_layout() plt.show()
可以看见，有些特征区别能力较强，能把3种花各自呈现出来；有的特征区别能力较弱，部分特征数据样本混杂在一起。
3、数据标准化
一般情况下，在进行训练前，数据经常需要进行标准化处理。
from sklearn.preprocessing import StandardScaler X_std = StandardScaler().fit_transform(X)
4、计算协方差矩阵
mean_vec = np.mean(X_std, axis=0) cov_mat = (X_std - mean_vec).T.dot((X_std - mean_vec)) / (X_std.shape[0]-1) print('协方差矩阵 \n%s' %cov_mat) # 利用numpy也可以 # print('NumPy 计算协方差矩阵: \n%s' %np.cov(X_std.T))
协方差矩阵 [[ 1.00675676 -0.10448539 0.87716999 0.82249094] [-0.10448539 1.00675676 -0.41802325 -0.35310295] [ 0.87716999 -0.41802325 1.00675676 0.96881642] [ 0.82249094 -0.35310295 0.96881642 1.00675676]]
5、求特征值与特征向量
cov_mat = np.cov(X_std.T) eig_vals, eig_vecs = np.linalg.eig(cov_mat) print('特征向量 \n%s' %eig_vecs) print('\n特征值 \n%s' %eig_vals)
特征向量 [[ 0.52308496 -0.36956962 -0.72154279 0.26301409] [-0.25956935 -0.92681168 0.2411952 -0.12437342] [ 0.58184289 -0.01912775 0.13962963 -0.80099722] [ 0.56609604 -0.06381646 0.63380158 0.52321917]] 特征值 [2.92442837 0.93215233 0.14946373 0.02098259]
6、按照特征值大小进行排序
# 把特征值和特征向量对应起来 eig_pairs = [(np.abs(eig_vals[i]), eig_vecs[:,i]) for i in range(len(eig_vals))] print (eig_pairs) print ('----------') # 把它们按照特征值大小进行排序 eig_pairs.sort(key=lambda x: x[0], reverse=True) # 打印排序结果 print('特征值又大到小排序结果:') for i in eig_pairs: print(i[0])
[(2.9244283691111126, array([ 0.52308496, -0.25956935, 0.58184289, 0.56609604])), (0.9321523302535072, array([-0.36956962, -0.92681168, -0.01912775, -0.06381646])), (0.14946373489813383, array([-0.72154279, 0.2411952 , 0.13962963, 0.63380158])), (0.020982592764270565, array([ 0.26301409, -0.12437342, -0.80099722, 0.52321917]))] ---------- 特征值又大到小排序结果: 2.9244283691111126 0.9321523302535072 0.14946373489813383 0.020982592764270565
7、计算累积结果
将特征向量累加起来，超过一定百分比时，就可以选择其为降维后的维度大小
# 计算累加结果 tot = sum(eig_vals) var_exp = [(i / tot)*100 for i in sorted(eig_vals, reverse=True)] print (var_exp) cum_var_exp = np.cumsum(var_exp) cum_var_exp
[72.62003332692029, 23.147406858644153, 3.711515564584534, 0.5210442498510144] array([ 72.62003333, 95.76744019, 99.47895575, 100. ])
可以发现，使用前两个特征值时，其对应的累积贡献率已经超过了95%，所以选择降到了二维。
# cumsum的用法例子 a = np.array([1,2,3,4]) print (a) print ('-----------') print (np.cumsum(a))
[1 2 3 4] ----------- [ 1 3 6 10]
画图可以更直接的展示
plt.figure(figsize=(6, 4)) plt.bar(range(4), var_exp, alpha=0.5, align='center', label='individual explained variance') plt.step(range(4), 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.show()
8、完成PCA降维
将前两个特征向量组合起来完成降维操作
matrix_w = np.hstack((eig_pairs[0][1].reshape(4,1), eig_pairs[1][1].reshape(4,1))) print('Matrix W:\n', matrix_w)
Matrix W: [[ 0.52308496 -0.36956962] [-0.25956935 -0.92681168] [ 0.58184289 -0.01912775] [ 0.56609604 -0.06381646]]
Y = X_std.dot(matrix_w) print("X.shape : ",X.shape) print("Y.shape : ",Y.shape)
X.shape : (149, 4) Y.shape : (149, 2)
可以看见将原来的数据从4维降到2维
9、可视化对比降维前后数据的分布
由于数据具有4个特征，无法在平面图中显示，因此只使用两维特征显示数据
plt.figure(figsize=(6, 4)) for lab, col in zip(('Iris-setosa', 'Iris-versicolor', 'Iris-virginica'), ('blue', 'red', 'green')): plt.scatter(X[y==lab, 0], X[y==lab, 1], label=lab, c=col) plt.xlabel('sepal_len') plt.ylabel('sepal_wid') plt.legend(loc='best') plt.tight_layout() plt.show()
降维后的结果
plt.figure(figsize=(6, 4)) for lab, col in zip(('Iris-setosa', 'Iris-versicolor', 'Iris-virginica'), ('blue', 'red', 'green')): plt.scatter(Y[y==lab, 0], Y[y==lab, 1], label=lab, c=col) plt.xlabel('Principal Component 1') plt.ylabel('Principal Component 2') plt.legend(loc='lower center') plt.tight_layout() plt.show()
