最近在学习STATA做分析的时候，发现这个软件很多功能很强大，但是背后的统计学知识要求也比较高，作为一边深入学习统计知识一遍用软件的小白，好多东西只是知其然不知其所以然，因此尝试自己把STATA的一些运算分解出来。因此这里记录一下学习内容。

在做STATA的主成分分析和因子分析的时候，觉得这两个东西很像，但是其中的原理好像也不太清楚，网上查了一些文章，花了不少时间才明白怎么做的，这里演示一下具体的做法。
对于这两个分析的说明，可以参考 [ https://www.zhihu.com/question/23685740

](https://www.zhihu.com/question/23685740) 。 [
https://wenku.baidu.com/view/33403b68ddccda38376baf6e.html

](https://wenku.baidu.com/view/33403b68ddccda38376baf6e.html)

个人感觉主成分分析是因子分析的一个特例，或者说简单版本吧。先上代码：

```code
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的基本效果一样，需要原始数据的朋友可以在CSDN里面搜 STATA统计分析与行业应用案例详解。

 

注意最后一步不是矩阵和向量的乘法，而是矩阵的每一项单独乘以一个系数，具体参见链接中的原理公式。

 

从刚才的主成分分析中我们可以看到，前两个因子可以解释96%的方差，已经完全满足需要了，因此我们就保留这两个因子，之后进行varimax旋转，过程如下

```code
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]])
[/code]

 

 

最后得出的即时varimax之后的f1 f2值。 在stata里面运行rotate能得到相同结果。

 

 

 

 

 

 

 

 

```