主成分分析(PCA)及其在MATLAB中的实现

主成分分析 Principal Component Analysis(PCA)

 1. 什么是主成分分析? What is PCA?

以下是Wikipedia英文Principal component analysis词条的叙述:

Principal component analysis (PCA) is a statistical procedure that uses an orthogonal transformation to convert a set of observations of possibly correlated variables into a set of values of linearly uncorrelated variables called principal components. The number of principal components is less than or equal to the number of original variables.

意为:

主成分分析(PCA)是一种统计过程,通过一种正交变换将一组可能相关联的参数的观测值转换成一组线性互不关联的参数值,转换后的参数即被称为主成分。主成分的个数小于等于原参数的数目。

This transformation is defined in such a way that the first principal component has the largest possible variance (that is, accounts for as much of the variability in the data as possible), and each succeeding component in turn has the highest variance possible under the constraint that it is orthogonal to (i.e., uncorrelated with) the preceding components.

意为:

还定义这种变换具有如下特点:第一个主成分具有可能的最大方差值(也就是,它尽可能地包含数据中的变化性),随后的主成分依次具有可能的最大方差,且满足与之前的成分正交这一限制。

The principal components are orthogonal because they are the eigenvectors of the covariance matrix, which is symmetric. PCA is sensitive to the relative scaling of the original variables.

意为:

主成分是正交的,因为它们是协方差矩阵的特征向量,而协方差矩阵是一个对称矩阵。主成分分析对原参数的相关度比较敏感。

2. 主成分分析能做什么?What can PCA do?

以下是Wikipedia中文 主成分分析 词条的叙述:

在多元统计分析中,主成分分析(英语:Principal components analysis,PCA)是一种分析、简化数据集的技术。主成分分析经常用于减少数据集的维数,同时保持数据集中的对方差贡献最大的特征。这是通过保留低阶主成分,忽略高阶主成分做到的。这样低阶成分往往能够保留住数据的最重要方面。但是,这也不是一定的,要视具体应用而定。由于主成分分析依赖所给数据,所以数据的准确性对分析结果影响很大。

关键词:简化数据,减少数据集的维数

PCA是最简单的以特征量分析多元统计分布的方法。通常情况下,这种运算可以被看作是揭露数据的内部结构,从而更好的解释数据的变量的方法。如果一个多元数据集能够在一个高维数据空间坐标系中被显现出来,那么PCA就能够提供一幅比较低维度的图像,这幅图像即为在讯息最多的点上原对象的一个‘投影’。这样就可以利用少量的主成分使得数据的维度降低了。

关键词:数据的内部结构,解释数据

3.怎么做主成分分析?How to conduct/perform PCA?

这里需要说明:在Matlab和Wikipedia的Principal component analysis词条中,都将原数据矩阵的行定义为不同的实验/观测,将列定义为不同的参数/变量,这与其他一些文献中的定义方法刚好相反,所以会产生误解,比如,Wikipedia中covariance matrix词条中将行定义为参数/变量,本文中以Matlab和Wikipedia的Principal component analysis词条中的定义为准。

在维基百科英文词条下有具体的推导过程,本文不作详述,只介绍具体的计算方法,但所用的符号和数据的设定与英文词条下一致。

Consider a data matrix, X, with column-wise zero empirical mean (the sample mean of each column has been shifted to zero), where each of the n rows represents a different repetition of the experiment, and each of the p columns gives a particular kind of datum (say, the results from a particular sensor).

首先,假设有一个矩阵X,列向经验均值为零(即每列的样本均值移动到零),其n行表示被重复的n次实验,共p列的每一列表示一种特别种类的数据(比如说一种传感器的输出)。

The full principal components decomposition of X can therefore be given as

T=XW

where W is a p-by-p matrix whose columns are the eigenvectors of XTX。

于是,X的完整的主成分分解可由 T=XW 得出,其中W是一个pxp的矩阵它的列是矩阵XTX的特征向量。注意,这里的XTX与其他情况下协方差矩阵算式中的XXT其实是同一个矩阵,只不过两种情况下X的行列定义刚好相反。

具体的,先定义源数据矩阵如下:

PCA-orignal-matrix-definition

注意,此处Xi为行向量。

此后,计算过程如下:

PCA-calculation-procedure

具体的:

  • 第一步,计算Rx,即XTX,或XT的协方差矩阵(无1/(n-1)的系数)。
  • 第二步,求解Rx的特征值和特征向量。
  • 第三步,将特征值按从大到小顺序排列,同时将特征矩阵中的各列按对应特征值的顺序重新排列得到矩阵W
  • 第四步,计算排序完成后的特征值的累积贡献率accumulative contribution rate(ACR),设定阈值,如0.9,得到满足阈值的主成分数目。
  • 第五部,以L个主成分代替了p个变量,实现降维,且计算出源数据转换后的主成分的值。

主成分分析在Matlab中的实现 Perform PCA in Matlab

1. PCA算法的检验

首先我们根据PCA的定义算法进行计算,所用的源数据13*4的double表,13个观测4个变量,该数据通过matlab的hald数据获得,为ingredients矩阵。

load hald;
ingredients =
     7    26     6    60
     1    29    15    52
    11    56     8    20
    11    31     8    47
     7    52     6    33
    11    55     9    22
     3    71    17     6
     1    31    22    44
     2    54    18    22
    21    47     4    26
     1    40    23    34
    11    66     9    12
    10    68     8    12

