主成分分析（PCA）详解

Posted 2023-03-30 生信小兔

主成分分析（PCA)是一种比较基础的数据降维方法，也是多元统计中的重要部分，在数据分析、机器学习等方面具有广泛应用。主成分分析目的是用较少的变量来代替原来较多的变量，并可以反映原来多个变量的大部分信息。

1.主成分分析（PCA）原理

对于一个含有n个数据，变量的个数为p的一个样本，我们可以用p维空间的n个点来表示这些数据。例如含有2个变量，3个数据（1，2），（2，2），（3，3）的样本，即可以表示为：

如果含有3个变量，就是三维空间中的散点。

通常情况下，我们在实验初会提出很多变量，并且采集这些数据，这些数据中各个变量往往会存在一定的相关性。而这些相关性便意味着可以进行数据的降维，用更少的变量来替代这些变量。事实上，在很多情况下数据的降维是十分必要的，一方面有利于问题的简化，另一方面便于计算机的计算：数据降维后变量的减少，会使计算机处理的数据大大减少，从而缩短数据处理时间。

主成分分析的直观理解，可以认为是旋转坐标轴，使得在旋转坐标轴后这些点在新的坐标系下在各个坐标轴（变量）方向投影的方差变大。其中如果在某坐标上的方差最大，那么这个坐标轴对应的这些散点的坐标就是第一主成分，其次就是第二主成分，依此类推。

 对于上图的8组数据，我们发现在x轴方向方差很大，在y轴方向方差为0，所以就可以用这些点的横坐标数据作为第一主成分，并且只选第一主成分便可以达到要求。

对于下图的情况，我们发现这些数据都几乎排列在一条直线上，并且在x轴方向和y轴方向的方差都比较大。但是如果把坐标轴旋转一定角度，使得这些数据在某个坐标轴的投影的方差比较大，便可以用新坐标系下方差较大的一个坐标轴坐标作为主成分。

 对于左图，数据为（1，2）、（2，4）......旋转坐标轴后，坐标为（,0)、（,0）......这样主成分就是新坐标系下变量x的数值：。

对于大多数情况，数据各个变量基本服从正态分布，所以变量为2的数据散点分布大致为一个椭圆，变量为3的散点分布大致为一个椭球，p个变量的数据大致分布在一个超椭圆。而通过旋转坐标系，使得超椭圆的长轴落在一个坐标轴上，其次超椭圆另一个轴也尽量落在坐标轴上。这样各个新的坐标轴上的坐标值便是相应的主成分。

 例如，对于上图的数据，在x轴和y轴的方差都很大，所以可以旋转坐标系，使得椭圆两个轴尽量落在坐标轴上。

 这样，我们便以散点在新坐标系下的x坐标作为第一主成分（因为x方向方差最大），y轴的坐标为第二主成分。

主成分分析的理论推导较为复杂，需要借助投影寻踪，构造目标函数等方法来推导，在多元统计的相关书籍中都有详细讲解。但是其结论却是十分简洁。所以，如果只是需要实际应用，了解主成分分析的基本原理与实现方法便足够了。

2.主成分分析（PCA）步骤

对于含有p个变量、n个数据的样本，可以计算出其协差阵。根据该协差阵可以求出p个特征值并从大到小排序得与相应的p个特征向量。则第i主成分为：

原样本中含有p个变量，通过主成分分析后变量个数会减少很多。选取主成分的个数需要依据主成分贡献率与累计贡献率。第k个主成分的贡献率为：

通常情况下，该主成分的贡献率越大，说明保存的原有数据的信息越多。样本前m个主成分的累计贡献率为：

通常情况下，如果累计贡献率达到80%以上，便可以认为选取前m个主成分能够很好地保留原来样本的信息。累计贡献率是判断选取主成分个数的标准，也反映了这些主成分对原有信息的保留情况。

事实上，由于各个变量单位的不同，利用样本协差阵进行主成分分析往往不够准确，所以可以用样本的相关阵来替代协差阵，之后的计算过程同上。

在实际应用中，从p个变量中提取m个主成分，需要对这新的m个变量做合理的解释。而进行解释往往需要根据具体问题来判断，并结合该主成分中保存原各个变量的比重，即求解该主成分与原始各个变量的因子载荷量。

其中表示第k主成分在原第i变量的因子载荷量（它表示变量Xi对主成分Yk的重要程度）。表示协差阵的第i行i列元素，表示第k个特征向量的第i行。如果是利用相关阵进行主成分分析，则因子载荷量为：

其中表示相关阵i行i列，即1。

