主成分与因子分析
人们往往希望能够找出少数具有代表性的变量来对复杂事物进行描述，这需要把反映该事物的很多变量或数据进行高度概括。主分成分析和因子分析便是如何利用复杂多样的数据来综合描述客观事物特征的分析方法和过程。
数据降维
基本问题
把反映一个事物特征的多个变量用较少且具有代表性的变量来描述，这个过程称之为数据降维。不同的变量往往是从不同的侧面或方面去描述事物特征的，这些不同方面称之为事物的维度。当反映事物方面太多的时候，过多的数据会对所描述对象造成混乱，很难的道正确结论。因此，应当把相关的维度进行总结概括，尽量降低数据为度，简要地对事物特征进行描述。
为了能够简要而不遗漏地反映事物特征，数据降维过程中应当解决如下几个基本问题
1.能否把数据的多个变量用较少的综合变量表示？
2.较少的综合变量包含有多少原来的信息？
3.能否利用找到的综合变量来对事物进行较为全面的分析？
解决好这些基本问题之后，就可以用简化的数据对事物进行描述或判定，从而得出统计分析的结论。
基本原理
二维数据，有两个变量，可由二维坐标轴上的横坐标和纵坐标来表示，因此每个观测值都有相应于这两个坐标轴的两个坐标值。在正态分布的假定下，这些数据在二维坐标轴上形成一个椭圆分布形状。
众所周知，椭圆有一个长轴一个短轴，且相互垂直。在短轴上，数据变化很少；长轴方向，数据变化的范围较大。长轴就是要找的主要综合变量。这样，由二维到一维的降维过程就完成了。
当坐标轴和椭圆的长短轴平行，那么代表长轴的变量就描述了数据的主要变化，而代表短轴的变量就描述了数据的次要变化。但是，坐标轴通常并不和椭圆的长短轴平行。因此，需要寻找椭圆的长短轴，即进行坐标平移或旋转变换，使得新变量和椭圆的长短轴平行。如果长轴变量代表了数据包含的大部分信息，就用变量在该轴上的变化代替原先的两个变量（舍去次要的另一个维度），降维就完成了。椭圆的长短轴相差得越大，降维效果就越好。
对于多维变量的降维，主要从高维椭球入手。首先把高维椭球的主轴找出来，再用代表大多数数据信息的最长的几个轴作为新变量。与二维椭圆分布形状类似，高维椭球的主轴也是互相垂直的，这些互相相交的新变量是原先变量的线性组合，可叫做主成分。
主成分选择的个数，有一定的选择标准，就是这些被选中主成分所代表的主轴的长度之和与主轴长度总和的比值，这个比值称之为阈值。根据相关文献，所选的主轴总长度之和占所有主轴长度之和约85%即可。但在实际应用过程中，要依据研究目的、研究对象和所手机的变量具体情况而定。
主轴越长，表示变量在该主轴方向上的变动程度越大，亦即方差越大。所以一般情况下，不去计算主轴的长度，而是计算其主轴方向的方差，根据所选取的主轴的方差之和与所有主轴方向上方差之和的比值，即方差贡献率的大小来判断应该取多少个主成分。
主成分分析
基本概念和原理
主成分分析（principal component）是数据降维的基本方法之一。主成分是由原始变量提取的综合变量，可以用如下式子来表示
\[\begin{align*}
& Y_1= \mu_{11}x_1 + \mu_{12}x_2 + \cdots + \mu_{1p}x_p \\
& Y_2= \mu_{21}x_1 + \mu_{22}x_2 + \cdots + \mu_{2p}x_p \\
& \vdots \\
& Y_p= \mu_{p1}x_1 + \mu_{p2}x_2 + \cdots + \mu_{pp}x_p
\end{align*}
\]
其中，\(Y\) 表示主成分，\(x\)为原始变量，\(\mu_{ij}\) 为系数，有约束条件：\(\mu_{k1}^2+\mu_{k2}^2+\cdots+\mu_{kp}^2=1\)，\(\mu_{ij}\) 可由原始数据协方差矩阵或相关系数矩阵确定。
在提取出来的各个主成分当中，\(Y_i,Y_j\)相互无关，且第一个主成分 \(Y_1\) 是 \(x_1,x_2,\cdots,x_p\) 的一切线性组合最大的；第二个主成分 \(Y_2\) 是 \(x_1,x_2,\cdots,x_p\) 的一切线性组合第二大的；以此类推，第 \(n\) 个主成分\(Y_n\) 是 \(x_1,x_2,\cdots,x_p\) 的一切线性组合第 \(n\) 大的。
由原始数据的协方差阵或相关系数矩阵，可计算出矩阵的特征值或特征根：
\[\lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_p
\]
其中，\(\lambda_1\) 对应 \(Y_1\) 的方差，\(\lambda_2\) 对应 \(Y_2\) 的方差，...，\(\lambda_p\) 对应 \(Y_p\) 的方差，因此有：
\[阈值 = \frac{被选择的主成分长度}{主轴成分总和}=\frac{选择的特征根的和}{特征根总和}=累积方差贡献率
\]
\(\lambda\) 对应的特征向量 \(\mu\) 就是主成分分析线性模型中对应的系数，如：\(\lambda_1\) 对应的特征向量为 \(\mu_{11},\mu_{12},\cdots,\mu_{1p}\) 为第1个主成分的线性组合系数，即 \(Y_1=\mu_{11}x_1 + \mu_{12}x_2 + \cdots + \mu_{1p}x_p\)。
这些系数称为主成分载荷(loading)，它表示主成分和相应原始变量的相关系数。相关系数绝对值越大，主成分对该变量的代表性也越大。根据上式计算出来的 \(Y\) 值称之为主成分得分。
在实际问题中，不同的变量往往有不同的量纲，为了不同量纲数据之间的可比性，保证所提取主成分与原始变量意义上的一致性，在进行主成分分析之前，可按照如下Z-Score公式将变量标准化或无量纲化：
\[x_i^* = \frac{x_i-E(x_i)}{\sqrt{Var(x_i)}}(i=1,2,\cdots,p)
\]
其中，\(E(x_i)\) 表示原始变量 \(x_i\) 的期望，\(Var(x_i)\) 表示 \(x_i\) 的方差。
基本步骤和过程
在主成分分析的过程中，通常要先把各变量进行无量纲化（即标准化）。把变量进行标准化后，可按照如下顺序进行主成分分析
1.选择协方差阵或相关阵计算特征根及对应的特征向量
2.计算方差贡献率，并根据方差贡献率阈值选取合适的主成分个数
3.根据主成分载荷大小对选取主成分进行命名
4.根据主成分载荷计算各个主成分得分
示例
import pandas as pd
import statsmodels.api as sm
import scipy.stats as stats
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
# 为评价地区综合发展水平，用人均GDP、人均可支配收入、人均消费指出等数据进行综合考察
live = pd.read_csv(\'./data/live.csv\', encoding=\'gb2312\')
print(live.head(5))
# District GDP Income Consumption Employment Education Health Life
# 0 北京 45444 17652.95 13244.20 0.3937 584.43 1295.76 76.10
# 1 山西 12495 8913.91 6342.63 0.2554 548.83 538.70 71.65
# 2 内蒙古 16331 9136.79 6928.60 0.2158 504.77 533.36 69.87
# 3 吉林 13348 8690.62 6794.71 0.1836 502.08 675.77 73.10
# 4 黑龙江 14434 8272.51 6178.01 0.2418 479.85 613.15 72.37
# 7个变量过于复杂，进行主成分降维
# # 方法一：matplotlib：使用SVD
# # matplotlib3中移除了PCA
# # from matplotlib.mlab import PCA as mlabPCA
# #
X = live.iloc[:, 1:8]
# live_pca1 = mlabPCA(X, standardize=True) # standardize表示是否将原始数据标准化V
# # mlab的PCA实例对象的fracs属性表示每个特征值占特征值总和的百分比，即每个主成分对应的方差贡献率
# # 为了更好地依据这些结果来决定应当取多少个主成分合适，将结果输出
# live_var = pd.DataFrame((live_pca1.s) / np.mean(live_pca1.s),
# columns=[\'Eigenvalue\'], index=list(range(1, 8)))
# s = 0
# p, c = [], []
# for i in range(0, len(live_pca1.fracs)):
# s += live_pca1.fracs[i]
# p.append(live_pca1.fracs[i])
# c.append(s)
# live_var[\'Proportion\'] = p
# live_var[\'Cumulative\'] = c
# print(live_var)
# # 相关系数矩阵计算的特征根Eigenvalue,对应方差贡献率Proportion,累计贡献率Cumulative
# # 分析过程的碎石图
# fig, (ax1, ax2) = plt.subplots(1, 2)
# fig.set_size_inches(6.8, 3)
# fig.subplots_adjust(wspace=0.4)
#
# ax1.plt(range(1, 8), live_var[\'Eigenvalue\'], \'o-\')
# ax1.set_title(\'Scree Plot\')
# ax1.set_xlabel(\'Principal Components\')
# ax1.set_ylabel(\'Eigenvalue\')
# ax1.grid()
#
# ax2.plt(range(1, 8), live_var[\'Proportion\'], \'o-\')
# ax2.plt(range(1, 8), live_var[\'Cumulative\'], \'bo-\')
# ax2.set_title(\'Variance Explained\')
# ax2.set_xlabel(\'Principal Components\')
# ax2.set_ylabel(\'Proportion\')
# ax2.grid()
# plt.show()
# # 展现出各特征根大小及其贡献
#
# # 提取主成分并计算主成分得分
# # 计算特征根对应的特征向量，将其作为原始变量线性组合的系数
# live_eigenvectors = pd.DataFrame(live_pca1.Wt,
# index=[\'Prin1\', \'Prin2\', \'Prin3\', \'Prin4\', \'Prin5\', \'Prin6\', \'Prin7\'],
# columns=(live.columns)[1:8]).T
# print(live_eigenvectors)
# # 根据各主成分对应的特征向量即系数，可以计算出各策划给你分的得分。
# # 可以根据主成分计算公式中的系数，即主成分载荷绝对值的大小来判定该主成分所主要代表原始变量的含义
# # 把原始变量标准化后，带入公式，计算出每个样本在对应主成分的得分
# print(live.iloc[20:21, :]) # 20号样本
# print(live_pca1.Y[20]) # 自动计算各样本在各个主成分上的得分
# print(live_pca1.project([...])) # 对新的原始数据计算其主成分得分
# 方法二：传统统计，与SAS、SPSS等一致
def PCA(x, components=None):
if components == None:
components = x.size / len(x) # 如果components参数未指定，则赋值为原始变量个数
average = np.mean(x, axis=0)
sigma = np.std(x, axis=0, ddof=1)
r, c = np.shape(x)
data_standardized = []
mu = np.tile(average, (r, 1))
data_standardized = (x - mu) / sigma # 使用Z-Score法对原始数据标准化
cov_matrix = np.cov(data_standardized.T) # 计算协方差矩阵
EigenValue, EigenVector = np.linalg.eig(cov_matrix) # 求解协方差矩阵的特征值和特征向量
index = np.argsort(-EigenValue) # 按照特征值大小降序排列
Score = []
Selected_Vector = EigenVector.T[index[:int(components)]]
Score = data_standardized * np.matrix(Selected_Vector.T)
return EigenValue[index], Selected_Vector, np.array(Score)
# 计算特征值贡献及累计贡献率
EigenValue, Vector, Score = PCA(np.asarray(X))
live_ev = pd.DataFrame((EigenValue), columns=[\'EigenValue\'], index=list(range(1, 8)))
prop = live_ev[\'EigenValue\'] / live_ev[\'EigenValue\'].sum()
s = 0
p, c = [], []
for i in range(1, len(prop) + 1):
s += prop[i]
p.append(prop[i])
c.append(s)
live_ev[\'Proportion\'] = p
live_ev[\'Cumulative\'] = c
print(live_ev)
# EigenValue Proportion Cumulative
# 1 4.725499 0.675071 0.675071
# 2 1.234341 0.176334 0.851406
# 3 0.448662 0.064095 0.915500
# 4 0.306114 0.043731 0.959231
# 5 0.213755 0.030536 0.989767
# 6 0.060574 0.008653 0.998421
# 7 0.011054 0.001579 1.000000
# 按照特征根累计贡献率85%的阈值，前两个特征根已经代表了原始数据的大部分信息，且前两个远远大于其余特征根贡献。故取2个主成分
# 特征值Eigenvalue,对应方差贡献率Proportion,累计贡献率Cumulative
# 特征值对应的特征向量
live_ev = pd.DataFrame(Vector,
index=[\'Prin1\', \'Prin2\', \'Prin3\', \'Prin4\', \'Prin5\', \'Prin6\', \'Prin7\'],
columns=(live.columns)[1:8]).T
print(live_ev)
# Prin1 Prin2 Prin3 ... Prin5 Prin6 Prin7
# GDP 0.441618 0.073883 0.083499 ... -0.325047 -0.796130 0.171583
# Income 0.447192 -0.029164 -0.193227 ... -0.358861 0.213764 -0.765499
# Consumption 0.435590 -0.016302 -0.394963 ... -0.276029 0.450039 0.611568
# Employment 0.122961 0.827743 0.098366 ... 0.231263 0.143879 -0.030374
# Education 0.365034 -0.397744 -0.255464 ... 0.718524 -0.101607 -0.056051
# Health 0.374018 0.307351 -0.001732 ... 0.345642 -0.062374 -0.010942
# Life 0.356365 -0.235799 0.851325 ... 0.011964 0.288153 0.079822
# 显示各主成分对应的特征向量的系数
# 各主成分得分
live_S = pd.DataFrame(live[\'District\'])
live_S[\'Prin1_Score\'] = Score[:, 0]
live_S[\'Prin2_Score\'] = Score[:, 1]
live_S[\'Score\'] = Score[:, 0] * 0.675071 + Score[:, 1] * 0.176334 # 计算公式
res = live_S.sort_values(by=\'Score\', ascending=False)
print(res)
# District Prin1_Score Prin2_Score Score
# 0 北京 5.697979 3.551375 4.472769
# 5 上海 6.049685 -1.500794 3.819326
# 7 浙江 3.817791 -1.245874 2.357590
# 6 江苏 1.154256 -1.300429 0.549895
# 8 福建 0.540295 0.011427 0.366753
# 9 山东 0.422840 -0.329803 0.227292
# 15 重庆 0.412800 -1.337303 0.042857
# 4 黑龙江 -0.571680 0.418921 -0.312054
# 3 吉林 -0.327915 -0.535612 -0.315813
# 1 山西 -0.581030 0.404957 -0.320829
# 19 陕西 -0.493203 0.063999 -0.321662
# 2 内蒙古 -0.667531 -0.004352 -0.451398
# 12 湖南 -0.463060 -0.841579 -0.460997
# 23 * -1.259577 2.027835 -0.492728
# 11 湖北 -0.879394 -0.332502 -0.652285
# 10 河南 -1.136978 0.384593 -0.699724
# 22 宁夏 -1.297820 0.535104 -0.781764
# 13 广西 -0.938053 -0.876802 -0.787862
# 20 甘肃 -1.503304 0.600191 -0.909003
# 16 四川 -1.273157 -0.600263 -0.965318
# 18 * -1.633711 0.663251 -0.985917
# 17 云南 -1.578860 0.373034 -1.000064
# 14 海南 -1.593914 -0.459780 -1.157080
# 21 青海 -1.896461 0.330405 -1.221984
# 绘图主成分载荷图
# 指定为黑体中文字体，防止中文乱码
plt.rcParams["font.sans-serif"] = ["Heiti TC"]
# 解决保存图像是负号\'-\'显示为方块的问题
plt.rcParams[\'axes.unicode_minus\'] = False
fig, ax = plt.subplots(1)
ax.plot(live_S[\'Prin1_Score\'], live_S[\'Prin2_Score\'], \'o\')
ax.set_xlabel(\'Component 1 (economy life)\')
ax.set_ylabel(\'Component 2 (employment)\')
ax.axvline(live_S[\'Prin1_Score\'].mean(), color=\'k\', ls=\'--\')
ax.axhline(live_S[\'Prin2_Score\'].mean(), color=\'k\', ls=\'--\')
dotxy = tuple(zip(live_S[\'Prin1_Score\'] - 0.2, live_S[\'Prin2_Score\'] + 0.15))
i = -1
for dot in dotxy:
i += 1
ax.annotate(live_S.iloc[i][\'District\'], xy=dot)
plt.show()
# 方法三：sklearn
from sklearn.decomposition import PCA as skPCA
from sklearn.preprocessing import scale
x = scale(X) # 数据标准化Z-Score转换函数
live_pca2 = skPCA(n_components=len(X.columns)).fit(x)
# 参数n_components可省略，此处保留全部主成分，拟合主成分模型时应采用标准化转换数据
# 特征值方差贡献率
res = live_pca2.explained_variance_ratio_
print(res) # [0.67507132 0.17633444 0.06409455 0.04373058 0.03053648 0.00865347 0.00157915]
# 各个主成分特征值对应的特征向量，即主成分载荷
res = live_pca2.components_
print(res)
# [[ 0.44161842 0.44719158 0.4355896 0.12296062 0.3650339 0.37401802
# 0.35636501]
# [ 0.07388346 -0.02916376 -0.01630157 0.82774347 -0.39774384 0.30735092
# -0.23579859]
# [-0.08349857 0.19322704 0.39496334 -0.09836587 0.25546362 0.00173219
# -0.8513253 ]
# [-0.15369971 -0.03697985 -0.03546836 -0.46361551 -0.33696562 0.80135152
# -0.05569072]
# [-0.32504735 -0.35886119 -0.27602861 0.23126301 0.71852416 0.34564199
# 0.01196397]
# [ 0.79612966 -0.21376358 -0.4500389 -0.14387884 0.10160667 0.06237445
# -0.28815288]
# [ 0.17158256 -0.76549902 0.61156775 -0.03037356 -0.05605106 -0.01094167
# 0.07982213]]
# 二维数组第一维度表示主成分，第而维度表示各变量在主成分上的载荷
# 各个样本主成分得分(第1和第2主成分)
res = live_pca2.transform(x)[:, :2]
print(res)
# [[ 5.82053058e+00 3.62775703e+00]
# [-5.93527114e-01 4.13666985e-01]
# [-6.81887725e-01 -4.44541795e-03]
# [-3.34967853e-01 -5.47131743e-01]
# [-5.83975125e-01 4.27930723e-01]
# [ 6.17980083e+00 -1.53307269e+00]
# [ 1.17908199e+00 -1.32839801e+00]
# [ 3.89990352e+00 -1.27267020e+00]
# [ 5.51915946e-01 1.16731070e-02]
# [ 4.31934797e-01 -3.36896398e-01]
# [-1.16143241e+00 3.92864715e-01]
# [-8.98307512e-01 -3.39653118e-01]
# [-4.73019121e-01 -8.59679414e-01]
# [-9.58228194e-01 -8.95659830e-01]
# [-1.62819568e+00 -4.69669275e-01]
# [ 4.21678287e-01 -1.36606580e+00]
# [-1.30053942e+00 -6.13173495e-01]
# [-1.61281735e+00 3.81056842e-01]
# [-1.66884829e+00 6.77516578e-01]
# [-5.03811182e-01 6.53758715e-02]
# [-1.53563722e+00 6.13099918e-01]
# [-1.93725027e+00 3.37511610e-01]
# [-1.32573357e+00 5.46612801e-01]
# [-1.28666790e+00 2.07144922e+00]]
# 结果与之前的有误差，原因是标准化时采用的标准差不同，skearn中的sacle采用分母为样本量n,
# 而通常使用满足无偏性的样本修正标准差，分母是样本量n-1.
# 对原始数据重新使用样本修正标准差进行标准化
x1 = (X - X.mean()) / X.std(ddof=1)
res = live_pca2.transform(x1)[:, :2]
print(res)
# [[ 5.69797937e+00 3.55137464e+00]
# [-5.81030407e-01 4.04957231e-01]
# [-6.67530587e-01 -4.35181972e-03]
# [-3.27915108e-01 -5.35611890e-01]
# [-5.71679535e-01 4.18920646e-01]
# [ 6.04968518e+00 -1.50079386e+00]
# [ 1.15425643e+00 -1.30042860e+00]
# [ 3.81779109e+00 -1.24587414e+00]
# [ 5.40295361e-01 1.14273298e-02]
# [ 4.22840415e-01 -3.29803048e-01]
# [-1.13697846e+00 3.84592953e-01]
# [-8.79393656e-01 -3.32501725e-01]
# [-4.63059708e-01 -8.41578873e-01]
# [-9.38052709e-01 -8.76801721e-01]
# [-1.59391403e+00 -4.59780393e-01]
# [ 4.12799855e-01 -1.33730330e+00]
# [-1.27315658e+00 -6.00263134e-01]
# [-1.57885950e+00 3.73033694e-01]
# [-1.63371071e+00 6.63251473e-01]
# [-4.93203443e-01 6.39993832e-02]
# [-1.50330439e+00 6.00191106e-01]
# [-1.89646148e+00 3.30405307e-01]
# [-1.29782027e+00 5.35103875e-01]
# [-1.25957713e+00 2.02783488e+00]]
# 方法四：mdp
import mdp
live_pca3 = mdp.pca(np.asarray(x1))
# mdp.pca可以直接得到各主成分得分，返回值是各主成分得分的一个数组
print(live_pca3)
# [[-5.69797937e+00 3.55137464e+00 -3.43926693e-01 7.28717798e-01
# 1.96285126e-01 -1.96635693e-01 6.91697964e-02]
# [ 5.81030407e-01 4.04957231e-01 -5.04362737e-01 -4.01002658e-01
# -4.76385921e-01 -2.03051981e-02 -1.28643669e-01]
# [ 6.67530587e-01 -4.35181972e-03 6.25705985e-02 -2.05617575e-02
# 9.56282728e-02 3.59643170e-01 2.88141106e-02]
# [ 3.27915108e-01 -5.35611890e-01 -7.44555275e-01 8.68521699e-01
# -1.73664358e-01 5.86590509e-02 1.47522757e-01]
# [ 5.71679535e-01 4.18920646e-01 -8.41221046e-01 1.30117401e-01
# -3.29124434e-01 1.56958573e-01 5.67624154e-02]
# [-6.04968518e+00 -1.50079386e+00 3.32669570e-01 -8.32409733e-01
# 1.78232749e-01 5.24695454e-01 8.32773600e-02]
# [-1.15425643e+00 -1.30042860e+00 -2.49331018e-01 8.91249982e-02
# 4.54254674e-01 2.56823534e-01 -9.22042392e-02]
# [-3.81779109e+00 -1.24587414e+00 7.95713358e-01 1.93303778e-01
# -3.23974153e-01 -4.16896633e-01 -1.39692990e-01]
# [-5.40295361e-01 1.14273298e-02 -1.29265357e-01 -6.57905645e-01
# 6.45894517e-01 -3.78516483e-01 -1.68357775e-01]
# [-4.22840415e-01 -3.29803048e-01 -7.60598855e-01 -3.65240648e-02
# 2.07225597e-01 9.48671820e-02 -8.35446095e-02]
# [ 1.13697846e+00 3.84592953e-01 -6.87148612e-01 -3.24032790e-01
# 8.06800612e-02 -7.13744266e-02 -1.22936692e-01]
# [ 8.79393656e-01 -3.32501725e-01 -2.38244357e-01 -5.80770466e-02
# -6.38442276e-02 -1.17962650e-02 2.09756299e-02]
# [ 4.63059708e-01 -8.41578873e-01 2.14159505e-01 5.62059412e-01
# -1.89396606e-01 -8.86742022e-02 1.78835907e-02]
# [ 9.38052709e-01 -8.76801721e-01 -1.23954896e-01 4.42103302e-02
# 9.27107304e-02 -2.42431559e-01 -4.29193402e-02]
# [ 1.59391403e+00 -4.59780393e-01 -1.11242583e+00 -3.12285729e-01
# 6.87391982e-01 -1.22858441e-01 6.35740076e-02]
# [-4.12799855e-01 -1.33730330e+00 4.39706056e-01 3.16880737e-01
# -7.15960465e-01 -2.97983286e-01 1.24109366e-01]
# [ 1.27315658e+00 -6.00263134e-01 -3.06044076e-01 2.40780176e-02
# 2.82708039e-01 -1.96890176e-01 1.63601179e-01]
# [ 1.57885950e+00 3.73033694e-01 1.10821963e+00 1.25862086e+00
# 4.14286204e-01 1.48662984e-01 -1.62175425e-01]
# [ 1.63371071e+00 6.63251473e-01 1.69542719e+00 -8.15264815e-01
# 6.59124471e-01 -2.19789032e-01 1.90806339e-01]
# [ 4.93203443e-01 6.39993832e-02 1.74264440e-01 -2.02023669e-01
# -1.22144178e+00 5.91490692e-02 -2.29008993e-03]
# [ 1.50330439e+00 6.00191106e-01 5.50107647e-01 -3.21001007e-01
# -4.05458393e-01 -1.12808351e-02 -3.15578384e-02]
# [ 1.89646148e+00 3.30405307e-01 7.50920481e-01 6.40276442e-01
# 2.99591705e-01 4.36988233e-01 -2.90900627e-02]
# [ 1.29782027e+00 5.35103875e-01 -3.25805459e-01 1.07894923e-01
# 7.64133609e-02 -4.05772976e-02 8.03365159e-02]
# [ 1.25957713e+00 2.02783488e+00 2.43125738e-01 -9.82717479e-01
# -4.71177150e-01 2.19562277e-01 -4.34203379e-02]]
因子分析
因子分析（factor analysis）是主成分分析的推广和发展，也是多元统计分析中国呢降维分析的一种方法。主成分分析通过线性组合将多个原始变量综合成若干主成分。在多变量分析中，某些变量间往往存在相关性。那么，是什么原因使变量间有关联呢？是否窜在不能直接观测到的，但影响可观测变量变化的公共因子？
因子分析就是寻找这些公共因子的分析方法，它是在综合原始变量信息的基础上构筑若干意义较为明确的公因子，以它们为框架分解原始变量，以此考察原始变量间的联系与区别。
基本原理
因子分析的基本目的就是用少数几个公共因子去买哦数许多指标或因素之间的联系，即将比较密切的几个变量归在同一类中，每一类变量就称为一个因子(之所以称为因子，是因为它往往是不可观测的，类似于隐变量)，以较少的几个因子反映原始资料的大部分信息。
主成分分析是因子分析的一个特例。通常情况下可采用主成分法估算出因子个数，二者却别和联系主要体现在如下几个方面:
1.因子分析是把原始变量表示成各因子的线性组合，而主成分分析则是把主成分表示成各原始变量的线性组合
2.主成分分析的重点在于解释原始变量的总方差，而因子分析则把重点放在解释原始变量之间的协方差
3.因子分析中的因子个数可根据研究者的需要而事先指定，指定因子数量不同可导致分析结果不同。在主成分分析中，有几个变量就有几个主成分
4.主成分分析中，当给定的协方差矩阵或者相关矩阵的特征值是唯一的时候，主成分一般是唯一的；而因子分析中因子不是唯一的，可以旋转得到不同的因子
- 因子分析模型
因子分析是从研究变量内部相互依存关系出发，把一些具有错综复杂关系的变量归结为少数几个综合因子的一种多元统计分析方法。它的基本思想是将原始变量进行分类，将相关性较高，即联系比较紧密的变量分在同一类中，而不同类变量之间的相关性则较低，那么每一类变量实际上就代表了一个基本结果，即公共因子。对于所研究的问题就是试图用最少个数的不可观测的所谓公共因子的线性函数来描述所研究的对象。
因子分析的一般模型如下
\[\begin{align*}
& x_1= a_{11}F_1 + a_{12}F_2 + \cdots + a_{1m}F_m + \varepsilon_1 \\
& x_2= a_{21}F_1 + a_{22}F_2 + \cdots + a_{2m}F_m + \varepsilon_2 \\
& \vdots \\
& x_p= a_{p1}F_1 + a_{p2}F_2 + \cdots + a_{pm}F_m + \varepsilon_p
\end{align*}
\]
其中，\(F_j(j=1,2,3,\cdots,m)\) 表示不可观测的因子或公因子组成的向量。利用 \(\alpha\) 因子提取法、H arris成分分析法、主成分法、最大似然法等方法均可进行因子提取。因子的含义必须结合具体问题的实际意义而定。
\(a_{ij}(i=1,2,\cdots,p;j=1,2,\cdots,m)\) 称为因子载荷。因子载荷就是第 \(i\) 变量与第 \(j\) 因子的相关系数，反映了第 \(i\) 变量在第 \(j\) 因子上的重要性，即表示变量 \(x_i\) 依赖于 \(F_j\) 的份量（比重）。
在因子分析模型中，把原始变量 \(x_i\) 的信息能够被 \(m\) 个公因子解释的程度称作共同度。其计算公式如下
\[Communality = \sum_{i=1}^{p}\sum_{j=1}^{m} {a_{ij}^2}
\]
由此可以判断公因子的解释能力。
在实际问题中，究竟取所少个因子进行分析？可根据提取出来的主成分方差贡献率来决定，方差贡献率越大，因子分析越有意义。通常情况下可参考累积方差贡献率85%阈值进行判定。
- 因子旋转
因子分析的目的不仅是找出因子，更重要的是知道每个因子的意义，以便对实际问题进行分析。如果因子的典型代表变量不很突出，为了对公因子 \(F\) 能够更好地解释，还需要进行因子旋转，通过适当的旋转得到比较满意的主因子。即使得每个原始变量仅在一个公因子上有较大的载荷，而在其余的公因子上的载荷比较小。
进行因子旋转，就是要使因子载荷矩阵中因子载荷的平方值向0和1两个方向分化，使大的载荷更大，小的载荷更小。旋转的方法有很多，因子旋转过程中，按照旋转坐标轴的位置不同，如果主轴相互正交，则称为正交旋转；如果主轴相互间不是正交的，则称为斜交旋转。可供选择的因子旋转方法有：方差最大化法、四分位最大法、平衡法、正交旋转法等。一般实际问题中常用方法是方差最大化、正交旋转法
- 因子得分
在因子分析中，人们往往更愿意用公因子反映原始变量，这样有利于描述研究对象的特征。因而往往将因子表示为原始变量的线性组合，即因子得分函数：
\[\begin{align*}
& f_1= \beta_{11}x_1 + \beta_{12}x_2 + \cdots + \beta_{1p}x_p \\
& f_2= \beta_{21}x_1 + \beta_{22}x_2 + \cdots + \beta_{2p}x_p \\
& \vdots \\
& f_p= \beta_{m1}x_1 + \beta_{m2}x_2 + \cdots + \beta_{mp}x_p
\end{align*}
\]
因子得分函数可计算每个样本的因子得分。但因子得分函数中方程的个数 \(m\) 小于变量的个数 \(p\)，所以并不能精确计算出因子得分，只能对因子得分进行估计。估计因子得分的方法较多，有：回归分析法、Bartlett估计法、Thomson估计法等。
基本步骤和过程
因子分析的核心问题有两个：一是如何构造因子变量；二是如何对因子变量进行命名解释。通常情况下，在进行因子分析之前亦要进行标准化或无量纲化，然后 可按照如下顺序进行因子分析
1.考察原始变量之间的相关性，如果各变量之间是独立的，那么可能不适用因子分析
2.计算变量之间的相关系数矩阵作为分析基础
3.确定提取公因子的方法并根据累积方差贡献率阈值进一步确认提取公共因子的数目
4.进行因子旋转，使公因子更具有解释性
5.计算各公因子得分
6.可根据各公因子得分对各样本进行综合考察
示例
import pandas as pd
import statsmodels.api as sm
import scipy.stats as stats
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
# 利用因子分析法对该网吧满意度进行综合评价
ic = pd.read_csv(\'./data/internet_cafe.csv\')
print(ic.head(5))
# No Switch Connection Speed ... Standard Settlement Success Efficiency
# 0 24 8 8 8 ... 2 2 2 2
# 1 14 3 10 4 ... 6 5 6 5
# 2 40 7 6 5 ... 5 6 5 5
# 3 25 4 5 2 ... 2 4 6 5
# 4 64 3 5 3 ... 8 8 5 6
# 14个变量对满意度进行分析，有些复杂，需提取公共因子。
# 构造用于分析的原始数据
X = ic.iloc[:, 1:len(ic.columns)] # 构造从Switch到Efficiency组成的分析对象
# 因子分析
from sklearn.decomposition import FactorAnalysis as skFA
x1 = (X - X.mean()) / X.std(ddof=1) # 对数据进行标准化
ic_fa = skFA(n_components=2).fit(x1)
# res = skFA(n_components=2).fit_transform(X) # 可以直接得到因子得分结果
# print(res)
# 查看所提取因子的因子载荷
ic_ev = pd.DataFrame(ic_fa.components_, index=[\'Factor1\', \'Factor2\'],
columns=(ic.columns)[1:len(ic.columns)]).T
print(ic_ev)
# Factor1 Factor2
# Switch -0.582475 0.528556
# Connection -0.429028 0.454938
# Speed -0.485041 0.692095
# Transformation -0.597188 0.726102
# Offline -0.446861 0.350302
# Timeliness -0.774073 -0.195317
# Initiative -0.720337 -0.341741
# Attitude -0.800199 -0.191607
# Skill -0.877507 -0.134740
# Consideration -0.713259 -0.313348
# Standard -0.818312 -0.336150
# Settlement -0.938716 -0.153299
# Success -0.912337 -0.095773
# Efficiency -0.903344 -0.052725
# 与Factor1相关性强的变量为Timeliness,Initiative,Attitude,Skill,Settlement,Success,Efficiency
# 与Factor2相关性强的变量为Switch,Connection,Speed,Transformation,Offline,Initiative,Consideration,Standard
# Factor1可命名为"人员服务质量"因子，Factor2可命名为"计数质量"因子
# 部分原始变量如Swith/Connection在两个因子的载荷差异较小，给因子含义确定带来困难，需要使用因子旋转的方法。
# 展示特征根及其对应贡献
def PCA(x, components=None):
if components == None:
components = x.size / len(x) # 如果components参数未指定，则赋值为原始变量个数
average = np.mean(x, axis=0)
sigma = np.std(x, axis=0, ddof=1)
r, c = np.shape(x)
data_standardized = []
mu = np.tile(average, (r, 1))
data_standardized = (x - mu) / sigma # 使用Z-Score法对原始数据标准化
cov_matrix = np.cov(data_standardized.T) # 计算协方差矩阵
EigenValue, EigenVector = np.linalg.eig(cov_matrix) # 求解协方差矩阵的特征值和特征向量
index = np.argsort(-EigenValue) # 按照特征值大小降序排列
Score = []
Selected_Vector = EigenVector.T[index[:int(components)]]
Score = data_standardized * np.matrix(Selected_Vector.T)
return EigenValue[index], Selected_Vector, np.array(Score)
EigenValue, Vector, Score = PCA(np.asarray(x1))
x1_ev = pd.DataFrame((EigenValue), columns=[\'EigenValue\'], index=list(range(1, 15)))
prop = x1_ev[\'EigenValue\'] / x1_ev[\'EigenValue\'].sum()
s = 0
p, c = [], []
for i in range(1, len(prop)+1):
s += prop[i]
p.append(prop[i])
c.append(s)
x1_ev[\'Proportion\'] = p
x1_ev[\'Cumulative\'] = c
print(x1_ev)
# EigenValue Proportion Cumulative
# 1 7.953348 0.568096 0.568096
# 2 2.443283 0.174520 0.742617
# 3 0.785253 0.056090 0.798706
# 4 0.652081 0.046577 0.845283
# 5 0.495501 0.035393 0.880676
# 6 0.427167 0.030512 0.911188
# 7 0.332092 0.023721 0.934909
# 8 0.243358 0.017383 0.952292
# 9 0.182532 0.013038 0.965330
# 10 0.143056 0.010218 0.975548
# 11 0.126175 0.009013 0.984560
# 12 0.104105 0.007436 0.991996
# 13 0.059469 0.004248 0.996244
# 14 0.052581 0.003756 1.000000
# 从各特征根贡献来看，前两个与其他的差别较大，且累计达到了74%,基本上可以较大程度上反映了原始数据的信息
# 因子得分
ic[\'Factor1\'] = ic_fa.transform(x1)[:, 0]
ic[\'Factor2\'] = ic_fa.transform(x1)[:, 1]
print(ic)
# No Switch Connection Speed ... Success Efficiency Factor1 Factor2
# 0 24 8 8 8 ... 2 2 3.527966 3.193777
# 1 14 3 10 4 ... 6 5 1.829194 -0.982942
# 2 40 7 6 5 ... 5 5 1.686458 0.433412
# 3 25 4 5 2 ... 6 5 2.385283 -0.346508
# 4 64 3 5 3 ... 5 6 1.248211 -1.841660
# .. .. ... ... ... ... ... ... ... ...
# 65 15 9 8 7 ... 10 10 -0.993967 0.210845
# 66 36 8 8 7 ... 10 10 -0.922337 -0.495290
# 67 48 9 9 8 ... 10 10 -1.071304 0.353790
# 68 30 8 8 7 ... 7 10 -0.208388 -0.073899
# 69 39 9 7 8 ... 10 10 -0.841607 0.923970
python中目前尚无成熟的用于因子旋转的包或模块。此外，python中现有计算载荷方法均采用SVD，与传统统计分析过程得到的结果或结论有差异。基于传统因子分析流程，定义了如下类FA
import numpy as np
import pandas as pd
# 基于方差最大化正交旋转的因子分析
class FA(object):
"""
该类用于因子分析，有5个方法分别用于计算：特征值和方差贡献率、旋转前的因子载荷、旋转后的因子载荷、因子得分系数和因子得分
"""
def __init__(self, component, gamma=1.0, q=20, tol=1e-8):
"""
Args:
component: 提取公因子的数目，在做因子分析时需事先指定
gamma: 最大方差正价参数
q: 最大方差正价参数
tol: 最大方差正价参数
"""
self.component = component
self.gamma = gamma
self.q = q
self.tol = tol
def var_contribution(self, data):
"""
该方法用于输出特征值、方差贡献率及累计方差贡献率
"""
# 将数据转存为数组形式方便操作
var_name = data.columns
data = np.array(data)
# 标准化数据
z = (data - data.mean()) / data.std(ddof=1)
# 按列求解相关系数矩阵，存入cor中(rowvar=0指定按列求解相关系数矩阵)
cor = np.corrcoef(z, rowvar=0)
# 求解相关系数矩阵特征值与特征向量，并按照特征值由大到小排序
# 注意：numpy中求出的特征向量是按列排列而非行，需要转置
eigvalue, eigvector = np.linalg.eig(cor)
eigdf = pd.DataFrame(eigvector.T).join(pd.DataFrame(eigvalue, dtype=float, columns=[\'eigvalue\']))
# 将特征向量按照特征值由大到小排列
eigdf = eigdf.sort_values(by=\'eigvalue\', ascending=False)
# 将调整好的特征向量存储在eigvector
eigvector = np.array(eigdf.iloc[:, :-1])
# 将特征值由大到小排序，存入eigvalue
eigvalue = list(np.array(eigvalue, dtype=float))
eigvalue.sort(reverse=True)
# 计算每个特征值的方差贡献率，存入varcontribution中
varcontribution = list(np.array(eigvalue / sum(eigvalue), dtype=float))
# 累计方差贡献率
leiji_varcontribution = []
for i in range(len(varcontribution)):
s = float(sum(varcontribution[: i + 1]))
leiji_varcontribution.append(s)
# 将特征值、方差贡献率、累计方差贡献率写入DataFrame
# 控制列的输出顺序
col = ["Eigvalue", "Proportion", "Cumulative"]
eig_df = pd.DataFrame(
{"Eigvalue": eigvalue, "Proportion": varcontribution, "Cumulative": leiji_varcontribution}, columns=col)
self.eigvalue = eigvalue
self.eigvector = eigvector
return eig_df
def loadings(self, data):
"""
用于输出旋转前的因子载荷阵
"""
factor_num = self.component
# 求解因子载荷阵
# 生成由前factor_num个特征值构成的对角阵，存入duijiao中用于计算因子载荷阵
eigvalue = self.var_contribution(data)["Eigvalue"]
duijiao = list(np.array(np.sqrt(eigvalue[:factor_num]), dtype=float))
eigmat = np.diag(duijiao)
zaihe = np.dot(self.eigvector[:factor_num].T, eigmat)
self.zaihe = zaihe
n = range(1, factor_num + 1)
col = []
for i in n:
c = \'Factor \' + str(i)
col.append(c)
zaihe = -pd.DataFrame(zaihe, columns=col)
zaihe.iloc[:, 1] = -zaihe.iloc[:, 1]
self.col = col
zaihe.index = data.columns
self.zaihe = zaihe
return zaihe
def varimax_rotation(self, data):
"""
对因子载荷阵进行最大方差正交旋转，返回旋转后的因子载荷阵
"""
zaihe = self.loadings(data)
m, n = zaihe.shape
R = np.eye(n)
d = 0
for i in range(self.q):
d_init = d
Lambda = np.dot(zaihe, R)
w, a, wa = np.linalg.svd(
np.dot(zaihe.T, np.asarray(Lambda) ** 3 - (self.gamma / m) * np.dot(Lambda, np.diag(
np.diag(np.dot(Lambda.T, Lambda))))))
R = np.dot(w, wa)
d = np.sum(a)
if d_init != 0 and d / d_init < self.tol:
break
orthogonal = np.dot(zaihe, R)
self.orthogonal = orthogonal
return pd.DataFrame(orthogonal, index=data.columns, columns=self.col)
def score_coef(self, data):
"""
用于计算因子得分函数
"""
# R为原始变量的相关阵
corr = np.corrcoef(data, rowvar=0)
A = self.varimax_rotation(data)
coefficient = pd.DataFrame(np.dot(np.array(A).T, np.mat(corr).I), columns=data.columns, index=self.col)
self.coefficient = coefficient
return coefficient
def score(self, data):
"""
用于计算因子得分
"""
data_scale = (data - data.mean()) / data.std(ddof=1)
F = np.dot(data_scale, self.coefficient.T)
F = pd.DataFrame(F)
col2 = []
n = range(1, self.component + 1)
for i in n:
c = "Score F" + str(i)
col2.append(c)
F.columns = col2
return F
if __name__ == \'__main__\':
ic = pd.read_csv(\'./data/internet_cafe.csv\')
print(ic.head(5))
# No Switch Connection Speed ... Standard Settlement Success Efficiency
# 0 24 8 8 8 ... 2 2 2 2
# 1 14 3 10 4 ... 6 5 6 5
# 2 40 7 6 5 ... 5 6 5 5
# 3 25 4 5 2 ... 2 4 6 5
# 4 64 3 5 3 ... 8 8 5 6
# 1. 构造原始数据
X = ic.iloc[:, 1:len(ic.columns)]
# 类实例化
ic_fa = FA(component=2)
# # 2. 变量特征根及对应的方差贡献率和累计方差贡献率
contribution = ic_fa.var_contribution(X)
# print(contribution)
# # Eigvalue Proportion Cumulative
# # 0 7.953348 0.568096 0.568096
# # 1 2.443283 0.174520 0.742617
# # 2 0.785253 0.056090 0.798706
# # 3 0.652081 0.046577 0.845283
# # 4 0.495501 0.035393 0.880676
# # 5 0.427167 0.030512 0.911188
# # 6 0.332092 0.023721 0.934909
# # 7 0.243358 0.017383 0.952292
# # 8 0.182532 0.013038 0.965330
# # 9 0.143056 0.010218 0.975548
# # 10 0.126175 0.009013 0.984560
# # 11 0.104105 0.007436 0.991996
# # 12 0.059469 0.004248 0.996244
# # 13 0.052581 0.003756 1.000000
# # 前两个特征根与其他相比，贡献率较大，且总和达到74%，基本上可以在较大程度上反映原始数据的信息
# # 3.提取2个因子，前2个特征根对应的特征向量为
loading = ic_fa.loadings(X)
# print(loading)
# # Factor 1 Factor 2
# # Switch 0.612655 0.629533
# # Connection 0.473459 0.549148
# # Speed 0.494931 0.677657
# # Transformation 0.595114 0.691249
# # Offline 0.508545 0.501489
# # Timeliness 0.816996 -0.239573
# # Initiative 0.751754 -0.385172
# # Attitude 0.833623 -0.185283
# # Skill 0.883247 -0.151843
# # Consideration 0.751921 -0.341353
# # Standard 0.824808 -0.365608
# # Settlement 0.932066 -0.163706
# # Success 0.909980 -0.109911
# # Efficiency 0.919745 -0.068276
# # 4.为了提高公因子的解释性，进行因子旋转
rotated_loadings = ic_fa.varimax_rotation(X)
# print(rotated_loadings)
# # Factor 1 Factor 2
# # Switch 0.254884 0.840650
# # Connection 0.168188 0.705294
# # Speed 0.128211 0.829299
# # Transformation 0.210946 0.887405
# # Offline 0.221250 0.679086
# # Timeliness 0.835723 0.162619
# # Initiative 0.844677 0.003322
# # Attitude 0.825545 0.218478
# # Skill 0.854255 0.270980
# # Consideration 0.824691 0.042318
# # Standard 0.900574 0.054266
# # Settlement 0.903065 0.282876
# # Success 0.858731 0.320507
# # Efficiency 0.848273 0.361974
# # 旋转后，各变量在两个因子上的载荷数值相差变大，解释性更好
# # 与Factor 1强相关的变量有：Timeliness,Initiative,Attitude,Skill,Consideration,Standard,Settlement,Success,Efficiency
# # 与Factor 2强相关的变量有：Switch,Connection,Speed,Transformation,Offline
# # 故Factor 1命名为"人员服务质量"因子；Factor 2命名为"计数质量"因子
# # 5.因子得分系数
coef = ic_fa.score_coef(X)
print(coef.T)
# Factor 1 Factor 2
# Switch -0.049975 0.264243
# Connection -0.050402 0.226980
# Speed -0.072171 0.274936
# Transformation -0.063540 0.285664
# Offline -0.037521 0.211682
# Timeliness 0.136292 -0.039889
# Initiative 0.156388 -0.096586
# Attitude 0.127939 -0.019193
# Skill 0.127192 -0.004170
# Consideration 0.148166 -0.080648
# Standard 0.160867 -0.085254
# Settlement 0.134875 -0.005662
# Success 0.122291 0.012618
# Efficiency 0.115552 0.028317
# # 6.因子得分
s = ic_fa.score(X)
s[\'No\'] = ic[\'No\']
res = s[s[\'Score F1\'] == s[\'Score F1\'].min()]
print("Cafe No.{0} is the min of Factor1: {1}".format(
s[s[\'Score F1\'] == s[\'Score F1\'].min()][\'No\'].values[0], s[\'Score F1\'].min()))
print("Cafe No.{0} is the min of Factor2: {1}".format(
s[s[\'Score F2\'] == s[\'Score F2\'].min()][\'No\'].values[0], s[\'Score F2\'].min()))
# Cafe No.24 is the min of Factor1: -4.72029769621459
# Cafe No.38 is the min of Factor2: -3.3789552361942445
# 编号为24的网吧在"人员服务质量"上得分最低，为-4.72分; 编号为38的网吧在"技术质量"上得分最低，为-3.37
由于因子得分可按照多种方法进行旋转，在估计因子得分函数的系数过程中也可使用多种分析方法，故因子得分函数的系数是不确定的。因此在实际问题中，一般不再像主成分分析过程中那样可根据各公因子得分及其贡献计算每个样本的综合得分。