1. 简介

主成分分析是指将数据中相关性很高的属性/变量转化成彼此相互独立或不相关的新属性/变量，利用较少的新属性/变量（主成分）去解释原来数据中的大部分属性/变量的一种降维方法。

2. 原理

2.1 基础思想

以学校课程为例，用 x1,x2,⋯ ,xp{x_1,x_2,\cdots,x_p} 表示 p 门课程的成绩变量， c1,c2,⋯ ,cp{c_1,c_2,\cdots,c_p} 表示各门课程的权重，则成绩的加权和为

s=c1x1+c2x2+⋯+cpxp\begin{array}{c} s = c_1x_1 + c_2x_2 + \cdots + c_px_p \end{array}

现在希望能够选择合适的权重来更好的区分学生的成绩 s1,s2,⋯ ,sn{s_1,s_2,\cdots,s_n}（ nn 为学生人数且 n>p{n \gt p}）。一般来说，如果这些值很分散，则表明区分的好。故需要寻找这样的加权，使得 s1,s2,⋯ ,sn{s_1,s_2,\cdots,s_n} 能够尽可能的分散。

而数据的分散程度可以利用方差的大小来体现，设 X1,X2,⋯ ,Xp{X_1,X_2,\cdots,X_p} 表示以 x1,x2,⋯ ,xp{x_1,x_2,\cdots,x_p} 为样本观测值的随机变量，规定（否则权值可以无穷大而失去意义） c12+c22+⋯+cp2=1{c_1^2 + c_2^2 + \cdots + c_p^2 = 1} 。在此约束下，若能找到 c1,c2,⋯ ,cp{c_1,c_2,\cdots,c_p} 使得 Var(c1X1+c2X2+⋯+cpXp){Var(c_1X_1 + c_2X_2 + \cdots + c_pX_p)} 的值达到最大，则表明此时学生的成绩区分达到最好。

可见，以上问题本质为一个求最优解的优化问题：

c12+c22+⋯+cp2=1max⁡ Var(c1X1+c2X2+⋯+cpXp)\begin{array}{c} c_1^2 + c_2^2 + \cdots + c_p^2 = 1 \\ \max \ Var(c_1X_1 + c_2X_2 + \cdots + c_pX_p) \end{array}

求出的最优解 c1,c2,⋯ ,cp{c_1,c_2,\cdots,c_p} 是一个单位向量，代表一个主成分方向，对应的主成分为 Z=c1X1+c2X2+⋯+cpXp{Z = c_1X_1 + c_2X_2 + \cdots + c_pX_p} 。
一个主成分不足以代表原来的 pp 个变量，故需要寻找第二、三、··· 个主成分。且后面的主成分不应该再包含前面已经求出来的主成分的信息，统计上理解就是主成分间两两之间的协方差为零，几何上理解就是主成分间两两正交。
设 Zi{Z_i} 表示第 ii 个主成分，i=1,2,⋯ ,p{i = 1,2,\cdots,p}，则

Z1=c11X1+c12X2+⋯+c1pXpZ2=c21X1+c22X2+⋯+c2pXp⋮Zp=cp1X1+cp2X2+⋯+cppXp\begin{array}{c} Z_1 = c_{11}X_1 + c_{12}X_2 + \cdots + c_{1p}X_p \\ Z_2 = c_{21}X_1 + c_{22}X_2 + \cdots + c_{2p}X_p \\ \vdots \\ Z_p = c_{p1}X_1 + c_{p2}X_2 + \cdots + c_{pp}X_p \end{array}

其中，对于每一个 ii，均有 ci12+ci22+⋯+cip2=1{c_{i1}^2 + c_{i2}^2 + \cdots + c_{ip}^2 = 1}。且 (c11,c12,⋯ ,c1p){(c_{11},c_{12},\cdots,c_{1p})} 使得 Var(Z1){Var(Z_1)} 的值达到最大；(c21,c22,⋯ ,c2p){(c_{21},c_{22},\cdots,c_{2p})} 不仅垂直于 (c11,c12,⋯ ,c1p){(c_{11},c_{12},\cdots,c_{1p})}，而且使得 Var(Z2){Var(Z_2)} 的值达到最大；(c31,c32,⋯ ,c3p){(c_{31},c_{32},\cdots,c_{3p})} 不仅垂直于 (c11,c12,⋯ ,c1p){(c_{11},c_{12},\cdots,c_{1p})}和(c21,c22,⋯ ,c2p){(c_{21},c_{22},\cdots,c_{2p})}，而且使 Var(Z3){Var(Z_3)} 的值达到最大；以此类推得到全部 pp 个主成分。

