期末速通主成分分析（PCA）

前言

       在数据分析的过程中，大部分情况我们会遇到很多变量，例如在分析英雄联盟的取胜之道以及各种因素对于游戏胜负的影响时，将会有golddiff（与对方整体经济差距）、champLevelDiff（与对方整体英雄等级差距）、wardsPlaced（队伍放置视野总数）等等诸多变量需要分析。在大数据时代，变量的个数更是只增不减。由于变量个数过多，所以在数据清洗后我们第一步需要做的就是降维，即用少数几个新的变量代替原有变量，合并重复信息，但不损失重要信息。多元统计分析中涉及两种经典降维方法，即主成分分析法以及因子分析。本文首先介绍主成分分析(PCA/Principal Component Analysis)。

一、主成分分析思想

（大白话 会的可以直接跳过）

       小学、初中、高中，各类周考、月考、期末考结束都会有一张令人深恶痛绝的成绩单，如果现在你是班主任，想看一下本班学生表现情况，你会怎么看捏？当然可以一门课一门课顺着看，但这必然是一种极其低效的办法，六门课你慢慢看，47门看到哪辈子去，通常正常情况下的想法大概是找一个或多个综合指标可以概括整体信息。
       而一个可以很大程度上将每个学生区分开的综合指标我们认为它是一个好的指标。回忆一下过去二十年被均分和总分支配的恐惧，这两种指标确实能为我们提供部分信息，但它仍具有局限性。设想一下，一个班里的同学总分差异并不大，但是有人语文、英语能考300，数学就考30；有的学生就反过来数学能考150，语文、英语两科一起就180，单看均分和总分来区分他们就很不合适了，所以此时就需要一些其他的指标将学生区分开。平均分本质上就是等权数计算，因此我们可以尝试通过给各科目赋予不同权数的方法，更好的区分学生。
       以上就是主成分分析直观上的理解。记语文成绩为$Y_1$、英语成绩为$Y_2$、数学成绩为$Y_3$，在所有可能的$Y_1$，$Y_2$，$Y_3$的线性组合中，找到一个或几个（通常小于3个），可以最大程度区分学生的线性组合或加权平均。

综上所述，主成分分析的思想可以简单概括为三点：
① PCA提供的是一个或几个综合指标
② 这些综合指标是由原来的变量通过线性组合或者加权平均构成的
③ PCA的目的是，最大程度地区分群体当中所有的个体

二、主成分分析（PCA）与线性判别分析（LDA）的区别

       LDA线性判别分析简单来说就是在寻找一个方向，可以使样本在这个方向的投影满足，不同种类的样本尽可能远，统一种类的样本尽可能靠近。如下图：

       从图中可以看出如果将样本点投影至y轴或x轴，两个类别都有大面积的重叠，而图中新找到的方向来看，两种颜色的数据分布较为集中，且类别之间的距离明显，显然更好的满足我们的要求。当然在实际案例分析中情况并不会这么简单，由于原始数据一般是超过二维的，所以投影后也大概率不是直线，而会是一个低维的超平面，具体情况具体分析。所以LDA的主要思想时最大化类间均值，最小化类内方差，即寻求最大化多个群体之间距离的线性组合。
       由第二部分的阐述可知PCA中的问题只有一个群体，我们的目标是找到使得这一个群体之间每个个体分的最开的变量组合。

三、总体主成分

理论推导

       由原始变量$Y_1,...,Y_p$到主成分$Z_1,...,Z_p$，需要满足什么条件?

① 线性组合 每一个$Z$是$Y$的线性组合或者说是加权平均
② 信息不重合 也就是$Z$之间是互不相关的，如果$Z_1$与$Z_2$之间相关性很高，就说明他们之间有信息重合的部分，冗余啊冗余，浪费空间要它干啥
③ 主成分按重要性排序 也就是每个主成分包含的信息含量大小，这里是用方差刻画的，方差越大包含的信息量越多，也就是越像原始数据，信息损失就越小。或者说方差越大，波动就越大，数据越分散，不确定性就越大，信息含量就越高。例如下面这种情况如果只能保留一个方向，必然选择$u_1$，很显然它和原始数据的分布最相似。

