主成分分析与因子分析-原理及python实现

下面两种算法一般都需标准化消除量纲影响

主成分分析(PCA)

目的

数据降维，将n维数据降为n’维数据。原数据 X : n × m , s a m p l e   p o i n t : ( x 1 , . . . , x n ) T , b a s e : { w 1 , . . . , w n } X:n\times m,sample\,point:(x_1,...,x_n)^T,base:\lbrace w_1,...,w_n\rbrace X:n×m,samplepoint:(x1​,...,xn​)T,base:{w1​,...,wn​}

转换到n‘维空间中，$x^{(i)}\to z^{(i)}=(z_i^{(i)},...,z_{n^{\prime }}^{(i)}{)}^T$对应n’维空间标准正交基: w ′ = { w ^ 1 ′ , . . . , w ^ n ′ ′ } = W : n × n ′ w'=\lbrace \hat w'_1,...,\hat w'_{n'}\rbrace=W:n\times n' w′={w^1′​,...,w^n′′​}=W:n×n′

有:$W^T{x}^{(i)}=z^{(i)}$,n’维空间中数据还原到n维空间中有$W{z}^{(i)}={\bar{x}}^{(i)}$

caution:$W^{T}W=E,W{W}^T\ne E,W^{T}X=Z$

优化函数

使得降维后的数据与原数据距离和最小(尽可能维持原来位置，保持尽可能多的数据)
$$\begin{gathered} \sum _{i=1}^m||{\bar{x}}^{(i)}-x^{(i)}|{|}_2^2 \\ =\sum _{i=1}^m||W{z}^{(i)}-x^{(i)}|{|}_2^2 \\ =\sum _{i=1}^m(W{z}^{(i)}{)}^T(W{z}^{(i)})-2\sum _{i=1}^m(W{z}^{(i)}{)}^T{x}^{(i)}+\sum _{i=1}^m{x}^{(i)T}{x}^{(i)} \\ =-\sum _{i=1}^m{z}^{(i)T}{z}^{(i)}+\sum _{i=1}^{,}{x}^{(i)T}{x}^{(i)}=-tr(X^{T}W{W}^{T}X+X^{T}X) \\ min\sum _{i=1}^m||{\bar{x}}^{(i)}-x^{(i)}|{|}_2^2\iff argmaxtr(X^{T}W{W}^{T}X),s.t.W^{T}W=E \end{gathered}$$

用Lagrange函数求最优解:

因为:$tr(X^{T}W{W}^{T}X)=tr(W^{T}X{X}^{T}W)$

$J=tr(W^{T}X{X}^{T}W-\Lambda W^{T}W)=tr(W^T(X{X}^T-\Lambda I)W)$

$\frac{\partial J}{\partial W}=2(X{X}^T-\Lambda I)W=0\Rightarrow X{X}^{T}W=\Lambda IW={\Lambda }^{\prime }W$ 微分具体过程见Appendix

${\Lambda }^{\prime }$是对角矩阵，对角线值为$X{X}^T$的特征值

${\sum }_{i=1}^m||{\bar{x}}^{(i)}-x^{(i)}|{|}_2^2=-tr(\Lambda W^{T}W)-tr(X^{T}X)$

从而要使得上式最小则取$X{X}^T$最大的n’个特征值作为${\Lambda }_{ii}$，对应n’特征向量$({\hat{x}}_1,...{\hat{x}}_{n^{\prime }})$即为W

计算方式

• 将X每个样本关于每个特征维度中心化得到$\hat{X}$，从而$\hat{X}{\hat{X}}^T=cov(X,X)$

• 计算协方差矩阵特征值选择最大的n’个，对应特征向量组成W

• X转化 ：$Z=W^{T}X$

Appendix

• 对矩阵的迹求导计算方法可以先拆开成求和形式再求

在上述优化过程中求$W^T(X{X}^T-\Lambda I)W=W^{T}AW=P$

$tr(P)={\sum }_{p=1}^{n^{\prime }}{\sum }_{j=1}^n{\sum }_{i=1}^n{w}_{jp}{a}_{ji}{w}_{ip}$

$\frac{\partial tr(P)}{\partial w_{jp}}={\sum }_i{a}_{ji}{w}_{ip}=[AW{]}_{jp}$

$\frac{\partial tr(P)}{\partial w_{ip}}={\sum }_j{w}_{jp}{a}_{ji}=[A^{T}W{]}_{ip}$

$\frac{\partial tr(P)}{\partial W}=AW+A^{T}W=2AW$

因子分析(FA)

目的

找到影响现有几个特征的影响因子，方便进行聚类。

以计量经济学对中国家庭金融调查CHFS中金融知识相关问题(是否了解利率，是否正确计算利率，对投资风险是否了解，是否了解投资相关概念)，四个问题的回答构成样本$(x_1,x_2,x_3,x_4)$，其中12，34问题有极大相关性，则可以通过因子分析得到两个变量$(F_1,F_2)$分布表示受访者对利率和投资的认识程度。因子分析将$x_1,...,x_4$转换为$F_1,F_2$，并得到相关信息。

过程

