PCA模型

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++c2pXpZp=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=(x11x12x1px21x22x2pxn1xn2xnp)\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~=(xijxjˉ)/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ˉ=1ni=1nzi=1ni=1nj=1pcjxij~=1nj=1pcji=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)=1ni=1n(zizˉ)2=1ni=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)=(z11z12z1pz21z22z2pzn1zn2znp)\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=(λ100λ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 试验中取值变化的大小。如果 λi0{\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=1ni=1naijsj=1n1i=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~n1  (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λm0{\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 (pm){p \ (p \leq m)} 个主成分

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

bj=λjk=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λkk=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);

本博客所有文章除特别声明外,均采用 CC BY-SA 4.0 协议 ,转载请注明出处!