【应用多元统计分析】-王学民Python主成分分析例题，特征值处理和可视化（2）
2024-05-08 07:45:24
title: “应用多元统计分析”
subtitle: “书上题目”
author: | OLSRR
由于字数限制，本文省去部分数据预览。
7.6
下表中给出的是美国 50 个州每 100000 个人中七种犯罪的比率数据, 这七种犯罪是: 杀人罪 ( x 1 ) \left(x_{1}\right) (x1) 、 强奸罪 ( x 2 ) \left(x_{2}\right) (x2) 、抢劫罪 ( x 3 ) \left(x_{3}\right) (x3) 、伤害罪 ( x 4 ) \left(x_{4}\right) (x4) 、夜盗罪 ( x 5 ) \left(x_{5}\right) (x5) 、盗穷罪 ( x 6 ) \left(x_{6}\right) (x6) 和汽车犯罪 ( x 7 ) \left(x_{7}\right) (x7) 。试对表中犯罪数据进行主成分分析。
答案
数据库准备：
import pandas as pd#不一定都会用到
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import matplotlib.pyplot as mp,seaborn
from pylab import mpl
from sklearn import preprocessing
from factor_analyzer.factor_analyzer import calculate_bartlett_sphericity
from factor_analyzer.factor_analyzer import calculate_kmo
数据准备：
df = pd.read_csv("/Users/che/Desktop/应用多元分析/《应用多元统计分析》（第五版，王学民 编著）配书资料/《应用多元统计分析》（第五版）文本数据（以逗号为间隔）/exec7.6.csv", encoding='gbk',index_col=0).reset_index(drop=True)
print(df)#显示数据
# x1 x2 x3 x4 x5 x6 x7
#0 14.2 25.2 96.8 278.3 1135.5 1881.9 280.7
#1 10.8 51.6 96.8 284.0 1331.7 3369.8 753.3
……
#49 5.4 21.9 39.7 173.9 811.6 2772.2 282.0
Bartlett’s球状检验:
chi_square_value, p_value = calculate_bartlett_sphericity(df)
print(chi_square_value, p_value)
#225.98725841755999 2.5997375924269063e-36
KMO检验：
kmo_all, kmo_model = calculate_kmo(df)
#[0.66737389 0.83457371 0.84650891 0.86291368 0.78956368 0.66450087
# 0.72385252]
检查了变量间的相关性和偏相关性，由此可以看出，数据均大于0.5，变量间的相关性较强，同时偏相关性较弱，由此进行因子分析的效果较好。
标准化：
df = preprocessing.scale(df)
print(df)
[[ 1.76493364e+00 -5.01338300e-02 -3.12049014e-01 6.75093885e-01-3.65336599e-01 -1.09848840e+00 -5.05748984e-01]……[-5.33973411e-01 -3.59949633e-01 -9.64914274e-01 -3.76843452e-01-1.12191907e+00 1.40426079e-01 -4.98958725e-01]]
计算相关系数:
e76_1=df.iloc[0:50,0:8]
df76 = pd.DataFrame(e76_1)
df76_corr = df76.corr()
print(df76_corr)
# x1 x2 x3 x4 x5 x6 x7
#x1 1.000000 0.601220 0.483708 0.648550 0.385817 0.101920 0.068814
……
#x7 0.068814 0.348902 0.590680 0.275843 0.557953 0.444180 1.000000
可视化结果：
seaborn.heatmap(df76_corr, center=0, annot=True, cmap='YlGnBu')
mp.show()
该相关矩阵表明，变量之间存在一定的相关性，即彼此之间信息有不少是重复的，从而有一定的降维空间。
特征值相关：
eig76_value,eig76_vector=np.linalg.eig(df76_corr)
df=pd.DataFrame({"eig76_value":eig76_value})
df=df.sort_values(by=["eig76_value"], ascending=False) # 获取累积贡献度
df["eig76_cum"] = (df["eig76_value"]/df["eig76_value"].sum()).cumsum()
eig76=df.merge(pd.DataFrame(eig76_vector).T, left_index=True, right_index=True)
print(eig76)
# eig76_value eig76_cum 0 ... 4 5 6
#0 4.114960 0.587851 0.300279 ... 0.440157 0.357360 0.295177
#1 1.238722 0.764812 0.629174 ... -0.203341 -0.402319 -0.502421
#2 0.725817 0.868500 0.178245 ... -0.209895 -0.539231 0.568384
#4 0.316432 0.913704 -0.232114 ... -0.057555 -0.234890 0.419238
#5 0.257974 0.950558 0.538123 ... 0.101033 0.030099 0.369753
#6 0.222039 0.982278 -0.259117 ... -0.535987 -0.039406 0.057298
#3 0.124056 1.000000 -0.267593 ... 0.648117 -0.601690 -0.147046
可视化：
plt.scatter(range(1, df76.shape[1] + 1), eig76_value)
plt.plot(range(1, df76.shape[1] + 1), eig76_value)
plt.title("Scree Plot")
plt.xlabel("Factors")
plt.ylabel("Eigenvalue")
plt.grid() # 显示网格
plt.show() # 显示图形
数据处理：
df76_std = (df76 - df76.mean())/df76.std()
loading = eig76.iloc[:2,2:].T
loading["vars"]=df76_std.columns
print(loading)
# 0 1 vars
#0 0.300279 0.629174 x1
……
#6 0.295177 -0.502421 x7
score = pd.DataFrame(np.dot(df76_std,loading.iloc[:,0:2]))
print(score)
0 1
0 -0.049880 2.096102
1 2.421515 -0.166523
……
49 -1.424635 -0.062683
可以认为，第一主成分是对所有犯罪率的度量，第二主成分是用于度量暴力犯罪在犯罪性质上占的比重，第三主成分很难给出明显的解释，因此只取前两个主成分。
mpl.rcParams['font.sans-serif'] = ['FangSong'] # 解决保存图像是负号'-'显示为方块的问题
mpl.rcParams['axes.unicode_minus'] = False
plt.plot(loading[0],loading[1], "o")
xmin ,xmax = loading[0].min(), loading[0].max()
ymin, ymax = loading[1].min(), loading[1].max()
dx = (xmax - xmin) * 0.2
dy = (ymax - ymin) * 0.2
plt.xlim(xmin - dx, xmax + dx)
plt.ylim(ymin - dy, ymax + dy)
plt.xlabel('first')
plt.ylabel('second')
for x, y,z in zip(loading[0], loading[1], loading["vars"]):
plt.text(x, y+0.1, z, ha='center', va='bottom', fontsize=13)
plt.grid(True)
plt.show()#用绝对值做比较
plt.plot(score[0],score[1], "o")
xmin ,xmax = score[0].min(), score[0].max()
ymin, ymax = score[1].min(), score[1].max()
dx = (xmax - xmin) * 0.2
dy = (ymax - ymin) * 0.2
plt.xlim(xmin - dx, xmax + dx)
plt.xlabel('first')
plt.ylabel('second')
for x, y,z in zip(score[0], score[1], score.index):plt.text(x, y+0.1, z, ha='center', va='bottom', fontsize=13)
plt.grid(True)
plt.show()
7.7
下表(完整数据可从作者网页上下载) 是纽约股票交易所的五只股票(阿莱德化学、杜邦、联合碳化物、埃克森和德士古)从 1975 年 1 月到 1976 年 12 月期间的周回报率。周回报率定义为
周回报率 = 本周五收盘价 − 上周五收盘价 上周五收盘价 \text{周回报率}=\frac{\text{本周五收盘价}-\text{上周五收盘价}}{\text{上周五收盘价}} 周回报率=上周五收盘价本周五收盘价−上周五收盘价
有拆股和支付股息时对收盘价进行调整, 试作主成分分析。
注:阿莱德化学、杜邦和联合碳化物属于化工类股票, 埃克森和德士古属于石油类股票。
答案
数据准备：
#不一定每个库都用到
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import matplotlib.pyplot as mp,seaborn
from factor_analyzer import calculate_bartlett_sphericity, calculate_kmo
from pylab import mpl
from sklearn import preprocessing
读取数据：
df1 = pd.read_csv("/Users/zz/Desktop/应用多元分析/《应用多元统计分析》（第五版，王学民 编著）配书资料/《应用多元统计分析》（第五版）文本数据（以逗号为间隔）/exec7.7.csv", encoding='gbk',index_col=0)
print(df1)
# x1 x2 x3 x4 x5#0.000000 0.000000 0.000000 0.039473 0.000000
#0.027027 -0.044855 -0.003030 -0.014466 0.043478
#0.122807 0.060773 0.088146 0.086238 0.078124
#0.057031 0.029948 0.066808 0.013513 0.019512
#0.063670 -0.003793 -0.039788 -0.018644 -0.024154
#... ... ... ... ...
#0.000000 -0.020080 -0.006579 0.029925 -0.004807
#0.021429 0.049180 0.006622 -0.002421 0.028985
#0.045454 0.046375 0.074561 0.014563 0.018779
#0.050167 0.036380 0.004082 -0.011961 0.009216
#0.019108 -0.033303 0.008362 0.033898 0.004566
#
#[100 rows x 4 columns]
Bartlett’s球状检验:
chi_square_value, p_value = calculate_bartlett_sphericity(df1)
print(chi_square_value, p_value)
#108.0670240489058 5.175016294414447e-21
标准化：
df1 = preprocessing.scale(df1)
print(df1)
[[-0.13526816 -0.13836418 -0.14405774 1.17734079 -0.13531234]……[ 0.34041039 -1.09296855 0.0689901 0.97952993 0.03128677]]
相关系数：
df77 = pd.DataFrame(df1)
# 计算相关系数
df77_corr = df77.corr()
print(df77_corr)
# x1 x2 x3 x4 x5
# x1 1.000000 0.576924 0.508656 0.386721 0.462178
……
# x5 0.462178 0.321953 0.425627 0.523529 1.000000
可视化：
seaborn.heatmap(df77_corr, center=0, annot=True, cmap='YlGnBu')
mp.show()
该相关矩阵表明，变量之间存在一定的相关性，即彼此之间信息有不少是重复的，从而有一定的降维空间。
特征值处理：
eig77_value,eig77_vector=np.linalg.eig(df77_corr)
eig77=pd.DataFrame({"eig77_value":eig77_value})
eig77=eig77.sort_values(by=["eig77_value"], ascending=False) # 获取累积贡献度
eig77["eig77_cum"] = (eig77["eig77_value"]/eig77["eig77_value"].sum()).cumsum()
eig77=eig77.merge(pd.DataFrame(eig77_vector).T, left_index=True, right_index=True)
print(eig77)
# eig77_value eig77_cum 0 1 2 3 4
#0 2.856487 0.571297 0.463541 0.457076 0.469980 0.421677 0.421329
#1 0.809118 0.733121 0.240850 0.509100 0.260577 -0.525265 -0.582242
#3 0.540044 0.841130 0.613357 -0.177900 -0.337036 -0.539018 0.433603
#4 0.451347 0.931399 -0.381373 -0.211307 0.664098 -0.472804 0.381227
#2 0.343004 1.000000 0.453288 -0.674981 0.395725 0.179448 -0.387467
画图：
plt.scatter(range(1, df77.shape[1] + 1), eig77_value)
plt.plot(range(1, df77.shape[1] + 1), eig77_value)
plt.title("Scree Plot")
plt.xlabel("Factors")
plt.ylabel("Eigenvalue")
plt.grid()
plt.show()
得分情况：
df77_std = (df77 - df77.mean())/df77.std()
loading = eig77.iloc[:2,2:].T
loading["vars"]=df77_std.columns
print(loading)
score = pd.DataFrame(np.dot(df77_std,loading.iloc[:,0:2]))
print(score)
# 0 1 vars
#0 0.463541 0.240850 0
#1 0.457076 0.509100 1
#2 0.469980 0.260577 2
#3 0.421677 -0.525265 3
#4 0.421329 -0.582242 4
0 1
0 0.244565 -0.676780
1 -0.203899 -1.105631
.. ... ...
99 0.116289 -0.984235
可以认为，第一主成分是对所有犯罪率的度量，第二主成分是用于度量暴力犯罪在犯罪性质上占的比重，第三主成分很难给出明显的解释，因此只取前两个主成分。
可视化：
mpl.rcParams['font.sans-serif'] = ['FangSong']
mpl.rcParams['axes.unicode_minus'] = False
plt.plot(loading[0],loading[1], "o")
xmin ,xmax = loading[0].min(), loading[0].max()
ymin, ymax = loading[1].min(), loading[1].max()
dx = (xmax - xmin) * 0.2
dy = (ymax - ymin) * 0.2
plt.xlim(xmin - dx, xmax + dx)
plt.ylim(ymin - dy, ymax + dy)
plt.xlabel('first')
plt.ylabel('second')
for x, y,z in zip(loading[0], loading[1], loading["vars"]):plt.text(x, y+0.1, z, ha='center', va='bottom', fontsize=13)
plt.grid(True)
plt.show()
plt.plot(score[0],score[1], "o")
xmin ,xmax = score[0].min(), score[0].max()
ymin, ymax = score[1].min(), score[1].max()dx = (xmax - xmin) * 0.2
dy = (ymax - ymin) * 0.2plt.xlim(xmin - dx, xmax + dx)
plt.ylim(ymin - dy, ymax + dy)plt.xlabel('first')
plt.ylabel('second')
for x, y,z in zip(score[0], score[1], score.index):plt.text(x, y+0.1, z, ha='center', va='bottom', fontsize=13)plt.grid(True)
plt.show()
【应用多元统计分析】-王学民Python主成分分析例题，特征值处理和可视化（2）相关推荐
- 【应用多元统计分析】上机四五——主成分分析因子分析
目录 一.主成分分析 1.princomp命令 2.screeplot命令 3.[例7.3.3]对[例6.3.3]中的数据从相关矩阵出发进行主成分分析 编辑(1)代码 (2)碎石图 (3)散点图 二 ...
- 多元统计分析—R与Python的实现配套代码
第2章 R和Python代码 set.seed(1010) (A=matrix(rnorm(12,-1,1),3,4)) #赋值表达式加括号意味着打印结果 (B=matrix(rnorm(12,1,1 ...
- 多元统计分析 小总结 python实现
元:数学上称变量为元 一元统计分析:一个变量(标志值或指标)的统计分析理论和方法(单个随机变量
- 【多元统计分析】因子分析——SPSS上机实验（过程+结果分析）
实验数据来自于何晓群<多元统计分析>第五版例题6.3 数据获取请关注公众号:321红绿灯 回复:例6.3 即可获得 区域公用事业的发展是地区綜合发展的重要组成部分,是促进社会发展的重要因素 ...
- 数学建模——多元统计分析例题及程序
多元统计分析例题及程序 主成分分析 简述 基本思想 计算步骤 例题 程序 程序理解 因子分析 相关性分析 回归分析 一元回归例题 聚类分析 主成分分析 简述 主成分分析(Principal Compo ...
- 【多元统计分析】Python实现对应分析
(一)题目要求 数据集包含地区生产总值的四个相关指标:x1劳动者报酬,x2生产税净额,x3固定资产折旧,x4营业盈余.对各个地区生产总值进行对应分析,揭示不同地区的生产总值构成特征.要求:画出对应分析 ...
- python 多元线性回归_多元统计分析之多元线性回归的R语言实现
多元统计分析之多元线性回归的R语言实现 多元统计分析--multivariate statistical analysis 研究客观事物中多个变量之间相互依赖的统计规律性.或从数学上说, 如果个体的观 ...
- 多元统计分析-主成分分析的原理与实现
目录 一.什么是主成分分析? 二.主成分分析的原理 三.主成分分析的应用 四.使用sklearn实现主成分分析 五.总结 一.什么是主成分分析? 主成分分析(Principal Component A ...
- 《多元统计分析与R语言》实验1【主成分分析】
<多元统计分析与R语言>实验1[主成分分析] 仅供参考 针对教材238页表格8-3(附件文档data_1_1.csv),进行主成分分析(代码和运行结果的截图). 链接:https://pa ...
最新文章
- FhqTreap的区间翻转
- linux unset命令,Linux unset命令
- 贪心 ---- C. Anu Has a Function位运算+贪心证明
- python读取数据库数据类型有哪些_数据库varchar 相当于python里的哪个类型
- arttemplate 不转义html,使用artTemplate模板引擎渲染错误
- Tensorflow【实战Google深度学习框架】TensorFlow模型的保存与恢复加载
- python发送文件到钉钉群_iOS python 自动化打包，并在钉钉群里发通知
- 施密特出售Google股票套现3.27亿美元
- 工作112:空字符串引发的问题
- matlab中结束脚本运行_MATLAB脱单指南
- mybatis开启字段自动映射为java驼峰命名规则
- （分层图）洛谷P4568[JLOI2011]飞行路线
- 用c语言输出英文字母表音标,26个英文字母表中文
- 数据分析实战（二）：流浪地球8W多条评论分析
- python根据经纬度画热力图_【python】使用python按照一定格式批量输出，地图热力图经纬度格式生成器...
- Manjaro安装与基本配置
- 交通管理PDA手持机
- TouchGFX 中文字体添加(直入正题）
- CSharp(C#)语言_第六章（表达式和运算符）
- 淘宝经典差评和爆笑回复