【注意事项】

  • 主成分分析的结果受量纲的影响，故实际应用中先将各变量的数据标准化，然后使用协方差矩阵或相关系数矩阵进行分析。
  • 使用相关系数矩阵求主成分时，Kaiser 主张将特征值小于 1 的主成分予以放弃（这也是 SPSS 软件的默认值）
  • 在实际研究中，由于主成分的目的是为了降维，减少变量个数，故一般选取少量的主成分（不超过 5 或 6 个），只要这些主成分的累计贡献率（定义见下文）能够达到 70%~80% 就行了。

设 pp 个变量 X1,X2,⋯ ,Xp{X_1,X_2,\cdots,X_p} 在第 ii 次试验中的取值为 xi1,xi2,⋯ ,xip(i=1,2,⋯ ,n){x_{i1},x_{i2},\cdots,x_{ip} (i = 1,2,\cdots,n)}，则原数据写成矩阵形式为

X=(x11x12⋯x1px21x22⋯x2p⋮⋮⋱⋮xn1xn2⋯xnp)\begin{array}{c} X = \left( \begin{matrix} x_{11} & x_{12} & \cdots & x_{1p} \\ x_{21} & x_{22} & \cdots & x_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ x_{n1} & x_{n2} & \cdots & x_{np} \end{matrix} \right) \end{array}

为了保证主成分分析不受量纲的影响，将 XX 进行标准化变换

xij~=(xij−xjˉ)/sj\begin{array}{c} \tilde{x_{ij} } = (x_{ij}-\bar{x_{j} })/s_j \end{array}

其中，xjˉ\bar{x_j}、sjs_j 分别为 XX 的第 jj 列的均值和标准差。
变换后的数据矩阵为

X~=(x11~x12~⋯x1p~x21~x22~⋯x2p~⋮⋮⋱⋮xn1~xn2~⋯xnp~)\begin{array}{c} \tilde{X} = \left( \begin{matrix} \tilde{x_{11} } & \tilde{x_{12} } & \cdots & \tilde{x_{1p} } \\ \tilde{x_{21} } & \tilde{x_{22} } & \cdots & \tilde{x_{2p} } \\ \vdots & \vdots & \ddots & \vdots \\ \tilde{x_{n1} } & \tilde{x_{n2} } & \cdots & \tilde{x_{np} } \end{matrix} \right) \end{array}

2.2 确定主成分变量

对于自变量的任意一个线性组合 ：

z=c1x1~+c2x2~+⋯+cpxp~,∑j=1pcj2=1\begin{array}{c} z = c_1\tilde{x_1} + c_2\tilde{x_2} + \cdots + c_p\tilde{x_p}, \sum_{j=1}^{p}c_j^2 = 1 \end{array}

其在第 ii 次试验中的取值为

zi=c1xi1~+c2xi2~+⋯+cpxip~ (i=1,2,⋯ ,n)\begin{array}{c} z_i = c_1\tilde{x_{i1} } + c_2\tilde{x_{i2} } + \cdots + c_p\tilde{x_{ip} }\ (i = 1,2,\cdots,n) \end{array}

由于 X~\tilde{X} 已经标准化了，因此

zˉ=1n∑i=1nzi=1n∑i=1n∑j=1pcjxij~=1n∑j=1pcj∑i=1nxij~=0\begin{array}{c} \bar{z}={1 \over n}\sum_{i=1}^{n}z_i = {1 \over n}\sum_{i=1}^{n}\sum_{j=1}^{p}c_j\tilde{x_{ij} } = {1 \over n}\sum_{j=1}^{p}c_j\sum_{i=1}^{n}\tilde{x_{ij} } = 0 \end{array}

记 w=(c1,c2,⋯ ,cp)T{w = (c_1,c_2,\cdots,c_p)^T}，则

M∗=Var(z)=1n∑i=1n(zi−zˉ)2=1n∑i=1nzi2=1n(X~w)T(X~w)\begin{array}{c} M^* = Var(z) = {1 \over n}\sum_{i=1}^{n}(z_i - \bar{z})^2 = {1 \over n}\sum_{i=1}^{n}z_i^2 = {1 \over n}(\tilde{X}w)^T(\tilde{X}w) \end{array}