sample point: $(x_1,...,x_n)$
$$\begin{gathered} \{\begin{array}{l}{x}_1=a_{11}{F}_1+...+a_{1n^{\prime }}{F}_{n^{\prime }}+{ϵ}_1\\ ...\\ x_n=a_{n1}{F}_1+...+a_{n{n}^{\prime }}{F}_{n^{\prime }}+{ϵ}_n\end{array} \\ x=AF+\mathrm{E} \end{gathered}$$

检验是否适合因子分析

• KMO检验

  $KM{O}_j=\frac{{\sum }_{i\ne j}{r}_{ij}^2}{{\sum }_{i\ne j}{r}_{ij}^2+{\sum }_{i\ne j}u}$

  $R=[r_{ij}]:correlationmatrix$

  $U=[u_{ij}]:partialcovariancematrix$

  $KMO=\frac{\sum {\sum }_{i\ne j}{r}_{ij}^2}{\sum {\sum }_{i\ne j}{r}_{ij}^2+\sum {\sum }_{i\ne j}u}$

  其中对应偏相关系数的计算参考

  u越大表示该变量在其它方向投影越分散，从而对应KMO值越小，一般总体KMO值大于0.7较好，计量中使用因子分析KMO较小的值为避免影响总体结果建议删去

• bartlett球形检验

  零假设为相关系数矩阵是单位矩阵，对应检验统计量(非常复杂)[参考](Bartlett’s test - Wikipedia)

因子的选择

因子分析中有很多确定因子变量的方法，如基于主成分模型的主成分分析和基于因子分析模型的主轴因子法、极大似然法、最小二乘法等，前者应用最为广泛

• 主成分分析中一般取$X{X}^T$前n’个最大的特征向量；在计量中，例如stata软件中默认取特征值大于1对应的因子
• 在中心化之后$X{X}^T$代表协方差,从而特征向量对应特征值即为方差贡献率，一般也取累计方差贡献率超过85%的前几个特征向量作为因子

因子旋转

• 特征向量表示原来n个变量在特征向量对应新因子中的构成，如果对于新因子F1，例子中是否了解利率，是否正确计算利率 两个变量占比最大，我们可以说F1可以对利率问题回答情况做出解释。
• 在上例中，如果F1,F2中四个问题占比没有明显差异，则为方便解释，使用旋转矩阵进行正交旋转，将F1,F2旋转为F1’,F2’使得对不同问题的解释性更突出。

因子旋转矩阵对应方差最大旋转矩阵，对应取值R为:
$$R_{Varmax}=argma{x}_R(\sum _{j=1}^{n^{\prime }}\sum _{i=1}^n(WR{)}_{ij}^4-\frac{\gamma }{n}\sum _{j=1}^{n^{\prime }}(\sum _{i=1}^n(WR{)}_{ij}^2{)}^2$$
即使得任何一个因子只在少数变量上有高载荷,而在其它变量上的载荷几乎为0。

因子得分

在类似主成分分析的因子选取后，得到新的n’个因子$z^{(i)}=W^T{x}^{(i)}$

通过n’个新因子的方差贡献率$\lambda$得到综合得分，即:

$scor{e}^i=\frac{{\sum }_{j=1}^{n^{\prime }}{\lambda }_j{z}_j}{{\sum }_{j=1}^{n^{\prime }}{\lambda }_j}$

如果进行因子旋转则按照旋转后对应的$\bar{W}=W\times R,R:n^{\prime }\times n^{\prime }$进行计算

程序

因子分析 stata程序

factor Q1_1 Q1_2 Q2_1 Q2_2 Q3_1 Q3_2, pcf
estat kmo
rotate
predict f1 f2 f3   #按照因子选择数
gen knowledge=(0.3083*f1+f2*0.2428+f3*0.1361)/(0.3083+0.2428+0.1361)#按照输出的因子载荷

主成分分析 python程序

参考

import numpy as np
from sklearn.decomposition import PCA
X = np.array([[-1, -1], [-2, -1], [-3, -2], [1, 1], [2, 1], [3, 2]])
pca = PCA(n_components=2)
newX = pca.fit_transform(X)     #=pca.fit(X) pca.transform(X)
invX = pca.inverse_transform(X)  #还原为Xbar
print(pca.explained_variance_ratio_)

因子分析 python程序

sklearn没有因子旋转相关函数，需要调用factor_analyzer这个包[具体API看文档](factor_analyzer package — factor_analyzer 0.3.1 documentation (factor-analyzer.readthedocs.io))

import factor_analyzer 
f=factor_analyzer.FactorAnalyzer
data=[[1,2,3,4],[23,33,44,55],[3,4,5,6],[12,44,56,77]]
f=f(rotation='varimax', n_factors=2, method='principal')
f.fit(data)

FactorAnalyzer(bounds=(0.005, 1), impute='median', is_corr_matrix=False,
               method='principal', n_factors=2, rotation='varimax',
               rotation_kwargs={}, use_smc=True)

 f.loadings_

array([[0.44773561, 0.89416599],
       [0.91066904, 0.41305279],
       [0.89686514, 0.44229149],
       [0.92472224, 0.38061092]])

f.get_communalities()#公因子方差:累计贡献率

array([0.99999999, 0.99993071, 0.99998885, 0.9999759 ])