因子分析用Python做的一個典型例子
一、實驗目的
采用合适的資料分析方法對下面的題進行解答
![](https://img.laitimes.com/img/9ZDMuAjOiMmIsIjOiQnIsISPrdEZwZ1Rh5WNXp1bwNjW1ZUba9VZwlHdsATOfd3bkFGazxCMx8VesATMfhHLlN3XnxCMwEzX0xiRGZkRGZ0Xy9GbvNGLpZTY1EmMZVDUSFTU4VFRR9Fd4VGdsYTMfVmepNHLrJXYtJXZ0F2dvwVZnFWbp1zczV2YvJHctM3cv1Ce-cmbw5COhdzYhRWMjFzY5MzYjNDMhVWN5ETN5kjZzYWMxQ2Mk9CXwMzLchDMxIDMy8CXn9Gbi9CXzV2Zh1WavwVbvNmLvR3YxUjL1M3Lc9CX6MHc0RHaiojIsJye.png)
二、實驗要求
采用因子分析方法,根據48位應聘者的15項名額得分,選出6名最優秀的應聘者。
三、代碼
import pandas as pd
import numpy as np
import math as math
import numpy as np
from numpy import *
from scipy.stats import bartlett
from factor_analyzer import *
import numpy.linalg as nlg
from sklearn.cluster import KMeans
from matplotlib import cm
import matplotlib.pyplot as plt
def main():
df=pd.read_csv("./data/applicant.csv")
# print(df)
df2=df.copy()
print("\n原始資料:\n",df2)
del df2['ID']
# print(df2)
# 皮爾森相關系數
df2_corr=df2.corr()
print("\n相關系數:\n",df2_corr)
#熱力圖
cmap = cm.Blues
# cmap = cm.hot_r
fig=plt.figure()
ax=fig.add_subplot(111)
map = ax.imshow(df2_corr, interpolation='nearest', cmap=cmap, vmin=0, vmax=1)
plt.title('correlation coefficient--headmap')
ax.set_yticks(range(len(df2_corr.columns)))
ax.set_yticklabels(df2_corr.columns)
ax.set_xticks(range(len(df2_corr)))
ax.set_xticklabels(df2_corr.columns)
plt.colorbar(map)
plt.show()
# KMO測度
def kmo(dataset_corr):
corr_inv = np.linalg.inv(dataset_corr)
nrow_inv_corr, ncol_inv_corr = dataset_corr.shape
A = np.ones((nrow_inv_corr, ncol_inv_corr))
for i in range(0, nrow_inv_corr, 1):
for j in range(i, ncol_inv_corr, 1):
A[i, j] = -(corr_inv[i, j]) / (math.sqrt(corr_inv[i, i] * corr_inv[j, j]))
A[j, i] = A[i, j]
dataset_corr = np.asarray(dataset_corr)
kmo_num = np.sum(np.square(dataset_corr)) - np.sum(np.square(np.diagonal(A)))
kmo_denom = kmo_num + np.sum(np.square(A)) - np.sum(np.square(np.diagonal(A)))
kmo_value = kmo_num / kmo_denom
return kmo_value
print("\nKMO測度:", kmo(df2_corr))
# 巴特利特球形檢驗
df2_corr1 = df2_corr.values
print("\n巴特利特球形檢驗:", bartlett(df2_corr1[0], df2_corr1[1], df2_corr1[2], df2_corr1[3], df2_corr1[4],
df2_corr1[5], df2_corr1[6], df2_corr1[7], df2_corr1[8], df2_corr1[9],
df2_corr1[10], df2_corr1[11], df2_corr1[12], df2_corr1[13], df2_corr1[14]))
# 求特征值和特征向量
eig_value, eigvector = nlg.eig(df2_corr) # 求矩陣R的全部特征值,構成向量
eig = pd.DataFrame()
eig['names'] = df2_corr.columns
eig['eig_value'] = eig_value
eig.sort_values('eig_value', ascending=False, inplace=True)
print("\n特征值\n:",eig)
eig1=pd.DataFrame(eigvector)
eig1.columns = df2_corr.columns
eig1.index = df2_corr.columns
print("\n特征向量\n",eig1)
# 求公因子個數m,使用前m個特征值的比重大于85%的标準,選出了公共因子是五個
for m in range(1, 15):
if eig['eig_value'][:m].sum() / eig['eig_value'].sum() >= 0.85:
print("\n公因子個數:", m)
break
# 因子載荷陣
A = np.mat(np.zeros((15, 5)))
i = 0
j = 0
while i < 5:
j = 0
while j < 15:
A[j:, i] = sqrt(eig_value[i]) * eigvector[j, i]
j = j + 1
i = i + 1
a = pd.DataFrame(A)
a.columns = ['factor1', 'factor2', 'factor3', 'factor4', 'factor5']
a.index = df2_corr.columns
print("\n因子載荷陣\n", a)
fa = FactorAnalyzer(n_factors=5)
fa.loadings_ = a
# print(fa.loadings_)
print("\n特殊因子方差:\n", fa.get_communalities()) # 特殊因子方差,因子的方差貢獻度 ,反映公共因子對變量的貢獻
var = fa.get_factor_variance() # 給出貢獻率
print("\n解釋的總方差(即貢獻率):\n", var)
# 因子旋轉
rotator = Rotator()
b = pd.DataFrame(rotator.fit_transform(fa.loadings_))
b.columns = ['factor1', 'factor2', 'factor3', 'factor4', 'factor5']
b.index = df2_corr.columns
print("\n因子旋轉:\n", b)
# 因子得分
X1 = np.mat(df2_corr)
X1 = nlg.inv(X1)
b = np.mat(b)
factor_score = np.dot(X1, b)
factor_score = pd.DataFrame(factor_score)
factor_score.columns = ['factor1', 'factor2', 'factor3', 'factor4', 'factor5']
factor_score.index = df2_corr.columns
print("\n因子得分:\n", factor_score)
fa_t_score = np.dot(np.mat(df2), np.mat(factor_score))
print("\n應試者的五個因子得分:\n",pd.DataFrame(fa_t_score))
# 綜合得分
wei = [[0.50092], [0.137087], [0.097055], [0.079860], [0.049277]]
fa_t_score = np.dot(fa_t_score, wei) / 0.864198
fa_t_score = pd.DataFrame(fa_t_score)
fa_t_score.columns = ['綜合得分']
fa_t_score.insert(0, 'ID', range(1, 49))
print("\n綜合得分:\n", fa_t_score)
print("\n綜合得分:\n", fa_t_score.sort_values(by='綜合得分', ascending=False).head(6))
plt.figure()
ax1=plt.subplot(111)
X=fa_t_score['ID']
Y=fa_t_score['綜合得分']
plt.bar(X,Y,color="#87CEFA")
# plt.bar(X, Y, color="red")
plt.title('result00')
ax1.set_xticks(range(len(fa_t_score)))
ax1.set_xticklabels(fa_t_score.index)
plt.show()
fa_t_score1=pd.DataFrame()
fa_t_score1=fa_t_score.sort_values(by='綜合得分',ascending=False).head()
ax2 = plt.subplot(111)
X1 = fa_t_score1['ID']
Y1 = fa_t_score1['綜合得分']
plt.bar(X1, Y1, color="#87CEFA")
# plt.bar(X1, Y1, color='red')
plt.title('result01')
plt.show()
if __name__ == '__main__':
main()
四、實驗步驟
(1)引入資料,資料标準化
因為資料是面試中的得分,量綱相同,并且資料的分布無異常值,是以資料可以不進行标準化。
(2)建立相關系數矩陣
計算皮爾森相關系數,從熱圖中可以明顯看出變量間存在的相關性。
進行相關系數矩陣檢驗——KMO測度和巴特利特球體檢驗:
KMO值:0.9以上非常好;0.8以上好;0.7一般;0.6差;0.5很差;0.5以下不能接受;巴特利球形檢驗的值範圍在0-1,越接近1,使用因子分析效果越好。
通過觀察上面的計算結果,可以知道,KMO值為0.783775605643526,在較好的範圍内,并且巴特利球形檢驗的值接近1,所有可以使用因子分析。
(3)求解特征值及相應特征向量
求公因子個數m,使用前m個特征值的比重大于85%的标準,選出了公共因子是五個。
(4)因子載荷陣
由上可以看出,選擇5個公共因子,從方差貢獻率可以看出,其中第一個公因子解釋了總體方差的50.092%,四個公共因子的方差貢獻率為86.42%,可以較好的解釋總體方差。
(5)因子旋轉
(6)因子得分
(7)根據應聘者的五個因子得分,按照貢獻率進行權重,得到最終各應試者的綜合得分,然後選出前六個得分最高的應聘者。
是以我們用因子分析産生的前六名分别是:40,39,22,2,10,23
![](https://img.laitimes.com/img/9ZDMuAjOiMmIsIjOiQnIsISPrdEZwZ1Rh5WNXp1bwNjW1ZUba9VZwlHdsATOfd3bkFGazxCMx8VesATMfhHLlN3XnxCMwEzX0xiRGZkRGZ0Xy9GbvNGLpZTY1EmMZVDUSFTU4VFRR9Fd4VGdsYTMfVmepNHLrJXYtJXZ0F2dvwVZnFWbp1zczV2YvJHctM3cv1Ce-cmbw5COhdzYhRWMjFzY5MzYjNDMhVWN5ETN5kjZzYWMxQ2Mk9CXwMzLchDMxIDMy8CXn9Gbi9CXzV2Zh1WavwVbvNmLvR3YxUjL1M3Lc9CX6MHc0RHaiojIsJye.png)
import pandas as pd
import numpy as np
import math as math
import numpy as np
from numpy import *
from scipy.stats import bartlett
from factor_analyzer import *
import numpy.linalg as nlg
from sklearn.cluster import KMeans
from matplotlib import cm
import matplotlib.pyplot as plt
def main():
df=pd.read_csv("./data/applicant.csv")
# print(df)
df2=df.copy()
print("\n原始資料:\n",df2)
del df2['ID']
# print(df2)
# 皮爾森相關系數
df2_corr=df2.corr()
print("\n相關系數:\n",df2_corr)
#熱力圖
cmap = cm.Blues
# cmap = cm.hot_r
fig=plt.figure()
ax=fig.add_subplot(111)
map = ax.imshow(df2_corr, interpolation='nearest', cmap=cmap, vmin=0, vmax=1)
plt.title('correlation coefficient--headmap')
ax.set_yticks(range(len(df2_corr.columns)))
ax.set_yticklabels(df2_corr.columns)
ax.set_xticks(range(len(df2_corr)))
ax.set_xticklabels(df2_corr.columns)
plt.colorbar(map)
plt.show()
# KMO測度
def kmo(dataset_corr):
corr_inv = np.linalg.inv(dataset_corr)
nrow_inv_corr, ncol_inv_corr = dataset_corr.shape
A = np.ones((nrow_inv_corr, ncol_inv_corr))
for i in range(0, nrow_inv_corr, 1):
for j in range(i, ncol_inv_corr, 1):
A[i, j] = -(corr_inv[i, j]) / (math.sqrt(corr_inv[i, i] * corr_inv[j, j]))
A[j, i] = A[i, j]
dataset_corr = np.asarray(dataset_corr)
kmo_num = np.sum(np.square(dataset_corr)) - np.sum(np.square(np.diagonal(A)))
kmo_denom = kmo_num + np.sum(np.square(A)) - np.sum(np.square(np.diagonal(A)))
kmo_value = kmo_num / kmo_denom
return kmo_value
print("\nKMO測度:", kmo(df2_corr))
# 巴特利特球形檢驗
df2_corr1 = df2_corr.values
print("\n巴特利特球形檢驗:", bartlett(df2_corr1[0], df2_corr1[1], df2_corr1[2], df2_corr1[3], df2_corr1[4],
df2_corr1[5], df2_corr1[6], df2_corr1[7], df2_corr1[8], df2_corr1[9],
df2_corr1[10], df2_corr1[11], df2_corr1[12], df2_corr1[13], df2_corr1[14]))
# 求特征值和特征向量
eig_value, eigvector = nlg.eig(df2_corr) # 求矩陣R的全部特征值,構成向量
eig = pd.DataFrame()
eig['names'] = df2_corr.columns
eig['eig_value'] = eig_value
eig.sort_values('eig_value', ascending=False, inplace=True)
print("\n特征值\n:",eig)
eig1=pd.DataFrame(eigvector)
eig1.columns = df2_corr.columns
eig1.index = df2_corr.columns
print("\n特征向量\n",eig1)
# 求公因子個數m,使用前m個特征值的比重大于85%的标準,選出了公共因子是五個
for m in range(1, 15):
if eig['eig_value'][:m].sum() / eig['eig_value'].sum() >= 0.85:
print("\n公因子個數:", m)
break
# 因子載荷陣
A = np.mat(np.zeros((15, 5)))
i = 0
j = 0
while i < 5:
j = 0
while j < 15:
A[j:, i] = sqrt(eig_value[i]) * eigvector[j, i]
j = j + 1
i = i + 1
a = pd.DataFrame(A)
a.columns = ['factor1', 'factor2', 'factor3', 'factor4', 'factor5']
a.index = df2_corr.columns
print("\n因子載荷陣\n", a)
fa = FactorAnalyzer(n_factors=5)
fa.loadings_ = a
# print(fa.loadings_)
print("\n特殊因子方差:\n", fa.get_communalities()) # 特殊因子方差,因子的方差貢獻度 ,反映公共因子對變量的貢獻
var = fa.get_factor_variance() # 給出貢獻率
print("\n解釋的總方差(即貢獻率):\n", var)
# 因子旋轉
rotator = Rotator()
b = pd.DataFrame(rotator.fit_transform(fa.loadings_))
b.columns = ['factor1', 'factor2', 'factor3', 'factor4', 'factor5']
b.index = df2_corr.columns
print("\n因子旋轉:\n", b)
# 因子得分
X1 = np.mat(df2_corr)
X1 = nlg.inv(X1)
b = np.mat(b)
factor_score = np.dot(X1, b)
factor_score = pd.DataFrame(factor_score)
factor_score.columns = ['factor1', 'factor2', 'factor3', 'factor4', 'factor5']
factor_score.index = df2_corr.columns
print("\n因子得分:\n", factor_score)
fa_t_score = np.dot(np.mat(df2), np.mat(factor_score))
print("\n應試者的五個因子得分:\n",pd.DataFrame(fa_t_score))
# 綜合得分
wei = [[0.50092], [0.137087], [0.097055], [0.079860], [0.049277]]
fa_t_score = np.dot(fa_t_score, wei) / 0.864198
fa_t_score = pd.DataFrame(fa_t_score)
fa_t_score.columns = ['綜合得分']
fa_t_score.insert(0, 'ID', range(1, 49))
print("\n綜合得分:\n", fa_t_score)
print("\n綜合得分:\n", fa_t_score.sort_values(by='綜合得分', ascending=False).head(6))
plt.figure()
ax1=plt.subplot(111)
X=fa_t_score['ID']
Y=fa_t_score['綜合得分']
plt.bar(X,Y,color="#87CEFA")
# plt.bar(X, Y, color="red")
plt.title('result00')
ax1.set_xticks(range(len(fa_t_score)))
ax1.set_xticklabels(fa_t_score.index)
plt.show()
fa_t_score1=pd.DataFrame()
fa_t_score1=fa_t_score.sort_values(by='綜合得分',ascending=False).head()
ax2 = plt.subplot(111)
X1 = fa_t_score1['ID']
Y1 = fa_t_score1['綜合得分']
plt.bar(X1, Y1, color="#87CEFA")
# plt.bar(X1, Y1, color='red')
plt.title('result01')
plt.show()
if __name__ == '__main__':
main()