対応分析（Correspondence Analysis）
A. Gifi (1991, reprinted with corrections), pp.265-271
対応分析のスクリプトを、Gifi(1991)の説明（pp.265-271）を基に作成してみた。以下に、備忘録として書き留めた。
対応表（Correspondence table；クロス表Cross table）
において、行
と行
の距離の２乗を、次式の
距離で与える。
式番号はGifi (1991)におけるものである。
は対応行列の
要素であり、
は第行の行和、
は第
列の列和である。
この距離は、次の値を第要素とする行列
の第行と第行のユークリッド距離の2乗である。
行列表記すると次式となる。
ここで、
はを対角要素とする対角行列であり、
はを対角要素とする対角行列である。
距離（8.15）に対応する距離を与える布置
は次式を満たすものとして与えることができる（Gifi, 1991, p.266）。
式（8.17）と式（a1）において、赤字で示した対応から、上式を満たすは次式の特異値分解から得ることが考えられる。
特異値分解については、＜岡本安晴、2009＞pp.47-49で説明している。
を次式で与える。
この値は式（a1）を満たすが、座標値が全て同じである次元も含まれているので、それを除いた解を得るために式（8.18）の左辺を次式のように変更したものの特異値分解を求める。
式（8.24）の左辺を
倍したときの各項の値は次のようになる。
ここでの要素の行と列が独立の仮定の下での期待値は次式で与えられることに注意。
従って、式（8.24）の左辺を
と置いたとき、
の要素は次式で与えられる。
行の座標に対して、列の座標
を、
の行の値による重み付きの平均値がになるように与える。すなわち
が成り立つようにする。このは次式で与えられる。
上の考え方に基づいて、リスト１に示す対応分析のスクリプト例を書いてみた。スクリプトファイルとサンプルデータは圧縮ファイルfiles.zipにまとめている。
なお、対応分析では対応表（クロス表）の行と列の扱いが対称ではないが、内積モデルでは対称に扱われている。
リスト１のスクリプトを、Gifi (1991)のデータ例（p. 269）に対して、以下のように実行してみた。
まず、Gifiのデータ例を図１のようにExcelファイルとして用意した。
図１
入力ファイルを用意した後、リスト１のスクリプトを実行する。
入力ファイル名を聞いてくるので、図１のファイルを入力ファイルとするときは、次のように入力する。
Excel Data
File(*.xlsx) = Gifi_p269.xlsx
入力ファイル名を設定して、Enterキーを押すと計算が始まり、図２のように求められた座標による布置が表示される。
図２
図２の布置は、Gifi (1991)のFigure 8.6のものと同じである。
図２のWindowを閉じると、次に表示する座標次元の入力が求められる。また、座標次元の入力メッセージの前に、布置のグラフがファイルに保存されたことも表示されているので、後で見ることができる。
Figdim1dim2.png was saved.
X-axis(<=2; if <= 0, then quit.)
= 2
Y-axis(<=2; if <= 0, then quit.)
= 1
上の設定例では、横軸が第２座標、縦軸が第１座標の設定になっているので、図３のような布置が表示される。
図３
表示座標として、０以下の値を設定すると、プログラムの実行終了である。
Figdim1dim2.png was saved.
X-axis(<=2; if <= 0, then quit.)
= 2
Y-axis(<=2; if <= 0, then quit.)
= 1
Figdim2dim1.png was saved.
X-axis(<=2; if <= 0, then quit.)
= -1
CA ended.
引用文献
Gifi, A. (1991). Nonlinear Multivariate Analysis, reprinted with corrections. John Wiley & Sons.
岡本安晴（2009） 統計学を学ぶための数学入門［下］．培風館
"""
Correspondence Analysis
cf. Gifi,1991. Nonlinear
Multivariate Analysis,
reprinted with correction
pp. 265-271
The script was written by Y.
Okamoto, 2022.12
"""
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
# ExcelファイルをPandasのDataFrame型として読み込む
fl_nm = input('Excel Data File(*.xlsx)
= ')
xlfile = pd.ExcelFile(fl_nm)
df = pd.read_excel(xlfile)
# １行目の値（列変数名）
col_names = df.keys()
# Excelの１行目
c00 = col_names[0]
# 行変数名の列のキー名
col_names = col_names[1:] # ２番目以降を取り出す
print(col_names)
# 行変数名
row_names = df[c00]
print(row_names)
# セル内の数値
df_data = df.values
# ２行目以降が取り出される
F = df_data[:,1:]
# 取り出した各行の２列目以降を取り出す
r = F.sum(axis = 1)
# 行単位の和
r = np.reshape(r, (len(r), 1)) # １次元配列から２次元配列（行列）への型変換
N = r.sum()
# データの総和
c = F.sum(axis = 0)
# 列単位の和
c = np.reshape(c, (len(c), 1)) # １次元配列から２次元配列（行列）への型変換
E = r @ c.T / N
# Expected
values
print('E =\n', E)
N_2_A = (F - E) / (E**0.5)
# The last line
in page 269..
print(N_2_A)
print('X^2 =', np.sum(N_2_A**2))
A = N_2_A / (N**0.5)
#
eq.(8.24)
print('A =\n', A)
# Aの特異値分解SVD
A = np.asfarray(A)
# Aの要素の型変換
K, Lmbd, Lh = np.linalg.svd(A,
full_matrices = False)
L = Lh.T
p = sum([1 if v > 1.0e-7 else 0 for
v in Lmbd]) # 布置の次元
print('p =', p)
Lmbd = Lmbd[:p]
# 0でない特異値
print('K =\n', K)
K = K[:,:p]
# 0でない特異値に対応する左特異ベクトル
print('K =\n', K)
print('Lmbd =', Lmbd)
print('Lmbd**2 =', Lmbd**2)
print('X^2 =', N * np.sum(Lmbd**2))
print('L =\n', L)
L = L[:,:p]
# 0でない特異値に対応する右特異ベクトル
print('L =\n', L)
X = ((N**0.5) *
np.diag(np.reshape(r**-0.5, (1, len(r)))[0])) \
@ K @
np.diag(Lmbd)
# eq.(8.19)
print('X =\n', X)
Y = ((N**0.5) * np.diag(np.reshape(c**-0.5,
(1, len(c)))[0])) @ L
print(Y)
XYp = X @ Y.T
# Projection of Y onto
X
print(XYp)
x_axis = 1
y_axis = 2 if p > 1 else 1
while True:
plt.plot(X.T[x_axis-1],
X.T[y_axis-1], 'o', label = 'Rows')
for i in
range(len(row_names)):
plt.text(X.T[x_axis-1][i]+0.01, X.T[y_axis-1][i]-0.01, row_names[i])
plt.plot(Y.T[x_axis-1],
Y.T[y_axis-1], 'o', label = 'Cols')
for i in
range(len(col_names)):
plt.text(Y.T[x_axis-1][i]+0.01, Y.T[y_axis-1][i]-0.01, col_names[i])
plt.legend()
plt.xlabel('Dim-{0:}'.format(x_axis), fontsize = 14)
plt.ylabel('Dim-{0:}'.format(y_axis), fontsize = 14)
plt.savefig(f'Figdim{x_axis}dim{y_axis}.png')
print(f'Figdim{x_axis}dim{y_axis}.png was saved.')
plt.show()
x_axis =
int(input('X-axis(<={0:}; if <= 0, then quit.) = '.format(p)))
if x_axis <= 0:
break
elif x_axis > p:
x_axis = p
y_axis =
int(input('Y-axis(<={0:}; if <= 0, then quit.) = '.format(p)))
if y_axis <= 0:
break
elif y_axis > p:
y_axis = p
print('CA ended.')