例如，对于人的身材的数据，通常包含身高、坐高、胸围、手臂长、肋围、腰围六个变量（指标）。通过主成分分析判断前三个主成分的累计贡献率大小可以达到要求，便选取这三个主成分。为了解释这三个主成分的意义，可以求各个主成分在原始的6个变量上的因子载荷量。如果主成分Y1对所有的变量的载荷都大致相等，就可以总结该主成分为身材大小成分；若第二主成分在胸围、肋围、腰围变量有相对较大或中等的正载荷量，在身高、坐高、手臂长有相对较大或中等的负载荷，就可以总结为胖瘦成分。所以，对提取的主成分的解释，除了借助因子载荷，也要根据实际问题进行总结，如果出现解释不了主成分，或者能够解释，但是主成分中的某些数据不符合实际情况，那也就意味着主成分分析的失败，此时只能换其它的方法了。

3.算法实现

兔兔这里以一种玉米的各个表现型数据为例，表现型分别为EH、SD、PH、FW、SFBN、IL、MEL、ELL、ELW、FLDW、ELA、TDW十二种，分别表示：穗位高（cm）、茎粗（cm）、株高（cm)、鲜质量（g）、雄花分枝数（个）、花序长度（cm）、雄穗长度（cm）、穗位叶长（cm）、穗位叶宽（cm）、功能叶干质量（g）、穗位叶面积（cm2）、总干质量（g）。

该数据集兔兔保存在corn.csv中，需要实际操作可以下载。我们这里就以相关阵的方法来进行主成分分析。

import numpy as np
import pandas as pd

df=pd.read_csv('corn.csv') #读取文件
df=pd.DataFrame(df)
R=df.corr() #样本相关阵
l,T=np.linalg.eig(R) #求特征值l与特征向量T，特征值默认从大到小排序
n,p=df.shape 
s=0
t=0
cr=[] #累计贡献率
for i in range(len(l)):
    t+=1
    contri=l[i]/np.sum(l) #第i主成分贡献率
    cr.append(contri)
    s+=contri #累计贡献率
    if s>=0.8: #累计贡献率达80%即停止
        break
pc=[] #主成分
for i in range(t):
    Y=np.dot(df,T[i].T) #计算第i主成分
    pc.append(Y)
factor_loading=[]
for i in range(t):
    a=[]
    for j in range(p):
        a.append(np.sqrt(l[i])*T[i][j]) #计算第i主成分在第j个变量上的载荷
    factor_loading.append(a)
factor_loading=np.array(factor_loading)
print('主成分个数：',t)
print('主成分：',np.mat(pc))
print('贡献率：',cr)
print('累计贡献率：',s)
print('因子载荷：',factor_loading)

通过上面的步骤，便可以求出各个数据。我们运行程序，发现只需要3个主成分便满足要求，并且累计贡献率接近88%，效果比较好。

通过最终的因子载荷，我们也能依据主成分在各个指标上的载荷大小来总结该主成分。不过我们发现，最终的主成分中含有负数值，但实际上，无论是叶片的面积，还是质量，或是总结成光合速率，都难以解释这些负数值，说明此时利用主成分分析是存在问题的，需要用其它的方法进行数据降维，如因子分析等。但如果不需要考虑实际意义，单纯就是进行数据降维，然后进一步做聚类、回归等，该方法也是可以使用的。

总结

主成分分析作为一种简便的数据降维方法，其应用还是十分广泛的，例如在基因组学、细胞学等生物学领域，或是其它学科，都可以找到主成分分析的实际应用。在使用主成分分析时需要了解其基本原理，掌握操作步骤，并能够对主成分进行合理解释，也需要了解该方法的缺陷。在实际的数据分析过程中，为了操作简便，也可采用SPSS等数据分析软件进行主成分分析，这样可以避免用代码实现该过程的繁琐。但如果需要真正理解主成分分析，深入理解数据分析、机器学习等知识，用代码实现这些公式的过程却是十分重要的。

详解主成分分析PCA

参考技术A

主成分分析（ Principal components analysis），简称PCA，是最主要的数据降维方法之一。本文从PCA的思想开始，一步一步推导PCA。

对于 , 。我们希望 从 维降到 维，同时希望信息损失最少。比如，从 维降到 ：

我们既可以降维到第一主成分轴，也可以降维到第二主成分轴。那么如何找到这这些主成分轴并且选择最优成分轴呢？

直观上，第一主成分轴 优于 第二主成分轴，即具有最大可分性。
下面解决一些基本概念。

欲获得原始数据新的表示空间，最简单的方法是对原始数据进行线性变换（基变换）：

其中 是原始样本， 是基向量， 是新表达。

数学表达：

其中 是行向量，表示第 个基， 是一个列向量，表示第 个原始数据记录.
当 时即 基的维度 < 数据维度时，可达到降维的目的。即：

