「統計学実践ワークブック」 演習問題etc Ch.26 「その他の多変量解析手法」

当記事は「統計学実践ワークブック(学術図書出版社)」の読解サポートにあたってChapter.26の「その他の多変量解析手法」に関して演習問題を中心に解説を行います。多次元尺度構成法や数量化法、対応分析は用いられることがあるので、演習を通して抑えておくと良いと思われました。

本章のまとめ

多次元尺度構成法

演習問題解説

問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