PCA模型

hotarugali

发布于 2022-03-16 15:31:27

8110

发布于 2022-03-16 15:31:27

文章被收录于专栏：hotarugaliの技术分享

1. 简介

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

2. 原理

2.1 基础思想

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

\begin{array}{c} s = c_1x_1 + c_2x_2 + \cdots + c_px_p \end{array}

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

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

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

\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}

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

\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}

其中，对于每一个 i ，均有 {c_{i1}^2 + c_{i2}^2 + \cdots + c_{ip}^2 = 1} 。且 {(c_{11},c_{12},\cdots,c_{1p})} 使得 {Var(Z_1)} 的值达到最大；{(c_{21},c_{22},\cdots,c_{2p})} 不仅垂直于 {(c_{11},c_{12},\cdots,c_{1p})} ，而且使得 {Var(Z_2)} 的值达到最大；{(c_{31},c_{32},\cdots,c_{3p})} 不仅垂直于 {(c_{11},c_{12},\cdots,c_{1p})} 和{(c_{21},c_{22},\cdots,c_{2p})} ，而且使 {Var(Z_3)} 的值达到最大；以此类推得到全部 p 个主成分。

【注意事项】

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

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

\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}

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

\begin{array}{c} \tilde{x_{ij} } = (x_{ij}-\bar{x_{j} })/s_j \end{array}

其中，\bar{x_j} 、s_j 分别为 X 的第 j 列的均值和标准差。变换后的数据矩阵为

\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 确定主成分变量

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

\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}

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

\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}

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

\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 = (c_1,c_2,\cdots,c_p)^T} ，则

\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} 来说，如果在n 次试验下它的取值变化不大，即 {z} 的 {M^*} 较小，则这个新变量 {z} 可以去掉。反之，如果 {z} 的 {M^*} 较大，则说明该新变量 {z} 的作用比较明显。因此，我们希望所选择的 {c_i\ (i=1,2,\cdots,p)}

{c_i\ (i=1,2,\cdots,p)}

能使 {M^*} 达到最大，即使得新变量 z 有较大作用。

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

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

其中 {\tilde{x} = (\tilde{x_1},\tilde{x_2},\cdots,\tilde{x_p})} 。

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

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

则第 {k+1} 个主成分 {Z_{k+1} = \tilde{x}^Tw} 可由以下两个条件决定：{w^T\eta_i = 0,\ i = 1,2,\cdots,k,\ w^Tw = 1} ；

在以上条件下，使得 {M^*} 达到最大。

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

2.3 选择主成分变量

确定了 p 个主成分后，将标准化后的数据 {\tilde{x_1},\tilde{x_2},\cdots,\tilde{x_p} } 转换成主成分 {z_1,z_2,\cdots,z_p} ，令

\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 = (\eta_1,\eta_2,\cdots,\eta_p)_{p \times p} } ，Q 为标准正交阵，且 {Z = XQ} ，则

\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}

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

【筛选方法】

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

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

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

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

\begin{array}{c} \rho_i = \sum_{j=1}^{r}r^2(z_j,x_i) \end{array}

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

3. 步骤

假设进行主成分分析的指标变量有 m 个：{x_1,x_2,\cdots,x_m} ；共有 n 个评价对象，第 i 个评价对象的第 j 个指标的取值为 {a_{ij} } ，评价矩阵为

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

3.1 对原始数据进行标准化

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

\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}

其中

\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}

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

\begin{array}{c} \tilde{x_i}={ {x_i-\mu_i} \over s_i} \ \ (i = 1,2,\cdots,m) \end{array}

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

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

3.2 计算相关系数矩阵

R相关系数矩阵 {R=(r_{ij})_{m \times 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}

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

3.3 计算特征值和特征向量

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

\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}

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

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

其中，{\tilde{X} = (\tilde(x_1),\tilde{x_2},\cdots,\tilde{x_m})^T} ，{Y = (y_1,y_2,\cdots,y_m)^T} ，{\mu = (\mu_1,\mu_2,\cdots,\mu_m)} 。标准化后的评价矩阵 {\tilde{A} } 转化为新评价矩阵 {\tilde{B} } ：

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

3.4 选择

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

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

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

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

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

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

3.5 计算综合评价值

定义综合得分为

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

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

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

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

4. 代码实现

4.1 MatLab 代码实现如下

利用 pcacov 函数：

% 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. 数据标准化
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);

本文参与腾讯云自媒体同步曝光计划，分享自作者个人站点/博客。

原始发表：2021-01-16，如有侵权请联系 cloudcommunity@tencent.com 删除

数据分析