以直角坐标系下的点(3,2)为例，欲将点(3,2)变换为新基上的坐标，就是用(3,2)与第一个基做内积运算，作为第一个新的坐标分量，然后用(3,2)与第二个基做内积运算，作为第二个新坐标的分量。

可以稍微推广一下，如果我们有m个二维向量，只要将二维向量按列排成一个两行m列矩阵，然后用“基矩阵”乘以这个矩阵，就得到了所有这些向量在新基下的值。例如(1,1)，(2,2)，(3,3)，想变换到刚才那组基上，则可以这样表示：

回顾一下，我们的目的是希望在降维过程中损失最少，换言之，我们希望投影后的数据尽可能分散开。这种分散程度可以用方差来表达， 方差 越大，数据越分散。

随机变量 表达了 的取值与其数学期望之间的偏离程度。若 较小，意味着 的取值主要集中在期望 也就是 的附近，反之，若 较大，意味着 的取值比较分散。

为了避免过于抽象，我们以一个具体的例子展开。假设我们5个样本数据，分别是 ，将它们表示成矩阵形式：

为了后续处理方便，我们首先将每个字段内所有值都减去字段均值，其结果是将每个字段都变为均值为0.

我们看上面的数据，设第一个特征为 ，第二个特征为 , 此时某一个样本可以写作：
且特征 的均值为2, 特征 的均值为3，所以变换后：

协方差 （Covariance）在 概率论 和 统计学 中用于衡量两个变量的总体 误差 。

比如对于二维随机变量 ，特征 除了自身的数学期望和方差，还需要讨论 之间互相关系的数学特征。

当 时，变量 完全独立，这也是我们希望达到的优化目标。

方差 是协方差的一种特殊情况，即当两个变量是相同的情况:

对于 二维 随机变量 ,

对于 n维 随机变量 ,

可见，协方差矩阵是 行 列的对称矩阵，主对角线上是方差，而协对角线上是协方差。

依然我们以一个具体的例子展开，还是这5个样本数据， , ，将它们去中心化后表示成矩阵形式：

那如果有 个样本的话，

对 做一些变换，用 乘以 的转置，并乘上系数1/m：

这不正是协方差矩阵嘛！

现在我们可以说：

回顾一下：

设 的协方差矩阵为 ， 的协方差矩阵为 ，且 。

我们要找的 不是别的，而是能让原始协方差矩阵对角化的 。

现在所有焦点都聚焦在了 协方差矩阵对角化 问题上。

由上文知道，协方差矩阵 是一个是对称矩阵，在线性代数上，实对称矩阵有一系列非常好的性质：

1）实对称矩阵不同特征值对应的特征向量必然正交。

2）设特征向量 重数为 ，则必然存在 个线性无关的特征向量对应于 ，因此可以将这 个特征向量单位正交化。

由上面两条可知，一个 行 列的实对称矩阵一定可以找到 个单位正交特征向量，设这 个特征向量为 ，我们将其按列组成矩阵：

则对协方差矩阵 有如下结论：

其中 为对角矩阵，其对角元素为各特征向量对应的特征值（可能有重复）。

结合上面的公式：

其中， 为对角矩阵，我们可以得到：

是协方差矩阵 的特征向量单位化后按行排列出的矩阵，其中每一行都是 的一个特征向量。如果设 按照 中特征值的从大到小，将特征向量从上到下排列，则用 的前 行组成的矩阵乘以原始数据矩阵 ，就得到了我们需要的降维后的数据矩阵 。

总结一下PCA的算法步骤：

设有 条 维数据。

1）将原始数据按列组成 行 列矩阵X

2）将 的每一行（代表一个特征）进行零均值化，即减去这一行的均值

3）求出协方差矩阵

4）求出协方差矩阵 的特征值及对应的特征向量

5）将特征向量按对应特征值大小从上到下按行排列成矩阵，取前 行组成矩阵

6） 即为降维到 维后的数据

这里以上文提到的：
，将它们表示成矩阵形式：

我们用PCA方法将这组二维数据其降到一维。

为了后续处理方便，我们首先将每个特征内所有值都减去字段均值，其结果是将每个字段都变为均值为0.

因为这个矩阵的每行已经是零均值，这里我们直接求协方差矩阵：

对于矩阵 :

和 分别是特征值和特征向量，
，则：

为了使这个方程式有非零解，矩阵 的行列式必须是 0 ：

即：

则：

分解得：

找到2个特征值， , ,

when :

即：

则：

和 可以取任意值，我们取归一化的 和 ，即： ,
此时 和

when :

即：

则：

和 可以取任意值，我们取归一化的 和 ，即：
此时 和

所以：

可以验证协方差矩阵C的对角化：

最后我们用 的第一行乘以数据矩阵，就得到了降维后的表示：

降维投影结果如下图：