初次发布于我的个人文档。（每次都是个人文档优先发布哦）
本文简要介绍一下主成分分析和因子分析这两个最常用的降维算法。
大题内容和之前发的主成分分析与因子分析没有区别，只是增加了R语言实现。
主成分分析PCA
如果你手上有一组数据，例如是大家的语文数学英语成绩。但是现在有一个问题，咱们的试卷出得有那么一点点不好，大家的成绩都集中在一起了，也就是试卷的区分度不大。现在，我们有没有办法补救呢？
注意：这只是一个例子而已，自然是不考虑我们进行各种变换之后的现实问题，例如这样搞成绩会不公平啊什么的。
总之，我们的核心问题是，有没有办法对现有数据进行变换，使得数据的每一个个体尽可能被分开。
这就是主成分分析的一个可以选择的切入点。
那我们要选择什么样的变换呢？以及有没有办法将一个群体之间的不同个体距离拉远。
以p维正态分布为例进行可行性探索
嗯，对我们先拿p维正态分布探索一下我们想法的可行性。
我们假设p维随机向量X服从协方差阵为$\Sigma$，均值向量为$\mu$的p维正态分布。
那么X的概率密度函数就是$P(x)=\frac{1}{(2\pi)^{\frac{p}{2}}|\Sigma|^{\frac{1}{2}}}e^{-\frac{1}{2}(x-\mu)^T\Sigma^{-1}(x-\mu)}$，我们来观察其概率密度等高线，显然，这里只有e的指数是变量，所以概率密度等高线满足：
$$-\frac{1}{2}(x-\mu)^T\Sigma^{-1}(x-\mu)=C_1$$也就是$(x-\mu)^T\Sigma^{-1}(x-\mu)=C_2$
我们可以对$\Sigma $进行谱分解。
谱分解说明：
根据线性代数的知识我们可以知道，任意实对称阵A可以正交相似对角化，即
$\forall 实对称矩阵A,\exist正交阵Q和对角阵\Lambda，使得A=Q\Lambda Q^T$
如果我们已知A的特征值$\lambda_1,\lambda_2,…,\lambda_n$和对应的特征向量$e_1,e_2,…,e_n$，则
Q=($e_1,e_2,…,e_n$)，$\Lambda =diag(\lambda_1,\lambda_2,…,\lambda_n)$
这意味着:
$$A=Q\Lambda Q^T$$$$=(e_1,e_2,...,e_n)diag(\lambda_1,\lambda_2,...,\lambda_n)(e_1,e_2,...,e_n)^T$$$$=(\lambda_1 e_1,\lambda_2 e_2,...,\lambda_n e_n)(e_1,e_2,...,e_n)^T$$$$=\lambda_1 e_1 e_1^T + \lambda_2 e_2 e_2^T +...+ \lambda_n e_n e_n^T$$$$=\sum_{i=1}^n\lambda_i e_i e_i^T$$这就是谱分解了。
我们假设$\Sigma$的特征值为$\lambda_1,\lambda_2,…,\lambda_p$，对应的特征向量为$e_1,e_2,…,e_p$，那么$\Sigma$就可以谱分解为
$\Sigma=\sum_{i=1}^p\lambda_i e_i e_i^T$，将他代入概率密度等高线方程就有，
$(x-\mu)^T{(\sum_{i=1}^p\lambda_i e_i e_i^T)}^{-1}(x-\mu)=C_2$
也就是，
$$\sum_{i=1}^p\frac{[e_i(x-\mu)]^T[e_i(x-\mu)]}{C_2}=1$$这是p维的类似椭圆的方程，当p=2时这就是椭圆方程。
这意味着，概率密度等高线是同样有着长轴和短轴。
因而，这意味着如果我们将原始变量X进行正交变换将坐标轴旋转到长轴上就可以达成我们的目标将一个群体之间的不同个体距离拉远。
演算
接下来，我们有了方向就可以进行推演了。
我们要将原始变量X进行正交变换得到新的一组变量，这从几何看就是进行坐标轴旋转。
总之，从代数角度看就是，设p维随机向量X=$(x_1,x_2,…,x_p)$的协方差阵为$\Sigma$。
那么我们就是要找一组新的变量Z=$(z_1,z_2,…,z_p)$使得（新变量被称为主成分）
$$z_1=a_{11}x_1+a_{12}x_2+...+a_{1p}x_p=a_1^TX$$$$z_2=a_{21}x_1+a_{22}x_2+...+a_{2p}x_p=a_2^TX$$…
$$z_p=a_{p1}x_1+a_{p2}x_2+...+a_{pp}x_p=a_p^TX$$而此时，$var(z_j)=a_j^T\Sigma a_j,cov(z_j,z_k)=a_j^T\Sigma a_k$
我们前面说了，我们希望将一个群体之间的不同个体距离拉远，也就是要最大化新变量Z的方差，与此同时我们自然希望各个新变量之间无关也就是：
最大化$var(z_j)$，希望$cov(z_j,z_k)=0$
对于$z_1$来说就是希望最大化$a_1^T\Sigma a_1$，但是显然我们可以通过无限扩大$a_1$的长度来实现最大化$z_1$的方差，这是我们不期望看到的。
所以我们再额外要求$z_1$的长度是1即$a_1^Ta_1=1$。
这样的话其实我们就是在最大化$a_1^T\Sigma a_1=\frac{a_1^T\Sigma a_1}{a_1^Ta_1}$（注意哦，现在分母为1所以除了等于没除）。
类似地，我们对$z_2$会要求$a_2^Ta_2=1$，并且$cov(z_2,z_1)=a_2^T\Sigma a_1=0$
最大化$a_2^T\Sigma a_2$。
以此类推，但是到最后一个变量$z_p$我们只能要求最小化$a_p^T\Sigma a_p$了，因为这个对应的是前面说的高维椭圆的短轴，是最小的。
那么，怎么进行最小化呢？
一般的教材这里就是上拉格朗日乘数法了，计算比较复杂我就不说了。给个结论吧。
$\forall a \in R^p,\Sigma \in M_p且\Sigma为对称矩阵。$
设$(\lambda_j,e_j)$为$\Sigma$的特征值、单位特征向量。
$$a\neq 0, a⊥e_1,e_2,...,e_{j-1},$$$$max \frac{a^T\Sigma a}{a^T a}=\lambda_j在a=e_j时取到最大。$$
利用上述结论就可以知道，X的第j个主成分$z_j=e_j^TX$
且$var(z_j)=e_j^T\Sigma e_j$
注意$(\lambda_j,e_j)$为$\Sigma$的特征值、单位特征向量。
所以$\Sigma e_j =\lambda_j e_j,e_j^Te_j=1$
因而$var(z_j)=e_j^T\Sigma e_j=e_j^T \lambda_j e_j=\lambda_j e^T_je_j=\lambda_j$。
所以，X的第j个主成分$z_j$是其协方差矩阵$\Sigma$的第j个单位特征向量乘以原始变量X，并且第j个主成分的方差就是$\Sigma$第j个特征值$\lambda_j$。
这就是主成分分析的结论。
更进一步的，如果我们对原始变量进行标准化然后再进行主成分分析，可以证明这相当于对原始变量的相关系数矩阵R进行对应的主成分分析。
降维
从上面的推导我们可以看到对p维向量进行主成分分析只能得到p个主成分，似乎不能降维啊。那么我们一般说的降维是怎么回事？
前面我们知道，第j个主成分的方差就是$\Sigma$第j个特征值$\lambda_j$。
如果我们将全部的主成分的方差求和，那就是对$\Sigma$全部的特征值的求和，也就是$\Sigma$的迹，也就是X各个分量的方差的和。
所以到现在为止我们还没有损失任何一点点方差。
如果降维的话就会损失方差了，这是因为所谓的降维就是将各个特征值从大到小排列，然后去掉比较小的特征值和对应的主成分。
这样的话就会损失方差了，也就损失了部分信息。这就是所谓的利用主成分分析进行降维。
R语言实现
R语言作为为统计而生的语言，实现主成分分析非常简单。
|
|
核心就是这里的prcomp函数。
当然，我们已经知道主成分分析的结果是样本协方差矩阵的特征向量特征值。
所以也可以用上面的公式手动求解。
|
|
因子分析
那么因子分析是什么？
还是看学生成绩数据吧，从学生的成绩上我们可以看到，优秀的学生似乎各科成绩都很好。也许你还会发现，各科成绩高度相关，这意味着他们可能由某一个潜在变量决定（智商）。
因子分析就是由原始数据寻找这样的潜变量。
由于潜变量的数量往往少于原始变量的数量，所以因子分析也是一种降维方法。
因子分析建立了一个因子模型，它认为原始变量Y是各个潜变量的线性组合，即
$$Y_i=l_{i1}F_1+l_{i2}F_2+...+l_{im}F_m+\epsilon_i$$其中，$F_j$是潜变量也叫公共因子，我们假设有m个，当然一般要求m不大于原始变量的个数p。
系数$l_ij$被称为因子载荷，$Y_i$则是原始变量而$\epsilon_i$是类似误差的特殊因子。
我们还对公共因子提了一些基础的要求，首先$F_i,F_j$不相关（正交），$F_i,\epsilon_j、\epsilon_i,\epsilon_j$不相关，来保持各个变量之间的独立性。
当然，因子模型用矩阵表示更简洁，就是
$Y=AF+\epsilon$
那前面说的那些要求用矩阵表示就是：
- $m \le p$
- cov(F,$\epsilon$)=0
- $D_F=var(F)=单位阵I_m$
- $D_\epsilon=var(\epsilon)=diag(\sigma_1^2,\sigma_2^2,…,\sigma_p^2) $
而因子模型最重要的是协方差阵的矩阵分解：
$var(X)=\Sigma=cov(AF+\epsilon,AF+\epsilon)=Acov(F,F)A^{-1}+cov(\epsilon,\epsilon)=AA^{-1}+D_\epsilon$
演算
那么如何求解因子模型中未知的A和$\epsilon$呢？
答案是利用协方差矩阵的矩阵分解：
$\Sigma=AA^{-1}+D_\epsilon$
而前面我们说过$\Sigma$可以分解为$Q\Lambda Q^T$，其中$Q=($e_1,e_2,…,e_n$)，$$\Lambda =diag(\lambda_1,\lambda_2,…,\lambda_n)$
我们再进行小小的变换，定义$\Lambda_2=(e_1\sqrt{\lambda_1},e_2\sqrt{\lambda_2},…,e_p\sqrt{\lambda_p})$，
则$\Sigma=\Lambda_2 \Lambda_2$。
其实就是在数值计算方法里说的，利用矩阵的谱分解对这个矩阵进行平方根分解。
当$\Sigma$的后p-m个特征值很小的时候，我们就可以忽略掉后面的项，用$\Lambda_2$的前m项估计A，从而$D_\epsilon=\Sigma-AA^{-1}$也就可以计算了。
这里需要注意一下，我们的因子分析也是一种降维的方法，你仔细看前面的因子模型，我们假设有m个潜变量，并且我们要求m小于原来的原始变量的数量n。
也就是说，这里的矩阵A只有m列而不是n列，因此其实这里也不完全是所谓的对$\Sigma$的平方根分解，而是只取了平方根分解的前m列。
你也可以理解为，我们暂时认为$D_\epsilon$是零，从而A就是$\Sigma$的平方根分解，或者说我们用$\Sigma$的平方根分解估计A。
而当我们估计出A后，又可以利用估计出的A迭代得到更好的对$D_{\epsilon}$的估计。只能说，这种玩法在数值计算方法里玩了很多次了。
主轴因子法
如果想进一步玩这种方法，那你会得到主轴因子法。
我们知道
$\Sigma=AA^{-1}+D_\epsilon$
前面说的主成分法实际上是先估计$D_{\epsilon}=0$，然后再估计A，再用估计的A得到对$D_{\epsilon}$的更好的估计。
而当你得到了对$D_{\epsilon}$更好的估计，你又可以倒回去估计A了，然后就循环了。你可以一直这样迭代，直到结果收敛。
在主轴因子法里，我们对$D_{\epsilon}$的初始估计未必是0，这也使得主轴因子法更泛用。
那如果我们对$D_{\epsilon}$进行了非0的估计，A就是$\Sigma-D_{\epsilon}$的平方根分解了。
不过，我们在讨论主轴因子法的时候往往先将原始变量进行标准化从而协方差矩阵$\Sigma=原始变量的相关系数矩阵R$
那么所谓的$\Sigma-D_{\epsilon}$就是约相关系数矩阵$R^*=R-D_\epsilon=AA^{-1}$，
至于这里$D_\epsilon$是怎么估计的，常见的有4种方法但是其实吧都有点无厘头的。
第一种方法就是取$D_{\epsilon}=0$，也就是前面的主成分法。
至于其他的方法他们都是从共同度的角度进行估计的。
所有你得先理解共同度是什么才行。
因子载荷矩阵的统计意义
共同度来自于对因子载荷矩阵的统计意义的研究。
观察因子模型，$Y=AF+\epsilon$
我们将它写开，观察A的每个元素的意义。
$$Y_i=l_{i1}F_1+l_{i2}F_2+...+l_{im}F_m+\epsilon_i$$计算看看$Y_i$的方差。
$D(Y_i)=D(l_{i1}F_1+l_{i2}F_2+…+l_{im}F_m+\epsilon_i)=l_{i1}^2D(F_1)+l_{i2}^2D(F_2)+…+l_{im}^2D(F_m)+D(\epsilon_i)$
前面我们在因子模型里假设，$D_F=var(F)=单位阵I_m$，从而$D(F_i)=1$，所以上式实际上就是
$D(Y_i)=l_{i1}^2+l_{i2}^2+…+l_{im}^2+\sigma_i^2$
为了方便呢，我们假设A的第i行元素的平方和，也就是$l_{i1}^2+l_{i2}^2+…+l_{im}^2$为共同度$h_i^2$
从而每一个原始变量$Y_i$的方差可以分解为第i共同度$h_i^2$加上第i特殊因子$\epsilon$的方差$\sigma_i^2$。
而如果竖着来看，我们计算A的第i列元素的平方和，那你会得到第i个因子对总的方差的贡献量，我们经常用它来衡量每个公共因子的重要性，不过这与本文无关。
本文想要的只是一个结果，$D(Y_i)=h_{i}^2+\sigma_i^2$
在上一部分，也就是主轴因子法，我们已经对原始变量进行了标准化，从而原始变量的方差都是1。
那么就有：$h_{i}^2+\sigma_i^2=1$
在主成分法里我们估计$\sigma_i^2=0$其实也就是估计$h_i^2=1$。
也就是说，我们给出$\sigma_i^2$的估计和给出共同度$h_i^2$的估计是等价的。
那在其他估计方法里，我们都是给出共同度的估计，其实也等价于同步给出了$\sigma_i^2$的估计。
例如，我们估计$h_{i}^2$为原始变量$X_i$与其他原始变量的复相关系数的平方，或者估计共同度为相关系数矩阵R的第i行的绝对值最大的非对角元的绝对值，或者估计共同度为$h_{i}^2=\frac{r_{ik}r_{il}}{r_{kl}}$,其中$r_{ik},r_{il}$为原始变量相关系数矩阵R的第i行上最大的两个非对角元。
但是你说要我给出这么估计的好处或者理由的话，其实也没有什么特别好的理由，而且这些估计方法，，，，
其实很有可能导致约相关系数矩阵有负的特征值，这样得出的结果很有可能是我们无法解释的。
所以我个人认为这些方法其实，没有那么的优秀。
因子分析的R语言实现
这个实现也非常简单，并且R语言内置的因子分析可以估计因子得分和进行因子旋转。
简单说一下这两个东西。这两个东西的理论非常复杂，结果也很难看，所以我也不打算细说。
前面我们的讨论只是得到了载荷矩阵A，但是有些时候我们还想得到因子模型中的F，这个F就被我们称为因子得分。
这玩意儿的算法一般可以使用Thomson回归法。最终的结果是$F=A^{T}R^{-1}Y$
而因子旋转，则是因为因子模型不唯一。
观察一下因子模型，$Y=AF+\epsilon$
显然我们可以对潜变量作正交变换T，使得$Y=AT^{-1}TF+\epsilon$
这时，我们可以得到另一组潜变量和载荷矩阵即$TF$和$AT^{-1}$。也就是说，因子分析的结果不唯一，他们直接可能差了一个正交变换。
而正交变换在几何上表示为旋转，也因此才有了因子旋转的名字。
所谓因子旋转就是将我们刚刚算出来的因子F，进行适当的正交变换，使得因子分析的结果解释性更好。
我们常用的因子旋转是所谓的最大方差旋转法，它在选择的时候保证能够最大化载荷矩阵中每一列载荷的平方的方差最大化。
当然，你确实也可以以其他目标进行旋转，但是对因子分析来说我们是要求最终分析的结果能够被解释的。也就是说，这个结果是要能够被人们理解的，而不是脱离现实的。
|
|