由上可知：（小写加粗是向量，矩阵乘法左行右列）

       $Z_1={\mathbf{a}}_1^{\prime }\mathbf{y}==a_{11}{Y}_1+a_{12}{Y}_2+...+a_{1p}{Y}_p$
       $Z_2={\mathbf{a}}_2^{\prime }\mathbf{y}==a_{21}{Y}_1+a_{22}{Y}_2+...+a_{2p}{Y}_p$
          $...$
       $Z_p={\mathbf{a}}_{\mathbf{p}}^{\prime }\mathbf{y}==a_{p1}{Y}_1+a_{p2}{Y}_2+...+a_{pp}{Y}_p$

       所以 $var(Z_j)=var({\mathbf{a}}_{\mathbf{j}}^{\prime }\mathbf{y})={\mathbf{a}}_{\mathbf{j}}^{\prime }\sum {\mathbf{a}}_{\mathbf{j}}$,   $cov(Z_j,Z_k)={\mathbf{a}}_{\mathbf{j}}^{\prime }\sum {\mathbf{a}}_{\mathbf{k}},k,j=1,...,p$.

       由于 主成分$Z_1,...,Z_p$需要按重要性排序，即按方差贡献度排序，目的是最大化方差，但是不想依赖于向量${\mathbf{a}}_{\mathbf{j}}$的长度使之最大化，所以增加了一个约束为 ${\mathbf{a}}_{\mathbf{j}}^{\prime }{\mathbf{a}}_{\mathbf{j}}=1$，即让$\frac{{\mathbf{a}}_{\mathbf{j}}^{\prime }\sum {\mathbf{a}}_{\mathbf{j}}}{{\mathbf{a}}_{\mathbf{j}}^{\prime }{\mathbf{a}}_{\mathbf{j}}}$最大化，${\mathbf{a}}_{\mathbf{j}}$在分母了，所以${\mathbf{a}}_{\mathbf{j}}\ne 0$.

（欸别害怕，只是符号比较多而已，本质思想其实非常简单，为了严谨才写那么多公式的，了解原理不记这些条条框框完全ok，耐心看完哈，向量$\mathbf{a}$只是个系数向量而已，就看做一个需要转置的常数好了）

       严谨的说就是：
第一主成分$Z_1={\mathbf{a}}_1^{\prime }\mathbf{y}$，在满足${\mathbf{a}}_1^{\prime }{\mathbf{a}}_1=1$情况下，最大化方差${\mathbf{a}}_1^{\prime }\sum {\mathbf{a}}_1$
第二主成分$Z_2={\mathbf{a}}_2^{\prime }\mathbf{y}$，在满足${\mathbf{a}}_2^{\prime }{\mathbf{a}}_2=1$，且${\mathbf{a}}_2^{\prime }\sum {\mathbf{a}}_1=0$情况下，最大化方差${\mathbf{a}}_2^{\prime }\sum {\mathbf{a}}_2$ （还要满足主成分之间不相关！！）
          $...$

第$j$主成分$Z_j={\mathbf{a}}_{\mathbf{j}}^{\prime }\mathbf{y}$，在满足${\mathbf{a}}_{\mathbf{j}}^{\prime }{\mathbf{a}}_{\mathbf{j}}=1$，且${\mathbf{a}}_{\mathbf{j}}^{\prime }\sum {\mathbf{a}}_{\mathbf{k}}=0,k=1,...,j-1$情况下，最大化方差${\mathbf{a}}_{\mathbf{j}}^{\prime }\sum {\mathbf{a}}_{\mathbf{j}}$

       当到第$p$个主成分时， ${\mathbf{a}}_{\mathbf{p}}^{\prime }\sum {\mathbf{a}}_{\mathbf{p}}$也就是最小的方差.

       所以所以，找主成分$Z_1,...,Z_p$，也就是找系数向量${\mathbf{a}}_1,...,{\mathbf{a}}_{\mathbf{p}}$.

       直接借助线性代数里的结论：对于任意向量$a$和$p\times p$对称矩阵$\sum$，记$({\lambda }_j,e_j),j=1,...,p$为$\sum$的特征对即特征值、特征向量，其中${\lambda }_1\ge ...\ge {\lambda }_p\ge 0$.

