当記事は「統計学実践ワークブック(学術図書出版社)」の読解サポートにあたってChapter.26の「その他の多変量解析手法」に関して演習問題を中心に解説を行います。多次元尺度構成法や数量化法、対応分析は用いられることがあるので、演習を通して抑えておくと良いと思われました。
本章のまとめ
多次元尺度構成法
下記で詳しく取り扱った。
https://www.hello-statisticians.com/explain-terms-cat/mds1.html
演習問題解説
問26.1
$[1]$
$$
\large
\begin{align}
B = – \frac{1}{2} \left( \mathit{I}_{n} – \frac{1}{n} \mathit{J}_n \right) D \left( \mathit{I}_{n} – \frac{1}{n} \mathit{J}_n \right)
\end{align}
$$
上記の式に基づいて$B$の計算を行う。
import numpy as np
I_4 = np.eye(4)
J_4 = np.ones([4,4])
D = np.array([[0., 5., 5., 16.], [5., 0., 4., 5.], [5., 4., 0., 5.], [16., 5., 5., 0.]])
B = -np.dot(np.dot(I_4-J_4/4., D), I_4-J_4/4.)/2.
lamb, eig_vec = linalg.eig(B)
print("B:")
print(B)
print("lambda: {}".format(lamb.real))
print("eigen vector:")
print(eig_vec.real)・実行結果
B:
[[ 4. -0. -0. -4.]
[-0. 1. -1. -0.]
[-0. -1. 1. -0.]
[-4. -0. -0. 4.]]
lambda: [ 8. 0. 2. 0.]
eigen vector:
array([[ 0.70710678, -0.70710678, 0. , 0. ],
[ 0. , 0. , -0.70710678, -0.70710678],
[ 0. , 0. , 0.70710678, -0.70710678],
[-0.70710678, -0.70710678, 0. , 0. ]])$[2]$
下記のようにそれぞれ計算できる。
import numpy as np
from scipy import linalg
I_4 = np.eye(4)
J_4 = np.ones([4,4])
D = np.array([[0., 5., 5., 16.], [5., 0., 4., 5.], [5., 4., 0., 5.], [16., 5., 5., 0.]])
B = -np.dot(np.dot(I_4-J_4/4., D), I_4-J_4/4.)/2.
lamb, eig_vec = linalg.eig(B)
x = eig_vec.real[:,0]*np.sqrt(lamb.real[0])
x_row = np.repeat(x,4).reshape([4,4])
x_col = np.repeat(x,4).reshape([4,4]).T
D = (x_row-x_col)**2
print("each x: {}".format(x))
print("D:")
print(D)・実行結果
each x: [ 2. 0. 0. -2.]
D:
[[ 0. 4. 4. 16.]
[ 4. 0. 0. 4.]
[ 4. 0. 0. 4.]
[ 16. 4. 4. 0.]]上記で計算したDでは元々の$D$を完全には再現されていないことが確認できる。
$[3]$
下記のようにそれぞれ計算できる。
import numpy as np
from scipy import linalg
I_4 = np.eye(4)
J_4 = np.ones([4,4])
D = np.array([[0., 5., 5., 16.], [5., 0., 4., 5.], [5., 4., 0., 5.], [16., 5., 5., 0.]])
B = -np.dot(np.dot(I_4-J_4/4., D), I_4-J_4/4.)/2.
lamb, eig_vec = linalg.eig(B)
x = np.dot(eig_vec.real[:,np.array([0,2])], np.array([[np.sqrt(lamb.real[0]), 0.], [0., np.sqrt(lamb.real[2])]]))
D = np.zeros([4,4])
for i in range(4):
for j in range(4):
D[i,j] = np.sum((x[i,:]-x[j,:])**2)
print("each x:")
print(x)
print("D:")
print(D)・実行結果
each x:
[[ 2. 0.]
[ 0. -1.]
[ 0. 1.]
[-2. 0.]]
D:
[[ 0. 5. 5. 16.]
[ 5. 0. 4. 5.]
[ 5. 4. 0. 5.]
[ 16. 5. 5. 0.]]xの$2$行目と$3$行目がワークブックの解答とは逆だが、固有ベクトルの向きが逆であるだけであり、該当する指定が特にないので、上記の結果も正しいと考えることができる。また、上記で計算したDでは元々の$D$が完全に再現されていることも確認できる。
参考
・準1級関連まとめ
https://www.hello-statisticians.com/toukeikentei-semi1