最近在学习STATA做分析的时候，发现这个软件很多功能很强大，但是背后的统计学知识要求也比较高，作为一边深入学习统计知识一遍用软件的小白，好多东西只是知其然不知其所以然，因此尝试自己把STATA的一些运算分解出来。因此这里记录一下学习内容。
在做STATA的主成分分析和因子分析的时候，觉得这两个东西很像，但是其中的原理好像也不太清楚，网上查了一些文章，花了不少时间才明白怎么做的，这里演示一下具体的做法。 对于这两个分析的说明，可以参考https://www.zhihu.com/question/23685740。https://wenku.baidu.com/view/33403b68ddccda38376baf6e.html
个人感觉主成分分析是因子分析的一个特例，或者说简单版本吧。先上代码：
df Out[32]: V2 V3 V4 V5 V6 V7 0 37066.0 26638.1 29218 11.2 7539 14210 1 52692.0 34634.4 30510 11.5 8394 14524 2 76909.0 46759.4 33261 12.4 9281 14608 3 91893.8 58478.1 35730 13.6 10077 15005 4 99595.3 67884.6 36454 14.0 10813 15733 5 113732.7 74462.6 38368 13.7 11356 16074 6 119048.0 78345.0 38046 12.5 11670 16100 7 126111.0 82067.0 40496 10.5 12393 16000 8 85673.7 89403.5 44452 10.0 13556 16300 df_uni = df.apply(lambda x: (x - x.mean()) / x.std()) df.corr() Out[34]: V2 V3 V4 V5 V6 V7 V2 1.000000 0.858753 0.735362 0.254789 0.783305 0.841657 V3 0.858753 1.000000 0.971990 -0.053640 0.988042 0.975661 V4 0.735362 0.971990 1.000000 -0.205558 0.992850 0.925418 V5 0.254789 -0.053640 -0.205558 1.000000 -0.183671 -0.009625 V6 0.783305 0.988042 0.992850 -0.183671 1.000000 0.953502 V7 0.841657 0.975661 0.925418 -0.009625 0.953502 1.000000 np.linalg.eig(df_uni.corr()) Out[35]: (array([ 4.62294858e+00, 1.15537221e+00, 1.64540586e-01, 5.39679774e-02, 2.44673331e-03, 7.23917617e-04]), array([[ 0.40431793, 0.33872688, 0.81928258, -0.14992877, -0.13851227, 0.09435506], [ 0.4645956 , 0.0020105 , -0.07957229, -0.0616411 , 0.43697537, -0.76358892], [ 0.45018537, -0.16110899, -0.34923608, -0.47114191, -0.65327173, -0.02593307], [-0.02844635, 0.91703669, -0.38331232, -0.09261985, 0.02411429, 0.04632759], [ 0.45863296, -0.12840259, -0.18572512, -0.18388754, 0.55460883, 0.63029293], [ 0.45481609, 0.04313279, -0.13804459, 0.84223154, -0.2343725 , 0.08911301]])) li =np.linalg.eig(df_uni.corr()) type(li) Out[37]: tuple li[0] Out[38]: array([ 4.62294858e+00, 1.15537221e+00, 1.64540586e-01, 5.39679774e-02, 2.44673331e-03, 7.23917617e-04]) li[1] Out[39]: array([[ 0.40431793, 0.33872688, 0.81928258, -0.14992877, -0.13851227, 0.09435506], [ 0.4645956 , 0.0020105 , -0.07957229, -0.0616411 , 0.43697537, -0.76358892], [ 0.45018537, -0.16110899, -0.34923608, -0.47114191, -0.65327173, -0.02593307], [-0.02844635, 0.91703669, -0.38331232, -0.09261985, 0.02411429, 0.04632759], [ 0.45863296, -0.12840259, -0.18572512, -0.18388754, 0.55460883, 0.63029293], [ 0.45481609, 0.04313279, -0.13804459, 0.84223154, -0.2343725 , 0.08911301]]) eig = li[0] vec = li[1] vec Out[42]: array([[ 0.40431793, 0.33872688, 0.81928258, -0.14992877, -0.13851227, 0.09435506], [ 0.4645956 , 0.0020105 , -0.07957229, -0.0616411 , 0.43697537, -0.76358892], [ 0.45018537, -0.16110899, -0.34923608, -0.47114191, -0.65327173, -0.02593307], [-0.02844635, 0.91703669, -0.38331232, -0.09261985, 0.02411429, 0.04632759], [ 0.45863296, -0.12840259, -0.18572512, -0.18388754, 0.55460883, 0.63029293], [ 0.45481609, 0.04313279, -0.13804459, 0.84223154, -0.2343725 , 0.08911301]]) sqr_eig = np.sqrt(eig) sqr_eig Out[47]: array([ 2.15010432, 1.07488242, 0.40563603, 0.23231009, 0.04946447, 0.02690572]) eig Out[48]: array([ 4.62294858e+00, 1.15537221e+00, 1.64540586e-01, 5.39679774e-02, 2.44673331e-03, 7.23917617e-04]) vec*sqr_eig Out[49]: array([[ 8.69325735e-01, 3.64091570e-01, 3.32330530e-01, -3.48299662e-02, -6.85143548e-03, 2.53869057e-03], [ 9.98929002e-01, 2.16105372e-03, -3.22773862e-02, -1.43198496e-02, 2.16147528e-02, -2.05449076e-02], [ 9.67945512e-01, -1.73173223e-01, -1.41662736e-01, -1.09451019e-01, -3.23137369e-02, -6.97747801e-04], [ -6.11626121e-02, 9.85706611e-01, -1.55485288e-01, -2.15165253e-02, 1.19280043e-03, 1.24647705e-03], [ 9.86108714e-01, -1.38017681e-01, -7.53367997e-02, -4.27189302e-02, 2.74334291e-02, 1.69584834e-02], [ 9.77902049e-01, 4.63626824e-02, -5.59958608e-02, 1.95658885e-01, -1.15931103e-02, 2.39764953e-03]]) 具体做法也很简单，对于给定的数据，先对其标准化处理，求其相关系数矩阵，之后求该矩阵的特征向量及特征值， 按照上面知乎答主的计算公式，计算出最终的结果。 该结果和在STATA中运行 factor V2-V7,pcf 出来的结果相同， 同时也是pca V2-V7 的结果，可以看出 pca 和 factor pcf的基本效果一样，需要原始数据的朋友可以在里面搜 STATA统计分析与行业应用案例详解。 注意最后一步不是矩阵和向量的乘法，而是矩阵的每一项单独乘以一个系数，具体参见链接中的原理公式。 从刚才的主成分分析中我们可以看到，前两个因子可以解释96%的方差，已经完全满足需要了，因此我们就保留这两个因子，之后进行varimax旋转，过程如下 f12 = [x[:2] for x in fs] f12 = np.array(f12) from numpy import eye, asarray, dot, sum, diag from numpy.linalg import svd def varimax(Phi, gamma = 1.0, q = 20, tol = 1e-6): p,k = Phi.shape R = eye(k) d=0 for i in xrange(q): d_old = d Lambda = dot(Phi, R) u,s,vh = svd(dot(Phi.T,asarray(Lambda)**3 - (gamma/p) * dot(Lambda, diag(diag(dot(Lambda.T,Lambda)))))) R = dot(u,vh) d = sum(s) if d/d_old < tol: break return dot(Phi, R) varimax(f12) array([[ 0.87227666, 0.35696404], [ 0.99891323, -0.00601514], [ 0.96649569, -0.18108995], [-0.05309267, 0.9861742 ], [ 0.98494602, -0.14608425], [ 0.97824877, 0.03835711]]) 最后得出的即时varimax之后的f1 f2值。 在stata里面运行rotate能得到相同结果。