在开始PCA前,首先需要将此矩阵列中心归零化,得到X

X = ingredients - repmat(mean(ingredients),13,1);
X =
   -0.4615  -22.1538   -5.7692   30.0000
   -6.4615  -19.1538    3.2308   22.0000
    3.5385    7.8462   -3.7692  -10.0000
    3.5385  -17.1538   -3.7692   17.0000
   -0.4615    3.8462   -5.7692    3.0000
    3.5385    6.8462   -2.7692   -8.0000
   -4.4615   22.8462    5.2308  -24.0000
   -6.4615  -17.1538   10.2308   14.0000
   -5.4615    5.8462    6.2308   -8.0000
   13.5385   -1.1538   -7.7692   -4.0000
   -6.4615   -8.1538   11.2308    4.0000
    3.5385   17.8462   -2.7692  -18.0000
    2.5385   19.8462   -3.7692  -18.0000

第一步,计算Rx

Rx = X'* X;
Rx =
   1.0e+03 *
    0.4152    0.2511   -0.3726   -0.2900
    0.2511    2.9057   -0.1665   -3.0410
   -0.3726   -0.1665    0.4923    0.0380
   -0.2900   -3.0410    0.0380    3.3620

第二步,求解特征值和特征向量:

[V,Lambda] = eig(Rx);
V =
    0.5062    0.5673    0.6460   -0.0678
    0.4933   -0.5440    0.0200   -0.6785
    0.5156    0.4036   -0.7553    0.0290
    0.4844   -0.4684    0.1085    0.7309
Lambda =
   1.0e+03 *
    0.0028         0         0         0
         0    0.1489         0         0
         0         0    0.8100         0
         0         0         0    6.2136

第三步,按特征值由大到小调整特征向量的列的顺序:

Lambda(4) > Lambda(3) > Lambda(2) > Lambda(1)
Lam(1) = Lambda(4,4);
Lam(2) = Lambda(3,3);
Lam(3) = Lambda(2,2);
Lam(4) = Lambda(1,1);
Lam =
 1.0e+03 *
 6.2136 0.8100 0.1489 0.0028
W(:,1) = V(:,4);
W(:,2) = V(:,3);
W(:,3) = V(:,2);
W(:,4) = V(:,1);
W =
   -0.0678    0.6460    0.5673    0.5062
   -0.6785    0.0200   -0.5440    0.4933
    0.0290   -0.7553    0.4036    0.5156
    0.7309    0.1085   -0.4684    0.4844

第四步,设定阈值为0.9,根据ACR值确定主成分数目:

ACR = cumsum(Lam)./sum(Lam);
ACR =
    0.8660    0.9789    0.9996    1.0000

由ACR的值可知,当主成分数目L=2时,ACR值>0.9的阈值,所以主成分数目为2即可。

第五步,计算源数据在主成分上的值:

T = X*W(:,1:2);
T =
   36.8218    6.8709
   29.6073   -4.6109
  -12.9818    4.2049
   23.7147    6.6341
   -0.5532    4.4617
  -10.8125    3.6466
  -32.5882   -8.9798
   22.6064  -10.7259
   -9.2626   -8.9854
   -3.2840   14.1573
    9.2200  -12.3861
  -25.5849    2.7817
  -26.9032    2.9310

至此,PCA计算过程结束。

2. Matlab函数的实现

在Matlab中进行PCA可通过princomp这一函数实现。

语法如下:

[COEFF,SCORE] = princomp(X)
[COEFF,SCORE,latent] = princomp(X)
[COEFF,SCORE,latent,tsquare] = princomp(X)
[...] = princomp(X,'econ')

其中COEFF表示变换系数,即W矩阵,SCORE表示源数据在主成分上的数值,即T矩阵,latent表示X的协方差矩阵的特征值,与之前的PCA定义中稍有不同,即相差1/(n-1)的倍数。

此外,princomp函数会自动对源数据矩阵进行列中心归零化,所以直接执行命令即可:

[pc,score,latent,tsquare] = princomp(ingredients);
pc =
   -0.0678   -0.6460    0.5673    0.5062
   -0.6785   -0.0200   -0.5440    0.4933
    0.0290    0.7553    0.4036    0.5156
    0.7309   -0.1085   -0.4684    0.4844
score =
   36.8218   -6.8709   -4.5909    0.3967
   29.6073    4.6109   -2.2476   -0.3958
  -12.9818   -4.2049    0.9022   -1.1261
   23.7147   -6.6341    1.8547   -0.3786
   -0.5532   -4.4617   -6.0874    0.1424
  -10.8125   -3.6466    0.9130   -0.1350
  -32.5882    8.9798   -1.6063    0.0818
   22.6064   10.7259    3.2365    0.3243
   -9.2626    8.9854   -0.0169   -0.5437
   -3.2840  -14.1573    7.0465    0.3405
    9.2200   12.3861    3.4283    0.4352
  -25.5849   -2.7817   -0.3867    0.4468
  -26.9032   -2.9310   -2.4455    0.4116
latent =
  517.7969
   67.4964
   12.4054
    0.2372
cumsum(latent)./sum(latent)
ans =
      0.86597
      0.97886
       0.9996
            1

可见,该结果与之前的五步法结果一致。

发表评论

电子邮件地址不会被公开。 必填项已用*标注

此站点使用Akismet来减少垃圾评论。了解我们如何处理您的评论数据