极值	极值点
${\max }_{x\ne 0}\frac{{\mathbf{a}}^{\prime }\sum \mathbf{a}}{{\mathbf{a}}^{\prime }\mathbf{a}}={\lambda }_1$	$\mathbf{a}=e_1$
${\max }_{x\ne 0,\mathbf{a}\perp e_1,...,e_{j-1}}\frac{{\mathbf{a}}^{\prime }\sum \mathbf{a}}{{\mathbf{a}}^{\prime }\mathbf{a}}={\lambda }_j$ 正交推导看下面	$\mathbf{a}=e_j$
${\min }_{x\ne 0}\frac{{\mathbf{a}}^{\prime }\sum \mathbf{a}}{{\mathbf{a}}^{\prime }\mathbf{a}}={\lambda }_p$	$\mathbf{a}=e_p$

（PCA本质还是坐标轴的旋转）
特征向量与特征值-3Blue1Brown

结论

记$({\lambda }_1,e_i),...,({\lambda }_p,e_p)$为协方差矩阵$\sum$的特征值-特征向量，${\lambda }_1\ge ...\ge {\lambda }_p\ge 0$，并且特征向量$e_1,....e_p$是标准化特征向量，则变量$Y_1,...,Y_p$的第$j$个主成分由下式给出：

                     $Z_j={\mathbf{e}}_{\mathbf{j}}^{\prime }\mathbf{y}={\mathbf{e}}_{\mathbf{j}1}{Y}_1+{\mathbf{e}}_{\mathbf{j}2}{Y}_2+...+{\mathbf{e}}_{\mathbf{j}\mathbf{p}}{Y}_p,j=1,...,p$

有$var(Z_j)={\mathbf{e}}_{\mathbf{j}}^{\prime }\sum {\mathbf{e}}_{\mathbf{j}}={\lambda }_j$       $cov(Z_j,Z_k)={\mathbf{e}}_{\mathbf{j}}^{\prime }\sum {\mathbf{e}}_{\mathbf{k}}=0,k,j=1,...,p$
并且${\sum }_{j=1}^{p}var(Z_j)={\sum }_{j=1}^{p}var(Y_j)$，即我们得到的主成分的所有方差的总的贡献度可以完全复原原始方差的维度.

${\sum }_{j=1}^{p}var(Z_j)={\sum }_{j=1}^{p}var(Y_j)$也就是可以还原原始变量的信息含量。原始变量Y的方差之和就是协方差矩阵的主对角线元素之和，也叫协方差矩阵$\sum$的迹。而特征值之和等于迹，于是上式成立。

系数向量$\mathbf{a}$为特征向量，方差为特征值.

第$k$个主成分贡献的方差，占总体方差的比例可表示如下:$$\frac{{\lambda }_k}{{\sum }_{j=1}^p{\lambda }_j},k=1,...,p$$
后续可以根据累计贡献率的值来选定主成分的个数.

标准化数据

       假如原始变量$Y_1,...,Y_p$的度量单位差别很大，所有导致每一个变量自身的波动性差异也很大，再通过主成分分析去寻找方差贡献率最大的方向的话，可能就会出现方差大的变量并不是它自身性质造成的，而是由于度量单位较小，数值变化较大所导致，此时由$\sum$生成的主成分会由方差大的做主导，显然不太何理。所以在主成分分析之前对原始变量分别做标准化、中心化，去除变量本身量纲的影响，并把坐标轴原点移动到原始变量中心的位置。
       需要注意的是，标准化后总体的协方差矩阵就是原总体$Y$的相关系数矩阵$\rho$，所以将上面结论中的对$\sum$求特征值、特征向量转化为对相关系数矩阵求特征值、特征向量即可。另外注意此时特征值之和即相关系数矩阵主对角线之和为$p$.
       如果在实际分析案例中变量的数值差异较大，推荐先对数据进行标准化之后获得一个较为均衡的数据集后用相关系数矩阵主成分分析。

