主成分分析与因子分析法
问题一
# 导包
import numpy as np
import pandas as pd
# 传入矩阵参数
R1 = np.array([[1, 0.577, 0.509, 0.387, 0.462], [0.577, 1, 0.599, 0.389, 0.322], [0.509, 0.599, 1, 0.436, 0.426], [0.387, 0.389, 0.436, 1, 0.523], [0.462, 0.322, 0.426, 0.523, 1]])
# 求矩阵的参数与特征向量
eigValues1, eigVects1 = np.linalg.eig(R1)
eigValues1, eigVects1
(array([2.85671099, 0.80916372, 0.34294995, 0.53967524, 0.4515001 ]),array([[ 0.46360519, 0.24033901, 0.45126217, 0.61170545, -0.38663456],[ 0.45710783, 0.50930473, -0.67622331, -0.17818949, -0.20647436],[ 0.47017564, 0.26044828, 0.40000721, -0.33505645, 0.66244469],[ 0.42145876, -0.5256649 , 0.17559859, -0.54076277, -0.47200602],[ 0.42122445, -0.58196989, -0.38502448, 0.43517554, 0.38243876]]))
设第iii主成分Yi=∑k=1nekXkY_i=\sum_{k=1}^{n} e_kX_kYi=∑k=1nekXk，根据数学推导，可令特征值由大到小所对应的正交单位化特征向量即为对应主成分的组合系数，得到的主成分满足Y1Y_1Y1尽可能包含所有信息，Yi(2≤i≤n)Y_i(2 \leq i \leq n)Yi(2≤i≤n)尽可能包含除了Yk(k≤i−1)Y_k(k \leq i-1)Yk(k≤i−1)之外的信息，而各主成分的方差等于对应的特征值。
根据相关矩阵做主成分分析，相当于对标准化的数据的协方差矩阵做主成分分析，得到特征值与特征向量为：
| 特征值 | 特征向量 |
|---|---|
| λ1\lambda_1λ1=2.857 | e1e_1e1=(0.467,0.457,0.470,0.421,0.421) |
| λ2\lambda_2λ2=0.809 | e2e_2e2=(0.240,0.509,0.260,-0.526,-0.581) |
| λ3\lambda_3λ3=0.540 | e3e_3e3=(0.612,-0.178,-0.335,-0.541,0.435) |
| λ4\lambda_4λ4=0.452 | e4e_4e4=(-0.387,-0.206,0.662,-0.472,0.382) |
| λ5\lambda_5λ5=0.343 | e5e_5e5=(0.451,-0.676,0.400,0.176,-0.385) |
因此X∗X^{*}X∗的五个主成分为：
Yk=ekX∗Y_k=e_kX^{*} Yk=ekX∗
# 计算贡献率
def contributionRate(values):v = list(values)l = len(v)result = np.array([v[i]*100 / sum(v) for i in range(l)])add = np.empty([l, ])add[0] = result[0]for i in range(1, l):add[i] = add[i-1] + result[i]return [result, add]
contribution = pd.DataFrame() # 利用变量名和特征值建立一个数据框
contribution['eigValues'] = eigValues1 # 特征值
contribution['贡献率'] = contributionRate(eigValues1)[0]
contribution['累计贡献率'] = contributionRate(eigValues1)[1]
contribution
| eigValues | 贡献率 | 累计贡献率 | |
|---|---|---|---|
| 0 | 2.856711 | 57.134220 | 57.134220 |
| 1 | 0.809164 | 16.183274 | 73.317494 |
| 2 | 0.342950 | 6.858999 | 80.176493 |
| 3 | 0.539675 | 10.793505 | 90.969998 |
| 4 | 0.451500 | 9.030002 | 100.000000 |
由于第kkk主成分的贡献率=λk∑i=15λi\frac{\lambda_k}{\sum_{i=1}^{5} \lambda_i}∑i=15λiλk，计算得到各主成分贡献率依次是57.1%，16.2%，6.86%，10.80%，9.03%。前两个主成分的累计贡献率为73.3%。
第一主成分意义是所有指标的加权和，当第一主成分较大时，可以推断各项反弹率均比较大，故可称为“总反弹率”因子。第二主成分 较大时， 反弹率比较小， 较大，因此第二主成分可称为“石油股票反弹率”因子。
问题二
字段信息：
| 指标 | 含义 |
|---|---|
| X1X_1X1 | 人均粮食支出 |
| X2X_2X2 | 人均副食支出 |
| X3X_3X3 | 人均烟酒茶支出 |
| X4X_4X4 | 人均其他副食支出 |
| X5X_5X5 | 人均衣着商品支出 |
| X6X_6X6 | 人均日用品支出 |
| X7X_7X7 | 人均燃料支出 |
| X8X_8X8 | 人均非商品支出 |
第一问：求相关系数矩阵R
# 数据预处理
data = pd.read_excel('4.2.xlsx',index_col=0)
# 0均值规范化
data = (data - data.mean()) / data.std()
#相关系数矩阵
R2 = data.corr()
R2
| X1 | X2 | X3 | X4 | X5 | X6 | X7 | X8 | |
|---|---|---|---|---|---|---|---|---|
| X1 | 1.000000 | 0.333647 | -0.054539 | -0.061254 | -0.289361 | 0.198796 | 0.348698 | 0.318677 |
| X2 | 0.333647 | 1.000000 | -0.022902 | 0.398931 | -0.156304 | 0.711134 | 0.413595 | 0.834952 |
| X3 | -0.054539 | -0.022902 | 1.000000 | 0.533329 | 0.496763 | 0.032830 | -0.139086 | -0.258358 |
| X4 | -0.061254 | 0.398931 | 0.533329 | 1.000000 | 0.698424 | 0.467917 | -0.171274 | 0.312757 |
| X5 | -0.289361 | -0.156304 | 0.496763 | 0.698424 | 1.000000 | 0.280129 | -0.208277 | -0.081234 |
| X6 | 0.198796 | 0.711134 | 0.032830 | 0.467917 | 0.280129 | 1.000000 | 0.416821 | 0.701586 |
| X7 | 0.348698 | 0.413595 | -0.139086 | -0.171274 | -0.208277 | 0.416821 | 1.000000 | 0.398868 |
| X8 | 0.318677 | 0.834952 | -0.258358 | 0.312757 | -0.081234 | 0.701586 | 0.398868 | 1.000000 |
从上表可以看出，变量间有较强的关系，因此适合做主成分分析和因子分析，以此达到降维的目的。
第二问：主成分分析
# 计算特征值和特征向量
eigValues2, eigVects2 = np.linalg.eig(R2)
eig = pd.DataFrame() # 利用变量名和特征值建立一个数据框
eig['names'] = data.columns # 列名
eig['eigValues'] = eigValues2 # 特征值
eig
| names | eigValues | |
|---|---|---|
| 0 | X1 | 3.096288 |
| 1 | X2 | 2.367223 |
| 2 | X3 | 0.919987 |
| 3 | X4 | 0.705925 |
| 4 | X5 | 0.498442 |
| 5 | X6 | 0.051470 |
| 6 | X7 | 0.130776 |
| 7 | X8 | 0.229888 |
# 确定公共因子个数
from math import sqrt
for k in range(1, 9): if eig['eigValues'][:k].sum() / eig['eigValues'].sum() >= 0.8: # 如果解释度达到80%, 结束循环print(k)break
4
由于第kkk主成分的贡献率=λk∑i=18λi\frac{\lambda_k}{\sum_{i=1}^{8} \lambda_i}∑i=18λiλk
# 计算贡献率
contribution = pd.DataFrame() # 利用变量名和特征值建立一个数据框
contribution['names'] = data.columns # 列名
contribution['eigValues'] = eigValues2 # 特征值
contribution['贡献率'] = contributionRate(eigValues2)[0]
contribution['累计贡献率'] = contributionRate(eigValues2)[1]
contribution
| names | eigValues | 贡献率 | 累计贡献率 | |
|---|---|---|---|---|
| 0 | X1 | 3.096288 | 38.703604 | 38.703604 |
| 1 | X2 | 2.367223 | 29.590288 | 68.293892 |
| 2 | X3 | 0.919987 | 11.499842 | 79.793734 |
| 3 | X4 | 0.705925 | 8.824067 | 88.617801 |
| 4 | X5 | 0.498442 | 6.230529 | 94.848330 |
| 5 | X6 | 0.051470 | 0.643369 | 95.491699 |
| 6 | X7 | 0.130776 | 1.634697 | 97.126396 |
| 7 | X8 | 0.229888 | 2.873604 | 100.000000 |
第三问：因子分析
# 构造初始因子载荷矩阵
col0 = list(sqrt(eigValues2[0]) * eigVects2[:,0]) #因子载荷矩阵第1列
col1 = list(sqrt(eigValues2[1]) * eigVects2[:,1]) #因子载荷矩阵第2列
col2 = list(sqrt(eigValues2[2]) * eigVects2[:,2]) #因子载荷矩阵第3列
col3 = list(sqrt(eigValues2[3]) * eigVects2[:,3]) #因子载荷矩阵第4列
A = pd.DataFrame([col0, col1, col2, col3]).T #构造因子载荷矩阵A
A.columns=['X1','X2','X3','X4'] #因子载荷矩阵A的公共因子
A
| X1 | X2 | X3 | X4 | |
|---|---|---|---|---|
| 0 | 0.439215 | 0.371164 | -0.665578 | 0.316560 |
| 1 | 0.913659 | 0.057862 | 0.068351 | 0.188935 |
| 2 | -0.032518 | -0.731499 | -0.554221 | -0.027204 |
| 3 | 0.447107 | -0.827879 | 0.020887 | 0.194140 |
| 4 | 0.038175 | -0.885373 | 0.046123 | -0.239764 |
| 5 | 0.866903 | -0.207209 | 0.139412 | -0.188390 |
| 6 | 0.558060 | 0.401080 | -0.274694 | -0.645366 |
| 7 | 0.896234 | 0.133982 | 0.260200 | 0.148706 |
# 建立因子模型
h = np.zeros(8) #变量共同度，反映变量对共同因子的依赖程度，越接近1，说明公共因子解释程度越高，因子分析效果越好
D = np.mat(np.eye(8))#特殊因子方差，因子的方差贡献度 ，反映公共因子对变量的贡献，衡量公共因子的相对重要性
A = np.mat(A) #将因子载荷阵A矩阵化
for i in range(8):a = A[i,:] * A[i,:].T #A的元的行平方和h[i] = a[0,0] #计算变量X共同度,描述全部公共因子F对变量X_i的总方差所做的贡献，及变量X_i方差中能够被全体因子解释的部分D[i,i] = 1 - a[0,0] #因为自变量矩阵已经标准化后的方差为1，即Var(X_i)=第i个共同度h_i + 第i个特殊因子方差
将因子表示成变量的线性组合
def varimax(Phi, gamma = 1.0, q =10, tol = 1e-6): #定义方差最大旋转函数p, k = Phi.shape #给出矩阵Phi的总行数，总列数R = np.eye(k) #给定一个k*k的单位矩阵d = 0for i in range(q):d_old = dLambda = np.dot(Phi, R)#矩阵乘法u, s, vh = np.linalg.svd(np.dot(Phi.T, np.asarray(Lambda) ** 3 - (gamma / p) * np.dot(Lambda, np.diag(np.diag(np.dot(Lambda.T,Lambda)))))) #奇异值分解svdR = np.dot(u,vh)#构造正交矩阵Rd = np.sum(s)#奇异值求和if d_old != 0 and d / d_old:return np.dot(Phi, R)#返回旋转矩阵Phi*Rrotation_mat = varimax(A) # 调用方差最大旋转函数
rotation_mat = pd.DataFrame(rotation_mat) # 数据框化
rotation_mat
| 0 | 1 | 2 | 3 | |
|---|---|---|---|---|
| 0 | 0.203150 | 0.062124 | -0.892652 | -0.178661 |
| 1 | 0.889045 | -0.029496 | -0.262154 | -0.135988 |
| 2 | -0.214454 | -0.871857 | -0.194565 | 0.008240 |
| 3 | 0.481798 | -0.792082 | 0.077324 | 0.240652 |
| 4 | 0.033408 | -0.819782 | 0.413444 | 0.029366 |
| 5 | 0.804396 | -0.274344 | 0.064507 | -0.350578 |
| 6 | 0.255416 | 0.127864 | -0.165170 | -0.924869 |
| 7 | 0.933336 | 0.109739 | -0.110305 | -0.125212 |
# 计算因子得分
data_ = np.mat(data) #矩阵化处理
factor_score = (data_).dot(A) #计算因子得分
factor_score = pd.DataFrame(factor_score)#数据框化
factor_score.columns = ['X1','X2','X3','X4'] #对因子变量进行命名
factor_score.index = data.index
factor_score['sum'] = factor_score.sum(axis=1)
结果排序
factor_score = factor_score.sort_values(by = ['sum'], ascending=False)
factor_score
| X1 | X2 | X3 | X4 | sum | |
|---|---|---|---|---|---|
| 地区 | |||||
| 广东 | 12.134223 | 3.505384 | -0.780327 | -1.320523 | 13.538757 |
| 海南 | 2.465522 | 4.890210 | -1.567251 | 1.625074 | 7.413554 |
| 广西 | 1.859072 | 1.900515 | 1.539035 | -0.063227 | 5.235395 |
| 上海 | 5.716003 | -3.940215 | 1.603553 | 1.666736 | 5.046078 |
| 福建 | 2.030434 | 2.114319 | -0.799141 | 0.606531 | 3.952143 |
| 湖南 | -0.911514 | 0.849711 | 1.107383 | -0.187768 | 0.857812 |
| 江西 | -3.450478 | 3.243735 | 0.989421 | -0.075377 | 0.707301 |
| 四川 | -0.237351 | -0.528700 | 1.121904 | 0.133922 | 0.489775 |
| 天津 | 0.766193 | -0.729420 | -0.647183 | 1.044525 | 0.434115 |
| 北京 | 3.153503 | -4.443551 | 1.185859 | 0.457701 | 0.353514 |
| 云南 | -1.173142 | 0.454616 | -0.215400 | 0.962859 | 0.028933 |
| 新疆 | -1.440261 | 0.648165 | 1.237420 | -0.432477 | 0.012847 |
| 江苏 | 0.269736 | -0.174160 | 0.491318 | -0.609942 | -0.023048 |
| 湖北 | -0.299919 | 0.907169 | -0.123593 | -0.584738 | -0.101081 |
| 陕西 | -1.078188 | 0.434750 | 0.037484 | 0.254255 | -0.351699 |
| 浙江 | 2.665959 | -2.113285 | -0.145554 | -0.905556 | -0.498436 |
| 贵州 | -2.215974 | 0.658934 | 0.294098 | 0.350311 | -0.912631 |
| 安徽 | -1.961910 | 0.677642 | -0.210394 | 0.521765 | -0.972897 |
| 河北 | -0.690123 | 0.454840 | -0.414944 | -0.635875 | -1.286102 |
| 吉林 | -2.276431 | 1.323683 | 0.105062 | -0.566823 | -1.414510 |
| 辽宁 | 0.079538 | -1.508676 | 0.178047 | -0.325763 | -1.576853 |
| 山东 | -0.248310 | -1.049270 | -0.064522 | -0.451012 | -1.813114 |
| 宁夏 | -0.757336 | -0.993462 | 0.136768 | -0.373849 | -1.987879 |
| 甘肃 | -2.080281 | 0.296880 | -0.726723 | 0.027185 | -2.482939 |
| 内蒙古 | -2.213952 | 0.203457 | -0.289416 | -0.190353 | -2.490265 |
| 河南 | -2.614718 | 0.540896 | -0.702427 | 0.131479 | -2.644770 |
| 山西 | -2.964064 | -0.257612 | 0.155919 | 0.194874 | -2.870882 |
| 黑龙江 | -2.332689 | 0.157537 | 0.222618 | -1.009044 | -2.961577 |
| 青海 | -1.959088 | 0.027927 | -1.477425 | 0.137006 | -3.271581 |
| 西藏 | -0.234453 | -7.552021 | -2.241588 | -0.381896 | -10.409958 |
第一主成分是所有指标除了人均副食支出的加权和，当第一主成分较大时，可以推断各项生活日用支出均比较大，故可称为“生活日用支出”因子。第二主成分 较大时人均粮食、副食、燃料、非商品支出比较小，人均烟酒茶、其他副食、衣着商品、日用品支出较大，因此第二主成分可称为“非必要支出”因子。
因子分析是研究如何以最少的信息丢失将众多原有变量浓缩成少数几个因子,如何使因子具有一定的命名解释性的多元统计分析方法.其核心是用较少的相互独立的因子反映原有变量的绝大部分信息。
其矩阵形式为X=AF+εX=AF+\varepsilonX=AF+ε
对数据进行因子分析前要先进行KMO球形检验，得到P值=0.000，非常显著，故满足因子分析的前提条件。
问题三
下表3是25个家庭的成年长子的头长(X1)、头宽(X2)与成年次子的头长(Y1)和头宽(Y2)的观测数据.试分别从样本协方差矩阵Σ和样本相关系数矩阵R出发做典型相关分析,求各典型变量对及典型相关系数,检验各典型变量对是否显著相关( ＝0.05).两种情况下的结果有何异同?
# 导包
df = pd.read_excel('4.3.xlsx')
# 输入
X = df[['X1', 'X2']]
y = df[['Y1', 'Y2']]
# 建立模型
import scipy
from sklearn.cross_decomposition import CCA
def cca(X, Y, k): l = len(k)reu = []for i in range(l):cca = CCA(n_components=k[i])# 训练数据cca.fit(X, y)# 降维操作X_train_r, y_train_r = cca.transform(X, y)reu.append(scipy.stats.pearsonr(X_train_r[:, k[i]-1], y_train_r[:, k[i]-1])) # 输出相关系数return reu
analysis = pd.DataFrame()
r = cca(X, y, [1, 2])
analysis['corr'] = r[0]
analysis['p-value'] = r[1]
analysis
| corr | p-value | |
|---|---|---|
| 0 | 0.788508 | 0.053740 |
| 1 | 0.000003 | 0.798626 |
无论是标准化的数据还是未标准化的数据，第一对典型相关变量的典型相关系数均为0.789，P值为0.056可以显著认为第一对典型变量能够提供X和Y的相关信息；第二对典型相关变量的典型相关系数为0.0003，P值为0.797，可以显著认为第二对典型变量不能够提供X和Y的相关信息。
主成分分析与因子分析法相关推荐
- 主成分分析、因子分析、聚类分析的比较与应用
听说这是一篇论文 不过我没详细看. 一.概述 主成分分析就是将多项指标转化为少数几项综合指标,用综合指标来解释多变量的方差- 协方差结构.综合指标即为主成分.所得出的少数几个主成分,要尽可能多地保留原 ...
- 主成分分析、因子分析和聚类分析的区别
主成分分析就是设法将原来众多具有一定相关性(比如P个指标),重新组合成一组新的互相无关的综合指标来代替原来的指标. 综合指标即为主成分.所得出的少数几个主成分,要尽可能多地保留原始变量的信息,且彼此不 ...
- 因子分析 factor analysis (七) ：因子分析法与主成分分析的异同
因子分析系列博文: 因子分析 factor analysis (一 ):模型的理论推导 因子分析 factor analysis (二 ) : 因子分析模型 因子分析 factor analysis ...
- 机器学习-特征抽取(主成分分析法/因子分析法/非负矩阵因子分解NMF算法)
1.特征抽取: 特征抽取是机器学习中另一种十分有用的方法,它与特性选择不同,特征抽取是对数据的特征进行概括和总结,而特性选择则主要是对数据中的不同特征进行比较和选取. 特征抽取是机器学习技术中的一个常 ...
- 主成分分析和因子分析及其在R中的…
1 主成分分析和因子分析比较 主成分分析和探索性因子分析是两种用来探索和简化多变量复杂关系的常用方法,它们之间有联系也有区别. 主成分分析(PCA)是一种数据降维方法,它能将大量相关变量转化为一组很 ...
- 主成分分析和因子分析十大不同点
主成分分析和因子分析无论从算法上还是应用上都有着比较相似之处,本文结合以往资料以及自己的理解总结了以下十大不同之处,适合初学者学习之用. 1.原理不同 主成分分析基本原理:利用降维(线性变换)的思想, ...
- 关于聚类分析、判别分析、主成分分析、因子分析等多元统计分析方法
转载自:http://blog.csdn.net/nieson2012/article/details/25408421 主成分分析与因子分析的区别 1. 目的不同: 因子分析把诸多变量看成由对每一个 ...
- 主成分分析与因子分析之比较及实证分析
一.问题的提出 在科学研究或日常生活中,常常需要判断某一事物在同类事物中的好坏.优劣程度及其发展规律等问题.而影响事物的特征及其发展规律的因素(指标)是多方面的,因此,在对该事物进行研究时,为了能更全 ...
- R语言学习笔记——高级篇：第十四章-主成分分析和因子分析
R语言 R语言学习笔记--高级篇:第十四章-主成分分析和因子分析 文章目录 R语言 前言 一.R中的主成分和因子分析 二.主成分分析 2.1.判断主成分的个数 2.2.提取主成分 2.3.主成分旋转 ...
- SPSS（十一）SPSS信息浓缩技术--主成分分析、因子分析（图文+数据集）
SPSS(十一)信息浓缩技术--主成分分析.因子分析(图文+数据集) 当我们的自变量存在多重共线性,表现为进行回归时候方程系数估计不正常以及方程检验结果不正常,也许我们可以使用变量挑选的办法(手动挑选 ...
最新文章
- DPM2012恢复单个Exchange2010用户邮箱
- 2017.8.17 开始了我的QT 学习。
- Flowable 数据库表结构 ACT_HI_PROCINST
- deepfakes怎么用_[mcj]deepfakesApp使用说明（2）
- one stage 与two stage解释
- Windows10/Servers2016应用商店恢复/安装
- java clone 深拷贝_Java clone() 浅拷贝 深拷贝
- Ext 2.0布局实例
- SVN回滚到指定旧版本操作指南
- 马哥linux作业,马哥linux0803作业内容
- 【2018滴滴】寻找丑数
- 数据分析师面试题目_数据分析师面试题目
- python中的wx模块
- 独木舟上的旅行java_独木舟上的旅行
- 网络安全认证与加密协议算法整合
- 备份恢复Lesson 08. Using RMAN-Encrypted Backups
- 美团运维SRE+运维开发一面面经汇总
- 思维导图的分类 利用思维导图绘制学习知识方法介绍
- 王兴是怎么看待共享单车这块业务的
- 软件开发过程大观——软件开发过程改进为什么能帮助软件质量提升？