案例：
?將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：行變量與列變量相互獨立
# 檢驗統(tǒng)計量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語法支持真是差勁虱咧！