四、样本主成分

       实际问题中的总体协方差矩阵$\sum$或相关系数矩阵$\rho$往往是未知的，需要通过样本数据估计。计算样本的协方差矩阵$S$和样本的相关系数矩阵$R$，再将总体主成分中的结论推广下来就好了。（方差是特征值，系数向量是特征向量）

主成分的命名：没有很明显的特征就不要随意命名，也不要随意把几个主成分合并在一起起一个名字。

五、主成分个数的选择

主成分分析的目的是在做降维，但是从$p$到$p$必然是没有达到降维的目的，所以为了有效归纳数据，达到降维的目的，怎么取、取几个个主成分比较合适，可以通过下面的方法：

① 百分比截点法 例如，累计方差贡献率第一次达到80%时截断
② 平均截点法 即保留特征值大于平均特征值所对应的主成分（如果是用相关系数矩阵做的话，这个平均值就是1，因为特征值之和等于矩阵的迹，即相关系数矩阵主对角线之和）
③ 碎石图 找拐点

六、R实例（结果解读）

       对52名学生的数学($x_1$)、物理($x_2$)、化学($x_3$)、语文($x_4$)、历史($x_5$)、英语($x_6$)，6科成绩进行主成分分析.

data = read.csv('score.csv', header = T);
R = round(cor(data),3); R
PCA_Cor = princomp(data, cor = T);PCA_Cor  # 基于相关系数矩阵
summary(PCA_Cor, loadings = T)
PCA_Cor$scores  # 主成分得分
screeplot(PCA_Cor, type = "line")  # 碎石图
load=loadings(PCA_Cor) #提取主成分载荷矩阵
plot(load[,1:2],xlim=c(-0.6,0.6),ylim=c(-0.6,0.6)) #作散点图
rnames=c("数学", "物理", "化学", "语文", "历史", "英语")
text(load[,1],load[,2],labels = rnames, adj=c(0.5, -0.5)) 
abline(h=0);  abline(v=0) #划分象限

由程序运行结果可知(保留三位小数)：
① 主成分的标准差（相关系数矩阵的6个特征值开方）分别为$\sqrt[2]{{\lambda }_1}=1.926$，$\sqrt[2]{{\lambda }_2}=1.124$，$\sqrt[2]{{\lambda }_3}=0.664$    … …
② 第一主成分的方差贡献率（信息量贡献率）为 $\frac{1.9261112^2}{6}=0.6183174$，第二主成分的方差贡献率（信息量贡献率）为 $\frac{1.1236019^2}{6}=0.2104135$   … …
③ 相关系数矩阵的特征向量分别为$\widehat{e_1}=(0.412,0.381,0.332,-0.461,-0.421,-0.430{)}^T$， $\widehat{e_2}=$   … …
④ 更具百分比截点法筛选了前两个主成分，或者绘制碎石图：

