前言:
因为spss不能直接得到主成分得分系数,参考csdn上其他博主写的文章,整理了一下用于计算主成分得分系数的代码
主成分分析原理
先略,后面再补
主成分分析代码
需要用到的库及文件读取,以下以读取csv文件为例,pandas还可以读取excel、sav(spss常用的数据集格式)等格式
案例数据:
全国重点水泥企业某年的经济效益分析,评价指标有:X1为固定资产利税率, X2为资金利税率, X3为销售收入利税率, X4为资金利润率, X5为固定资产产值率, X6-流动资金周转天数, X7-万元产值能耗, X8-全员劳动生产率;现有15家水泥企业的数据如下,试利用主成分法综合评价其效益。
| X1 | X2 | X3 | X4 | X5 | X6 | X7 | X8 |
企业A | 16.68 | 26.75 | 31.84 | 18.4 | 53.25 | 55 | 28.83 | 1.75 |
企业B | 19.7 | 27.56 | 32.94 | 19.2 | 59.82 | 55 | 32.92 | 2.87 |
企业C | 15.2 | 23.4 | 32.98 | 16.24 | 46.78 | 65 | 41.69 | 1.53 |
企业D | 7.29 | 8.97 | 21.3 | 4.76 | 34.39 | 62 | 39.28 | 1.63 |
企业E | 29.45 | 56.49 | 40.74 | 43.68 | 75.32 | 69 | 26.68 | 2.14 |
企业F | 32.93 | 42.78 | 47.98 | 33.87 | 66.46 | 50 | 32.87 | 2.6 |
企业G | 25.39 | 37.85 | 36.76 | 27.56 | 68.18 | 63 | 35.79 | 2.43 |
企业H | 15.05 | 19.49 | 27.21 | 14.21 | 56.13 | 76 | 35.76 | 1.75 |
企业I | 19.82 | 28.78 | 33.41 | 20.17 | 59.25 | 71 | 39.13 | 1.83 |
企业J | 21.13 | 35.2 | 39.16 | 26.52 | 52.47 | 62 | 35.08 | 1.73 |
企业K | 16.75 | 28.72 | 29.62 | 19.23 | 55.76 | 58 | 30.08 | 1.52 |
企业L | 15.83 | 28.03 | 26.4 | 17.43 | 61.19 | 61 | 32.75 | 1.6 |
企业M | 16.53 | 29.73 | 32.49 | 20.63 | 50.14 | 69 | 37.57 | 1.31 |
企业N | 22.24 | 54.59 | 31.05 | 37 | 67.95 | 63 | 32.33 | 1.57 |
企业O | 12.92 | 20.82 | 25.12 | 12.54 | 51.07 | 66 | 39.18 | 1.83 |
1. 库及数据集文件读取
# 数据处理
import pandas as pd
import numpy as np
# 绘图
import seaborn as sns
import matplotlib.pyplot as plt
#读取数据集文件
df = pd.read_csv("hw91.csv")
'''怎么把数据变成csv格式?
pd.read_csv默认读取UTF-8的编码模式
可以把数据复制到excel选择另存为CSV UTF-8(逗号分隔)格式'''
2. 相关性检验
进行主成分分析前,需要对数据进行相关性检验
检验方法一般是
(1)KMO检验
(2)巴特利球形检验
# KMO检验
''' 检查变量间的相关性和偏相关性,取值在0-1之间;
KMO统计量越接近1,变量间的相关性越强,偏相关性越弱,因子分析的效果越好
'''
#记得检查一下是否安装了 factor_analyzer的库
from factor_analyzer.factor_analyzer import calculate_kmo
#kmo_all, kmo_model = calculate_kmo(df)
kmo = calculate_kmo(df)
print(f'KMO检验结果为:{kmo[1]}')
#进行球状检验
from factor_analyzer.factor_analyzer import calculate_bartlett_sphericity
chi_square_value, p_value = calculate_bartlett_sphericity(df)
print(f'卡方={chi_square_value},P值={p_value}')
输出:
附:spss结果
(3)数据标准化,求相关系数矩阵
# 数据标准化
from sklearn import preprocessing
dfz = preprocessing.scale(df)
print(f'数据标准化后结果如下:\n{dfz}')
数据标准化结果
利用标准化后的数据求相关系数矩阵
# 相关系数矩阵
covX = np.around(np.corrcoef(dfz.T),decimals=3)
print(f'相关系数矩阵如下:\n{covX}')
输出结果:
附:spss结果
3. 主成分选择
(1)特征值和特征向量的求解
# 求解特征值和特征向量
featValue, featVec = np.linalg.eig(covX.T)
print(f'解得,特征值为:\n{featValue}, \n特征向量为:\n{featVec}')
featValue = sorted(featValue)[::-1]
featValue
输出结果:
绘制散点图和折线图
# 绘制散点图和折线图
plt.scatter(range(1, df.shape[1] + 1),featValue)
plt.plot(range(1, df.shape[1]+1), featValue)
#显示图的标题和xy轴的名字
plt.title('Scree Plot')
plt.xlabel('Factors')
plt.ylabel('Eigenvalue')
plt.grid()
plt.show()
输出散点图(也就是常说的碎石图)
进行特征值贡献度和累计贡献度的计算
#特征值贡献度
gx = featValue/np.sum(featValue)
#累计贡献度
lg = np.cumsum(gx)
print(f'特征值累计贡献度:\n{lg}')
输出结果
这里可以加一个绘制累计贡献度图像的代码,下回再加()
选择主成分:
关于主成分的选择,一般根据特征值累计贡献度或者根据特征值是否大于1,本例参考spss的评判标准,选择用特征值是否大于1
#主成分选择
'''参考spss判定:特征值>1的选为主成分'''
k = [i for i in range(0, len(featValue)) if featValue[i]>1]
k = list(k)
print(f'选择的主成分为:{k}')
'''第0个成分即第1个成分'''
输出结果:
附:spss 主成分选择的结果
关于主成分选择地判定条件当然还可以这样选择:
'''下行代码为按照贡献度进行判断,但哪种才是最通用的方法,目前我也不清楚'''
k = [i for i in range (len(lg)) if lg[i] < 0.85]
'''下行代码适用于主成分分析回归不降维时的主成分得分系数计算
也就是有多少个自变量就有多少个主成分,这样能够保证数据最大程度地被提取'''
k= [i for i in range(0, len(featValue))]
4. 主成分得分系数计算
选择主成分对应的特征向量矩阵,即为主成分得分系数(?)与spss得到的结果一致
#选择出主成分对应的特征向量矩阵
selectVec = np.matrix(featVec.T[k]).T
selectVec = selectVec*(-1)
print(f'主成分对应的特征向量矩阵为:\n{selectVec}')
以下为spss计算得到的结果:但此处我们求出的是公因子的得分而非主成分得分,还需将因子得分乘以相应的方差才可得到主成分得分,计算结果如下表所示。
可以发现python得到的结果和spss处理后的结果一致