对于新变量 z{z} 来说，如果在 nn 次试验下它的取值变化不大，即 z{z} 的 M∗{M^*} 较小，则这个新变量 z{z} 可以去掉。反之，如果 z{z} 的 M∗{M^*} 较大，则说明该新变量 z{z} 的作用比较明显。因此，我们希望所选择的 ci (i=1,2,⋯ ,p){c_i\ (i=1,2,\cdots,p)} 能使 M∗{M^*} 达到最大，即使得新变量 zz 有较大作用。

如果 XTX{X^TX} 的特征值 λ1≥λ2≥⋯≥λp{\lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_p}，它们对应的标准化正交特征向量为 η1,η2,⋯ ,ηp{\eta_1,\eta_2,\cdots,\eta_p}，则 M∗=(Xw)T(Xw)/n{M^* = (Xw)^T(Xw)/n} 的最大值在 w=η1{w = \eta_1} 时达到，且最大值为 λ1/n{\lambda_1 / n} 。此时新变量 z{z} 为

z=x~Tη1\begin{array}{c} z = \tilde{x}^T\eta_1 \end{array}

其中 x~=(x1~,x2~,⋯ ,xp~){\tilde{x} = (\tilde{x_1},\tilde{x_2},\cdots,\tilde{x_p})}。

常记 Z1=x~Tη1{Z_1 = \tilde{x}^T\eta_1} ，并称 Z1{Z_1} 为变量的第一主成分。
一般的，如果已经确定了 kk 个主成分：

Zi=x~Tηi  (i=1,2,⋯ ,k)\begin{array}{c} Z_i = \tilde{x}^T\eta_i\ \ (i = 1,2,\cdots,k) \end{array}

则第 k+1{k+1} 个主成分 Zk+1=x~Tw{Z_{k+1} = \tilde{x}^Tw} 可由以下两个条件决定：

  • wTηi=0, i=1,2,⋯ ,k, wTw=1{w^T\eta_i = 0,\ i = 1,2,\cdots,k,\ w^Tw = 1}；
  • 在以上条件下，使得 M∗{M^*} 达到最大。

由二次型的条件极值可知，第 k+1{k+1} 个主成分就是 Zk+1=x~Tηk+1  (i=1,2,⋯ ,p){Z_{k+1} = \tilde{x}^T\eta_{k+1}\ \ (i = 1,2,\cdots,p)} 。如此总共可以取到 pp 个主成分 Zi=x~Tηi (i=1,2,⋯ ,p){Z_i = \tilde{x}^T\eta_i \ (i = 1,2,\cdots,p)} 。

2.3 选择主成分变量

确定了 pp 个主成分后，将标准化后的数据 x1~,x2~,⋯ ,xp~{\tilde{x_1},\tilde{x_2},\cdots,\tilde{x_p} } 转换成主成分 z1,z2,⋯ ,zp{z_1,z_2,\cdots,z_p}，令

Z=X~(η1,η2,⋯ ,ηp)=(z11z12⋯z1pz21z22⋯z2p⋮⋮⋱⋮zn1zn2⋯znp)\begin{array}{c} Z = \tilde{X}(\eta_1,\eta_2,\cdots,\eta_p) = \left( \begin{matrix} z_{11} & z_{12} & \cdots & z_{1p} \\ z_{21} & z_{22} & \cdots & z_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ z_{n1} & z_{n2} & \cdots & z_{np} \end{matrix} \right) \end{array}

记 Q=(η1,η2,⋯ ,ηp)p×p{Q = (\eta_1,\eta_2,\cdots,\eta_p)_{p \times p} }，QQ 为标准正交阵，且 Z=XQ{Z = XQ}，则

ZTZ=QTXTXQ=QT(XTX)Q=(λ10⋱0λp)\begin{array}{c} Z^TZ = Q^TX^TXQ = Q^T(X^TX)Q = \left( \begin{matrix} \lambda_1 & & 0 \\ & \ddots & \\ 0 & & \lambda_p \end{matrix} \right) \end{array}

由上式可知，XTX{X^TX} 的特征值 λi{\lambda_i} 度量了第 ii 个主成分 zi{z_i} 在第 nn 试验中取值变化的大小。如果 λi≈0{\lambda_i \approx 0}，则说明该主成分 zi{z_i} 在 nn 次实验中的取值变化很小，可以选择剔除该主成分。