第一和第二主成分的方差较大，三之后都很小而且差不多大，所以只选取了第一、第二主成分，其他忽略，所以第一主成分和第二主成分分别为：
    $z_1^{\ast }=0.412x_1^{\ast }+0.381x_2^{\ast }+0.332x_3^{\ast }-0.461x_4^{\ast }-0.421x_5^{\ast }-0.430x_6^{\ast }$
    $z_2^{\ast }=0.376x_1^{\ast }+0.357x_2^{\ast }+0.563x_3^{\ast }+0.279x_4^{\ast }+0.415x_5^{\ast }+0.407x_6^{\ast }$

       结合案例分析，第一主成分前三个系数（数理化）为正，后三个系数（语数英）为负，且系数的绝对值都是0.3、0.4左右，所以第一主成分刻画的是理科和文科课程的成绩差异，即偏科的情况，第一主成分数值越大说明理科越好，反之则文科越好。
       第二主成分的系数都是正的，绝对值相差也不太多，所以可以近似的看作一个均衡的表现（类似于平均分），所以第二主成分数值越大说明整体成绩越好，反之则越不好。
       那每个学生的成绩在这两个指标上的表现如何，或者叫主成分得分能告诉我们的信息有：


       6号、7号学生的第一个指标也就是第一主成分得分很大，说明他们的理科成绩比文科成绩好很多，返回去看原始数据，6号的数学83、物理100、化学79、语文41、历史67、英语50，确实理科成绩更高一些；反之，49号学生的第一主成分得分很小，所以他的文科成绩要优于理科成绩，返回原始数据，49号的数学49、物理52、化学62、语文100、历史96、英语100，文科的分数是远远高于理科的。

       3号、5号同学的第二主成分得分较小，33号同学得分较大，通过对$z_2\ast$的分析可知，得分越高，整体成绩越好，返回原始数据，3号的数学67、物理63、化学49、语文65、历史67、英语57；33号的数学86、物理78、化学92、语文87、历史87、英语77. 33号同学的均衡表现优于3号.
       用主成分载荷矩阵的前两列数据绘制主成分载荷散点图如下：

可以更加直观的看出第一主成分确实具有明显的文理科差异特征。

注意

       在实际操作过程中会发现不是所有结果都这么完美，其实数据分析的第一步数据清洗是非常重要的，所以最费时间的还是能找出一组不错的数据。其次，princomp并不适用于那种样本量$n$小于变量个数$p$的情况，所以在吴喜之教授的《复杂数据统计方法——基于R的应用》中有更加贴近定义的code，感兴趣可以试一试，改起来不会很复杂。

data = read.csv('data/score.csv')

a = eigen(cor(scale(data[-1])))
round(a$values,3)
cca = (a$values)/sum(a$values);cca  # 贡献率
ca = cumsum(a$values /sum(a$values));ca  # 累计贡献率
plot(1:length(a$va),a$va,type = 'o', pch = 5,col=4,main = 'Scree Plot', xlab = 'Component Number', ylab='Eigen Value')
plot(1:length(a$va),ca,type = 'o', pch = 5,col=4,main = 'Cumulative', xlab = 'Component Number', ylab='Cumulative Contribution')
par(mfrow=c(1,1))
round(cca,3)
round(ca,3)

loadings=sweep(a$vectors,2,sqrt(a$values),'*')
b = loadings ;b
plot(b[,1:2],type='n',xlab = 'Component 1', ylab = 'Component 2')
text(b[,1],b[,2],names(data[-1]),cex=0.7)
abline(h=(max(b[,2])+min(b[,2]))/2,col=2)
abline(v=(max(b[,1])+min(b[,1]))/2,col=2)

write.table(round(b,3),file='载荷.csv', row.names=TRUE,col.names=TRUE);

# 主成分得分
sc = as.matrix(scale(data[,-1])%*%a$vectors[,1:4])  # 计算得分
plot(sc[,1],sc[,2],type='n',xlab = 'Component 1', ylab = 'Component 2')
text(sc[,1],sc[,2],data[,1],cex=1)
abline(v=(max(sc[,1])+min(sc[,1]))/2,col=2)
abline(h=(max(sc[,2])+min(sc[,2]))/2,col=2)

七、应用

       综上所述，主成分分析主要用于构建综合指标俩区分目标群体，例如构建顾客的各种消费行为的综合指标来对客户分级.
       高维的数据可视化其实也可以采用第一 、第二主成分散点图，通过二维的平面，直观表述数据的聚类信息等特征.
       通过主成分分析分析去处数据中的噪音，用重要的部分（主成分）重构原始数据，突出数据特征，所以可以用PCA做人脸识别.

（我滴天 终于结束了 公式好难敲 因子分析有缘再见 大概率无缘了属于是）