【机器学习】PCA原理

声明：写此博客是为了日后方便复习，加深印象，交流学习。如有错误还请指教。
更多学习可以参考：刘建平Pinard：主成分分析（PCA）原理总结，

PCA算法

** PCA算法描述**
主成分分析（Principal components analysis，以下简称PCA）是最常用的一种降维算法。
PCA是一种分析、简化数据集的技术。是数据维数压缩，尽可能降低原数据的维度（复杂度），**损失少量信息。可以削减回归分析或者聚类分析中的特征数量。
** 举例
简单的降维情况就是从二维降到一维，数据如下图。我们希望找到某一个维度方向，它可以代表这两个维度的数据。图中列了两个向量方向${\mu }_1$,${\mu }_2$。那么哪个向量可以更好的代表原始数据集呢？从直观上也可以看出，${\mu }_1$比${\mu }_2$要好。
为什么${\mu }_1$比${\mu }_2$要好？可以有两种解释，第一种解释是样本点到这个直线的距离足够近，第二种解释是样本点在这个直线上的投影能尽可能的分开。假如我们把n’从1维推广到任意维，则我们的希望降维的标准为：样本点到这个超平面的距离足够近,或者说样本点在这个超平面上的投影能尽可能的分开。
向量的表示
假设有A,B两个向量$A=(a_1，a_2...a_n)，B=(b_1，b_2...b_n)$
内积：$(a_1，a_2...a_n{)}^T\ast (b_1，b_2...b_n{)}^T=a_1{b}_1+a_2{b}_2...a_n{b}_n$;
解释：A*B=|A||B|cos(a)（设向量B的模为1，则A和B的内积值等于A向B所在直线投影的矢量长度。）

向量可以表示为（x，y），实际上表示线性组合$x(1,0{)}^T+y(0,1{)}^T$,(1,0),(0,1)称为二维空间的一组基。
基变换：
$P=(p_1，p_2...p_n)$代表一组基，$(a_1，a_2...a_m)$表示一条数据。两个矩阵相乘的意义是将右边矩阵中的每一列列向量变换到左边矩阵中每一行行向量为基所表示的空间中去。
协方差矩阵：
方向：如何选择基才能尽量保留最多的原始信息？一种直观的说法是：希望投影后的投影值尽可能的分散。
方差：$Var(a)=\frac{1}{m}{\sum }_{i=1}^m(a_i-\mu {)}^2$。
协方差(假设均值为0）：$Cov(a,b)=\frac{1}{m}{\sum }_{i=1}^m{a}_i{b}_i$
1.寻找一个一维基，使得所有数据变换为这个基上的坐标表示后，方差最大。
优化目标：
将一组N维向量降为K维（K大于0，小于N），目标是选择K个单位正交基，使原始数据变换到这组基上后；各字段两两之间协方差为0，字段的方差尽可能的大。
协方差矩阵：$X=\left(\begin{array}{cccc}{a}_1& a_2& ...& a_m\\ b_1& b_2& ...& b_m\end{array}\right)$

矩阵对角线上的两个元素分别是两个字段的方差，而其它元素是a和b的协方差。可以看出矩阵是一个实对称阵，可以找到一个n个单位正交特征向量$E=(e_1，e_2...e_n)$
进行对角化。
$E^{T}CE=A=\begin{array}{cccc}{\lambda }_1& & & \\ & {\lambda }_2& & \\ & & ...& \\ & & & {\lambda }_n\end{array}$
根据特征值的从大到小，将特征向量从上到下排列，则用K行组成的矩阵乘以原始数据矩阵X，就可以得到我们需要的降维后的矩阵Y。
PCA算法流程
(1）对所有的样本进行中心化：$x^{(i)}=x^{(i)}-\frac{1}{m}{\sum }_{j=1}^m{x}^{(j)}$;
(2)计算协方差矩阵$X{X}^T$;
(3)对矩阵$X{X}^T$进行特征值分解；
(4)取出最大的n’个特征值对应的特征向量$(w_1，w_2...w_{n^{\prime }})$ ,将所有的特征向量标准化后，组成特征向量矩阵W。
(5)对样本集中的每一个样本$x^{(i)}$,转换为新的样本$z^{(i)}=W{x}^{(i)}$;
(6)得到输出样本集$D=(z^{(1)},z^{(2)},...z^{(m)})$;
PCA实例

代码
时间：2020-3-3
实验采用的是iris数据集。

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()

实验数据拥有四种属性，三种类别，使用以下代码将属性类别详细内容可视化；

X = df.ix[:,0:4].values
y = df.ix[:,4].values
from matplotlib import pyplot as plt
import math

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):
    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()

对数据进行标准化处理

from sklearn.preprocessing import StandardScaler
X_std = StandardScaler().fit_transform(X)
#计算协方差矩阵
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)
#计算特征值、特征向量，并对特征值做归一化处理
cov_mat = np.cov(X_std.T)
eig_vals, eig_vecs = np.linalg.eig(cov_mat)
print('Eigenvectors \n%s' %eig_vecs)
print('\nEigenvalues \n%s' %eig_vals)
# Make a list of (eigenvalue, eigenvector) tuples
eig_pairs = [(np.abs(eig_vals[i]), eig_vecs[:,i]) for i in range(len(eig_vals))]
print (eig_pairs)
print ('----------')
# Sort the (eigenvalue, eigenvector) tuples from high to low
eig_pairs.sort(key=lambda x: x[0], reverse=True)

# Visually confirm that the list is correctly sorted by decreasing eigenvalues
print('Eigenvalues in descending order:')
for i in eig_pairs:
    print(i[0])
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    
a = np.array([1,2,3,4])
print (a)
print ('-----------')
print (np.cumsum(a))
#特征向量重要程度用柱状图表示出来
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()

#取出前两个特征值，并进行组合成4x2矩阵
matrix_w = np.hstack((eig_pairs[0][1].reshape(4,1),
                      eig_pairs[1][1].reshape(4,1)))
print('Matrix W:\n', matrix_w)
Y = X_std.dot(matrix_w)
Y
#查看原始属性分布（未作PCA之前）
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()