【筛选方法】

  • 首先将 XTX{X^TX} 的特征值按由大到小的次序排序。
  • 然后剔除这些特征值 λr+1,λr+2,⋯ ,λp{\lambda_{r+1},\lambda_{r+2},\cdots,\lambda_{p} }对应的主成分，这些特征值满足

∑i=r+1pλi<15%∑i=1pλi\begin{array}{c} \sum_{i=r+1}^{p}\lambda_i \lt 15\%\sum_{i=1}^{p}\lambda_i \end{array}

即筛选留下的主成分对应的特征值之和所占比重（定义为累计贡献率）要超过 85%{85\%}。

  • 单纯考虑累计贡献率有时是不够的，还需考虑选择的主成分对原始变量的贡献值，定义贡献值为相关系数的平方和。如果选择的主成分为 z1,z2,⋯ ,zp{z_1,z_2,\cdots,z_p} ，则它们对原始变量 xi~{\tilde{x_i} } 的贡献值为

ρi=∑j=1rr2(zj,xi)\begin{array}{c} \rho_i = \sum_{j=1}^{r}r^2(z_j,x_i) \end{array}

这里 r(zj,xi){r(z_j,x_i)} 表示 zj{z_j} 与 xi{x_i} 的相关系数。

3. 步骤

假设进行主成分分析的指标变量有 mm 个：x1,x2,⋯ ,xm{x_1,x_2,\cdots,x_m}；共有 nn 个评价对象，第 ii 个评价对象的第 jj 个指标的取值为 aij{a_{ij} }，评价矩阵为

A=(aij)n×m\begin{array}{c} A = (a_{ij})_{n \times m} \end{array}

3.1 对原始数据进行标准化

将各指标值转换成标准指标 aij~{\tilde{a_{ij} }} ：

aij~=aij−μjsj (i=1,2,⋯ ,n; j=1,2,⋯ ,m)\begin{array}{c} \tilde{a_{ij} } = { {a_{ij} - \mu_j} \over s_j}\ (i=1,2,\cdots,n;\ j=1,2,\cdots,m) \end{array}

其中

μj=1n∑i=1naijsj=1n−1∑i=1n(aij−μj)2j=1,2,⋯ ,m\begin{array}{c} \mu_j = {1 \over n}\sum_{i=1}^na_{ij} \\ s_j = {1 \over {n-1} }\sum_{i=1}^n(a_{ij}-\mu_j)^2 \\ j = 1,2,\cdots,m \end{array}

即 μj{\mu_j}，sj{s_j} 为第 jj 个指标的标准均值和样本标准差。对应地，称

xi~=xi−μisi  (i=1,2,⋯ ,m)\begin{array}{c} \tilde{x_i}={ {x_i-\mu_i} \over s_i} \ \ (i = 1,2,\cdots,m) \end{array}

为标准化指标变量。标准化后的评价矩阵为

A~=(aij~)n×m\begin{array}{c} \tilde{A} = (\tilde{a_{ij} })_{n \times m} \end{array}

3.2 计算相关系数矩阵 RR

相关系数矩阵 R=(rij)m×m{R=(r_{ij})_{m \times m} }，其中

rij=∑k=1naki~⋅akj~n−1  (i,j=1,2,⋯ ,m)\begin{array}{c} r_{ij} = { {\sum_{k=1}^n\tilde{a_{ki} } \cdot \tilde{a_{kj} } } \over {n-1} }\ \ (i,j=1,2,\cdots,m) \end{array}

易知式中 rii=1{r_{ii} = 1}，rij=rji{r_{ij} = r_{ji} }，rij{r_{ij} } 是第 ii 个指标与第 jj 个指标的相关系数。

3.3 计算特征值和特征向量

计算相关系数矩阵 RR 的特征值 λ1≥λ2≥⋯≥λm≥0{\lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_m \geq 0}，及对应的特征向量 μ1,μ2,⋯ ,μm{\mu_1,\mu_2,\cdots,\mu_m}，其中 μj=(μ1j,μ2j,⋯ ,μmj)T{\mu_j = (\mu_{1j},\mu_{2j},\cdots,\mu_{mj})^T}，由特征向量组成 mm 个新的指标变量

