案例：
将1660个人组成的样本按照心理健康状况[good->slight->moderate->damaged]和社会经济状况[Level_A(高)-->LevelE(低)]进行交叉分类如下表所示
原案例传送门
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
# 频数矩阵
X = np.array([[121,57 ,72 ,36 ,21 ],
[188,105,141,97 ,71 ],
[112,65 ,77 ,54 ,54 ],
[86 ,60 ,94 ,78 ,71 ]
])
df = pd.DataFrame(X,
columns=['Level_A','Level_B','Level_C','Level_D','Level_E'],
index=['good','slight','moderate','damaged'])
print(df)
Level_A Level_B Level_C Level_D Level_E
good 121 57 72 36 21
slight 188 105 141 97 71
moderate 112 65 77 54 54
damaged 86 60 94 78 71
# 频数总和
X_Sum = np.sum(df.values)
# 频率矩阵
P = df.values/X_Sum
# 行边缘频率
Pr_margin = np.sum(P,1)
# 列边缘频率
Pc_margin = np.sum(P,0)
# 行边缘频率-1/2次方对角方阵
Dr_1_2 = np.diag(np.power(Pr_margin,-0.5))
# 列边缘频率-1/2次方对角方阵
Dc_1_2 = np.diag(np.power(Pc_margin,-0.5))
# 过渡矩阵
Z = np.zeros((P.shape))
for i in range(P.shape[0]):
for j in range(P.shape[1]):
Z[i][j] = (P[i][j] - Pr_margin[i]*Pc_margin[j])/np.power(Pr_margin[i]*Pc_margin[j],0.5)
print(Z)
[[ 0.06903397 0.01321385 0.00286337 -0.04560926 -0.07412414]
[ 0.00748675 0.0022116 0.00362348 0.00224728 -0.02129074]
[ 0.00335509 0.007487 -0.01807692 -0.01223391 0.02382773]
[-0.07387793 -0.02171255 0.01038691 0.04952401 0.06934967]]
卡方检验：
检验行变量和列变量之间(心理健康状况和社会经济状况)是否相互独立，若两者之间相互独立则不用进行对应分析说明两变量之间无关
.原假设H0：行变量与列变量相互独立
# 检验统计量chi2
chi2 = X.sum() * np.sum(np.power(Z,2))
print('chi2:{}'.format(chi2))
# 显著性水平alpha
alpha = 0.05
# 自由度freedom
freedoom = (Z.shape[0]-1) * (Z.shape[1]-1)
# 显著性水平alpha的临界值
chi2_alpha = stats.chi2.ppf(1-alpha*0.5,freedoom)
print('chi2_alpha:{}'.format(chi2_alpha))
chi2:45.594052238116525
chi2_alpha:23.33666415864534
因为chi2>chi_alpha，所以拒绝原假设，行和列变量之间不相互独立，可以进行对应分析
# 奇异值分解
U,Sigma,VT = np.linalg.svd(Z)
# 主惯量
inertia = np.power(Sigma,2)
# 主惯量贡献
inertia_contri = np.power(Sigma,2)/np.sum(inertia)
# 打印主惯量贡献表
table1 = pd.DataFrame(np.c_[Sigma,inertia,inertia_contri].T,
columns=['dim1','dim2','dim3','dim4'],
index=['Sigma','inertia','inertia_contri'])
print(table1)
# 主惯量贡献碎石图
plt.plot([1,2,3,4],inertia_contri,'r*-',linewidth = 2)
plt.xlabel('Strnge Num')
plt.ylabel('Inertia Contribute')
dim1 dim2 dim3 dim4
Sigma 0.161318 0.037088 0.008195 9.435155e-17
inertia 0.026024 0.001376 0.000067 8.902214e-33
inertia_contri 0.947475 0.050080 0.002445 3.241141e-31
从运碎石图可以看出第一维(dim1)贡献了94%的主惯量，第二维(dim2)贡献了5%的主惯量，使用两维作对应分析已经包含了原有数据99%的信息量，且第一维(dim1)占主要地位，对应分析时主要关注dim1。
# 心理健康状况荷载矩阵(取两维)
F = Dr_1_2.dot(U[:,0:2]).dot(np.diag(Sigma[0:2]))
# 心理健康状况标签
mental = ['good','slight','moderate','damaged']
# 社会经济水平荷载矩阵(取两维)
G = Dc_1_2.dot(VT.T[:,0:2]).dot(np.diag(Sigma[0:2]))
# 社会经济水平标签
economic = ['Level_A','Level_B','Level_C','Level_D','Level_E']
# 绘制对应分析图
plt.scatter(F[:,0],F[:,1],c='r',marker='*',linewidths=3)
for i in range(F.shape[0]):
plt.text(F[i,0],F[i,1],mental[i])
pass
plt.scatter(G[:,0],G[:,1],c='b',marker='o',linewidths=3)
for j in range(G.shape[0]):
plt.text(G[j,0],G[j,1],economic[j])
pass
plt.xlabel('dim1')
plt.ylabel('dim2')
print(F)
print(G)
[[-0.25966303 0.01327174]
[-0.02945593 0.02257353]
[ 0.01420067 -0.06995721]
[ 0.2372966 0.01969363]]
[[-0.1828982 -0.01552113]
[-0.05902637 -0.02244335]
[ 0.00888363 0.04233544]
[ 0.16540302 0.04332697]
[ 0.28768131 -0.06188021]]
细心的小伙伴可能已经看出来我们所求得荷载矩阵F等价于原案例中的-X1，我们所求的荷载矩阵G等价于原案例中的-Y1;
这是因为原案例中的U做了一个变换使U的每一维度(列)的数值之和>1(若U的某一列数值之和小于1则这一维度的数值取负，对应的V的列也取负，以使其满足Z=U.SIGMA.VT)
U是否做变换不影响对应分析，因为我们只关注相对位置的大小，等比例的放缩不影响相对关系
# 对U和VT进行变换
U_SUM = U.sum(axis=0)
for k in range(U_SUM.shape[0]):
if U_SUM[k] < 0:
U[:,k] = -U[:,k]
VT.T[:,k] = -VT.T[:,k]
# 心理健康状况荷载矩阵(取两维)
F = Dr_1_2.dot(U[:,0:2]).dot(np.diag(Sigma[0:2]))
# 心理健康状况标签
mental = ['good','slight','moderate','damaged']
# 社会经济水平荷载矩阵(取两维)
G = Dc_1_2.dot(VT.T[:,0:2]).dot(np.diag(Sigma[0:2]))
# 社会经济水平标签
economic = ['Level_A','Level_B','Level_C','Level_D','Level_E']
# 绘制对应分析图
plt.scatter(F[:,0],F[:,1],c='r',marker='*',linewidths=3)
for i in range(F.shape[0]):
plt.text(F[i,0],F[i,1],mental[i])
pass
plt.scatter(G[:,0],G[:,1],c='b',marker='o',linewidths=3)
for j in range(G.shape[0]):
plt.text(G[j,0],G[j,1],economic[j])
pass
plt.xlabel('dim1')
plt.ylabel('dim2')
print(F)
print(G)
[[ 0.25966303 -0.01327174]
[ 0.02945593 -0.02257353]
[-0.01420067 0.06995721]
[-0.2372966 -0.01969363]]
[[ 0.1828982 0.01552113]
[ 0.05902637 0.02244335]
[-0.00888363 -0.04233544]
[-0.16540302 -0.04332697]
[-0.28768131 0.06188021]]
简书的markdown语法支持真是差劲！