y1=μ11x1~+μ21x2~+⋯+μm1xm~y2=μ12x1~+μ22x2~+⋯+μm2xm~⋮ym=μ1mx1~+μ2mx2~+⋯+μmmxm~\begin{array}{c} y_1 = \mu_{11}\tilde{x_1} + \mu_{21}\tilde{x_2} + \cdots + \mu_{m1}\tilde{x_m} \\ y_2 = \mu_{12}\tilde{x_1} + \mu_{22}\tilde{x_2} + \cdots + \mu_{m2}\tilde{x_m} \\ \vdots \\ y_m = \mu_{1m}\tilde{x_1} + \mu_{2m}\tilde{x_2} + \cdots + \mu_{mm}\tilde{x_m} \end{array}

其中 yi{y_i} 表示第 ii 主成分。即

Y=X~Tμ\begin{array}{c} Y = \tilde{X}^T\mu \end{array}

其中，X~=((~x1),x2~,⋯ ,xm~)T{\tilde{X} = (\tilde(x_1),\tilde{x_2},\cdots,\tilde{x_m})^T}，Y=(y1,y2,⋯ ,ym)T{Y = (y_1,y_2,\cdots,y_m)^T}，μ=(μ1,μ2,⋯ ,μm){\mu = (\mu_1,\mu_2,\cdots,\mu_m)} 。
标准化后的评价矩阵 A~{\tilde{A} } 转化为新评价矩阵 B~{\tilde{B} }：

B~=A~μ\begin{array}{c} \tilde{B} = \tilde{A}\mu \end{array}

3.4 选择 p (p≤m){p \ (p \leq m)} 个主成分

计算特征值 λj (j=1,2,⋯ ,m){\lambda_j \ (j = 1,2,\cdots,m)} 的信息贡献率和累计贡献率。称

bj=λj∑k=1mλk  (j=1,2,⋯ ,m)\begin{array}{c} b_j = {\lambda_j \over {\sum_{k=1}^{m}\lambda_k} } \ \ (j = 1,2,\cdots,m) \end{array}

为主成分 yj{y_j} 的信息贡献率；称

αp=∑k=1pλk∑k=1mλk\begin{array}{c} \alpha_p = { {\sum_{k=1}^{p}\lambda_k} \over {\sum_{k=1}^m\lambda_k} } \end{array}

为主成分 y1,y2,⋯ ,yp{y_1,y_2,\cdots,y_p} 的累计贡献率，当 αp{\alpha_p} 接近1（αp=0.85/0.90/0.95{\alpha_p = 0.85 / 0.90 / 0.95}）时，选择前 pp 个指标变量 y1,y2,⋯ ,yp{y_1,y_2,\cdots,y_p} 作为 pp 个主成分，代替原来的 mm 个指标变量。选择该 pp 个主成分后，新评价矩阵 C~\tilde{C} 为：

C~=B~(n×p)\begin{array}{c} \tilde{C} = \tilde{B}_{(n \times p)} \end{array}

3.5 计算综合评价值

定义综合得分为

z=∑j=1pbjyj\begin{array}{c} z = \sum_{j=1}^pb_jy_j \end{array}

其中 bj{b_j} 为第 jj 个主成分的信息贡献；则 nn 个评价对象的综合评价值向量为

Z=C~b\begin{array}{c} Z = \tilde{C}b \end{array}

其中，b=(b1,b2,⋯ ,bp)T{b = (b_1,b_2,\cdots,b_p)^T} 。
然后便可以依据该评价值向量对这 nn 个对象进行综合评价。

4. 代码实现

4.1 MatLab 代码实现如下

  • 利用 pcacov 函数：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
% 1. 数据标准化
data=zscore(data);
% 2. 计算相关系数矩阵
R=corrcoef(data);
% 3. PCA 分析
% x：特征向量矩阵；y：特征值向量；z：主成分贡献率向量（总和为 100 ）
[x,y,z]=pcacov(R);
% 4. 选择 5 个主成分
p = 5;
% 5. 计算新评价矩阵
C = data*x;
C = C(:,1:p);
% 6. 计算综合评价值
Z = C*z(1:p)/100;
  • 利用 pca 函数：
1
2
3
4
5
6
7
8
9
10
11
% 1. 数据标准化
data=zscore(data);
% 2. PCA 分析
% x：特征向量矩阵；C：新评价矩阵；y：特征值向量
[x,C,y]=pca(data);
% 3. 选择 5 个主成分
p = 5;
% 4. 计算主成分贡献率（总和为 1 ）
z = y/sum(y);
% 5. 计算综合评价值
Z = C*z(1:p);