個人差尺度法（INDSCAL）
INDSCAL (Carroll &
Wish, 1974)では、個人
の空間における事物
と
の距離は次式で与えられる。
Carroll & Wish (1974)では、式（1）をCANDECOMP (CANonical DECOMPosition)の形
に変形して解の求め方が説明されている。CANDECOMPは、PARAFACあるいはCP (Canonical Polyadic decomposition)とも呼ばれている。式（２）の解法として、パラメータ群ごとに交互に最小2乗解を求める（ALS: Alternative Least Squaraes）交互最小２乗法が説明されている。INDSCALの解法について、Takane et al.(1977)において詳しく論じられている。
多次元尺度構成法において、直接、距離との差を極値探索法により最小にする布置を求める方法が試みられている（計量モデル、非計量モデル）。式（１）のモデルに対して交互にウェイト
と布置
の値を極値探索法により求めてみた。心理学の場合、距離データが直接与えられることは期待し難いので、距離と単調関係にあるデータ
が与えられたものとする。また、距離（１）を直接対象とする方法であるので、ユークリッド距離以外に簡単に拡張できるが、ここではユークリッド距離の場合について考える。極値探索法は、Powell法を用いたので、偏導関数などは不要である。
アルゴリズムは、以下のようなものを設定した。
与えられたデータをと置く。
（１）と単調関係にある距離の推定値
をKruskal’s up-and-down-block algorithm (Kruskal, 1964, pp. 127-128)と呼ばれている方法で求める。
（２）を固定して、式（１）の値
の2乗との2乗との差の２乗和が最小になるを極値探索法（Powell法）で求める。
（３）を固定して、式（１）の値の2乗との2乗との差の２乗和が最小になるを極値探索法（Powell法）で求める。
は、中心が原点、各座標の絶対値の最大値が１となるようにスケーリングを行う。のスケーリングに合わせてのスケーリングを行う。
（４）値が収束しているか調べる。収束していないときは、ステップ（２）に戻る。
（５）ステップ（４）において収束していると判定されたときは、現在のとの値に基づいての更新を行う。更新値に変化がある（収束していない）と判定されたときはステップ（２）に戻る。更新値に変化がない（収束している）と判定されたときは終了である。
上のアルゴリズムによるメインスクリプトをリスト１に示した。非計量単調変換（Kruskal’s up-and-down-block algorithm）のスクリプトはリスト２に、極値探索法の関数を利用するための1次元ベクトル関数を用意するクラス型のスクリプトはリスト３に示した。これらのファイルは圧縮ファイルindscalnonmetricfiles.zipとしてまとめた。
リスト3では、極値探索法の関数minimizeを用いるため、目的関数は１次元の配列（ベクトル）の関数としているが、１次元配列を２次元の配列に変換して、式（１）により距離を計算している。いる。この変換はリスト３に宣言されているクラス型class for_wでは、関数
def f(self, z):
"""
１次元配列（ベクトル）zを2次元配列wに並べ換える
"""
w =
np.resize(z, (self.p, self.r))
により、リスト３のクラス型class for_xでは、関数
def f(self, z):
"""
１次元配列（ベクトル）zを２次元配列xに並べ換える
"""
x =
np.resize(z, (self.n, self.r))
により、1次元配列zを２次元配列に変換している。
リスト１によるINDSCAL分析の入力データファイルは、CSVファイルとして用意する。例catdist.csvを図１に示す。
図１ catdist.csv
図１のデータは、事物A、B、C、D、E、F、G、H、Iの距離を、人P1、P2、・・・について表形式で設定したものである。
１行目の第２列から事物の名前を入力する。
第２行目には、１列目に人の名前を設定する。
第３行目から、先頭に事物の名前を記入した後、その事物と各事物との距離を設定していく。
すべての事物についてのデータを設定後、次の行の１列目に次の人の名前を設定する。名前の設定行の次の行から事物の距離を設定していく。
上の手順で、すべての人の事物についての距離のデータを設定する。
設定したデータの分析は、リスト１のスクリプトの実行時に入力データファイルとして指定する。
リスト１のスクリプトを実行すると、入力データファイル名の設定が求められる。次の例は、図１のファイル名を設定したものである。
Input File
Name (*.csv) = catdist.csv
入力データファイル名を設定すると、次に出力ファイル名の設定が求められるので、適当なテキストファイル名（ファイル拡張子が.txt）を設定する。
Input File Name (*.csv) = catdist.csv
Output File Name (*.txt) = Results.txt
出力ファイル名を設定すると、入力ファイルが読み込まれて分析が始まる。
まず、
の初期値が設定される。これは、距離データの平均値に対して２重中心化（岡本、2014、pp.197-200）により得られる座標値を初期値としている。の初期値はすべて１としている。
２重中心化された配列の固有分解の固有値が、正の値のみスクリープロットの形式で表示される（図２）。
図２
このスクリープロットから、布置の次元数の推測を行うことができる。図２のグラフの場合、次元数は２が適当であると思われる。
図２のWindowを閉じると、次元数を聞いてくる。次の例は、次元数を２と設定したものである。
Number of
dimensions(<=4) = 2
次元数を設定すると極値探索法の計算が始まる。解が収束するとShepard diagramが表示される（図３）。
図３ Shepard diagram
横軸に入力データ値、縦軸に推測距離をとり、個人別に折れ線で繋いだものである。
図３のWindowを閉じると、事物の布置が表示される（図４）。
図４
図４のWindowを閉じると、ウェイトの布置が表示される（図５）。
図５
図５のWindowを閉じると、端末に次の布置の描画における横軸（xdim）と縦軸（ydim）の設定を求めるメッセージが表示される。
Choose dimensions for plot. Set a
negative integer when you want to quit.
xdim(<=2) = -1
上の例では、１次元と２次元の布置の描画で充分であるので、布置の描画を終了するために「―１」を設定している。
―１を設定すると、スクリプトの実行は終了して、端末にはSTRESSの値が表示される。
Choose dimensions for plot. Set a
negative integer when you want to quit.
xdim(<=2) = -1
STRESS = 0.00000
スクリプトの実行終了後、出力ファイルを開くと以下のような内容になっている。
Input Data File = catdist.csv
Distances =
[[[1.0 2.0 4.0 2.0 3.0 4.0 4.0 4.0
5.0]
[2.0 1.0 2.0 3.0 2.0 3.0 4.0 4.0 4.0]
[4.0 2.0 1.0 4.0 3.0 2.0 5.0 4.0 4.0]
[2.0 3.0 4.0 1.0 2.0 4.0 2.0 3.0 4.0]
[3.0 2.0 3.0 2.0 1.0 2.0 3.0 2.0 3.0]
[4.0 3.0 2.0 4.0 2.0 1.0 4.0 3.0 2.0]
[4.0 4.0 5.0 2.0 3.0 4.0 1.0 2.0 4.0]
[4.0 4.0 4.0 3.0 2.0 3.0 2.0 1.0 2.0]
[5.0 4.0 4.0 4.0 3.0 2.0 4.0 2.0 1.0]]
[[1.0 3.0 5.0 1.0 3.0 5.0 1.0 3.0 5.0]
[3.0 1.0 3.0 3.0 1.0 3.0 3.0 1.0 3.0]
[5.0 3.0 1.0 5.0 3.0 1.0 5.0 3.0 1.0]
[1.0 3.0 5.0 1.0 3.0 5.0 1.0 3.0 5.0]
[3.0 1.0 3.0 3.0 1.0 3.0 3.0 1.0 3.0]
[5.0 3.0 1.0 5.0 3.0 1.0 5.0 3.0 1.0]
[1.0 3.0 5.0 1.0 3.0 5.0 1.0 3.0 5.0]
[3.0 1.0 3.0 3.0 1.0 3.0 3.0 1.0 3.0]
[5.0 3.0 1.0 5.0 3.0 1.0 5.0 3.0 1.0]]
[[1.0 1.0 1.0 3.0 3.0 3.0 5.0 5.0 5.0]
[1.0 1.0 1.0 3.0 3.0 3.0 5.0 5.0 5.0]
[1.0 1.0 1.0 3.0 3.0 3.0 5.0 5.0 5.0]
[3.0 3.0 3.0 1.0 1.0 1.0 3.0 3.0 3.0]
[3.0 3.0 3.0 1.0 1.0 1.0 3.0 3.0 3.0]
[3.0 3.0 3.0 1.0 1.0 1.0 3.0 3.0 3.0]
[5.0 5.0 5.0 3.0 3.0 3.0 1.0 1.0 1.0]
[5.0 5.0 5.0 3.0 3.0 3.0 1.0 1.0 1.0]
[5.0 5.0 5.0 3.0 3.0 3.0 1.0 1.0 1.0]]]
ck_diff_nm =
2.041340950954587e-21, i_nm_t = 2
Coordinates
dim-1
dim-2
A -0.96659 0.96650
B -0.00000 1.00000
C 0.96658 0.96649
D -1.00000 -0.00000
E -0.00000 -0.00001
F 0.99999 -0.00001
H -0.96663 -0.96649
I 0.00001 -1.00000
J 0.96664 -0.96649
Weights
dim-1
dim-2
P1 0.13382 0.13378
P2 0.24722 0.02036
P3 0.02044 0.24719
STRESS = 0.00000
先頭に、入力データファイル名が出力され、続いて読み込まれたデータ値が出力されている。
「Coordinates」に、事物の座標値xが出力され、続いて、重みwの値が出力されている。
最後に、STRESSの値が示されている。STRESSは、次式で算出されている。
Carroll, J. D. & Wish, M. (1974) Models and Methods for Three-way Multidimensional Scaling. In Krantz, D. H., Atkinson, R. C., Luce, R. D.,& Suppes, P. (1074) Contemporary Developments in Mathematical Psychology, Volume II: Measurement, psychophysics, and neural information processing. W. H. Freeman and Company.
Kruskal, J. B. (1964). Nonmetric multidimensional scaling: A numerical method. Psychometrika, 29, 115-129.
Takane, Y., Young, F. W., & deLeeuw, J. (1977) Nonmetric individual differences multidimensional scaling: An alternating least squares method with optimal scaling features. Psychometrika, 42, 7-67.
リスト１ INDSCALメインスクリプト（indscal_nonmetric.py：圧縮ファイルindscalnonmetricfiles.zipにまとめた）
"""
By Yasuharu Okamoto, 2024.11
"""
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import scipy.optimize as so
from myopt import *
# 極値探索用関数のクラス型
from BlockAverage_v1 import * # 非計量単調変換用
from copy import copy
fin_nm = input('Input File Name
(*.csv) = ')
with open(fin_nm, 'r') as f:
raw_data = pd.read_csv(f)
fout_nm = input('Output File Name
(*.txt) = ')
fout = open(fout_nm, 'wt')
fout.write('Input Data File = ' +
fin_nm + '\n')
obj_names = [v for v in
raw_data.keys()[1:]]
print(obj_names)
v_data = raw_data.values
print(v_data)
v_data = v_data[:,1:]
n = len(obj_names)
# 対象物数
n_ind = len(v_data) // (n + 1)
# 個人数
print('n =', n, ' n_ind =',n_ind)
pers_names = []
for i in range(n_ind):
pers_names.append(raw_data.values[i * (n + 1)][0])
print('Persons = ', pers_names)
del_ijk = []
for i in range(n_ind):
del_ijk.append(v_data[1 + (n
+ 1)*i:(n + 1)*(i+1)])
del_ijk = np.array(del_ijk)
print('Distances =\n', del_ijk)
fout.write(f'Distances
=\n{del_ijk}\n')
#
# 初期値設定用の平均距離の配列
#
m_dist = np.mean(del_ijk, axis=0)
D_2 = m_dist**2
J = np.diag([1]*n) - np.ones((n,n))/n
B = J @ (-D_2 / 2) @ J
# 2重中心化
B = np.array(B, dtype=float)
eig_val, eig_vec =
np.linalg.eigh(B)
# 内積行列の固有分解
ck_idx = eig_val > 1.0e-7
# 正の固有値の選択用配列添字
max_r = np.sum(ck_idx)
# 正の固有値の数
sortidx = np.argsort(-eig_val)
# 固有値の大きいものから降順に並べるための添字配列
eig_val = eig_val[sortidx]
eig_vec = eig_vec[:, sortidx]
#
# 正の固有値のスクリープロット
#
plt.title('Eigen Values(>1.0e-7)')
plt.plot(np.arange(1, max_r+1),
eig_val[:max_r])
plt.xticks(np.arange(1, max_r+1))
plt.ylabel('Eigen value')
plt.show()
r = int(input(f'Number of
dimensions(<={max_r}) = '))
if r > max_r:
r = max_r
print(f'Number of dimensions
replaced by {r}')
#
# 座標値xと重みw0の初期値
#
x = eig_vec[:,:r] * (eig_val[:r]**0.5)
w0 = np.ones((n_ind, r))
x = x - np.mean(x, axis=0)
# 中心化
s_x = np.max(np.abs(x), axis=0)
x = x / s_x
# 標準化
x_1 = x
# 初期値
w_1 = w0
# 初期値
d_ijk = copy(del_ijk)
i_nm_t = 0
# 非計量変換の回数
ck_d_hat = [0]*n_ind
ck_del = [0]*n_ind
x_g1 = x_1
w_g1 = w_1
while True:
i_nm_t += 1
x_g0 = x_g1
w_g0 = w_g1
i_trial = 0 # 非計量変換値に対する極値探索の回数
d_ijk1 = copy(d_ijk)
#
#
非計量変換
#
for s in range(n_ind):
d_delta = np.empty((n * (n - 1)) // 2)
t_pos = 0
for
i in range(1, n):
for j in range(i):
d_delta[t_pos] = del_ijk[s][i][j]
t_pos += 1
d_hat = np.empty((n * (n - 1) // 2))
i_pos = np.empty((n * (n - 1) // 2))
j_pos = np.empty((n * (n - 1) // 2))
t_pos = 0
for
i in range(1, n):
for j in range(i):
d_hat[t_pos] = d_ijk[s][i][j]
i_pos[t_pos] = i
j_pos[t_pos] = j
t_pos += 1
d_pair = np.array([(u, v, ipos, jpos) for u, v,ipos, jpos in
zip(d_delta, d_hat, i_pos, j_pos)],
dtype = ([('delta', float), ('dhat', float),
('ipos', int), ('jpos', int)]))
idx
= np.argsort(d_pair,order=['delta', 'dhat'])
dhat
= d_pair[idx]['dhat']
i_pair = d_pair[idx]['ipos']
j_pair = d_pair[idx]['jpos']
transf_d = nm_transf(dhat)
ck_del[s] = d_pair[idx]['delta']
ck_d_hat[s]
= transf_d
d_ijk1v = np.zeros((n, n))
for
p in range((n * (n - 1)) // 2):
d_ijk1v[i_pair[p]][j_pair[p]] = transf_d[p]
d_ijk1v[j_pair[p]][i_pair[p]] = transf_d[p]
d_ijk1v /= np.max(d_ijk1v)
for
j in range(n):
for k in range(n):
d_ijk[s][j][k] = d_ijk1v[j][k]
#
# 非計量変換値d_ijkに対する極値探索
#
while True:
i_trial +=1
x_0
= x_1
# 極値探索の初期値
w_0
= w_1
# 極値探索の初期値
opt_w = for_w(d_ijk**2, n, n_ind, r, x_0) # wの極値探索のためのクラス型
#
# wの極値探索
#
results = so.minimize(opt_w.f, np.resize(w_0, (n_ind*r)),
method='Powell')
print(results)
w = (np.reshape(results.x,
(n_ind, r))**2) # wの設定
opt_x = for_x(d_ijk**2, n, n_ind, r, w) # xの極値探索のためのクラス型
#
# xの極値探索
#
results = so.minimize(opt_x.f, np.resize(x_0, (n * r)), method='Powell')
print(results)
x =
np.reshape(results.x, (n, r))
# xの設定
x =
x - np.mean(x, axis=0)
# xの中心化
s_x
= np.max(np.abs(x), axis=0)
# 標準化のためのスケール値
x =
x / s_x
# 標準化
w =
w * (s_x**2)
# xの標準化に合わせた調整
print('x =\n', x)
print('w =\n', w)
x_1
= x # 極値探索による更新値
w_1
= w # 極値探索による更新値
ck_diff = np.sum((x_1 - x_0)**2) + np.sum((w_1 - w_0)**2) # 探索前後の変化量
print(f'ck_diff(at trial {i_trial}) =', ck_diff)
if
ck_diff < (1.0e-7): # 探索後の変化量が小さいので極値探索を終了する
break
if
i_trial > 50:
# 十分に探索を行ったとして終了する
break
#
# 更新値に対する距離の計算
#
for s in range(n_ind):
v =
x_1 * (w_1[s]**0.5)
for
j in range(n):
for k in range(n):
d_ijk1[s][j][k] = (np.sum((v[j] - v[k])**2))**0.5
v_max = np.max(d_ijk1[s])
d_ijk1[s] /= v_max
# 最大値を１に標準化
x_g1 = x_1
w_g1 = w_1
ck_xw = np.sum((x_g1 -
x_g0)**2) + np.sum((w_g1 - w_g0)**2)
print('ck_xw =', ck_xw,
' i_nm_t =', i_nm_t)
if ck_xw < 1.0e-5:
fout.write(f'\nck_xw = {ck_xw}, i_nm_t = {i_nm_t}\n')
break
ck_diff_nm = np.sum((d_ijk -
d_ijk1)**2)
print('\nck_diff_nm =',
ck_diff_nm, ' i_nm_t =', i_nm_t)
if ck_diff_nm < 1.0e-5:
fout.write(f'\nck_diff_nm = {ck_diff_nm}, i_nm_t = {i_nm_t}\n')
break
if i_nm_t > 20:
print('i_nm_t =', i_nm_t)
fout.write(f'\ni_nm_t = {i_nm_t}\n')
break
d_ijk = d_ijk1
#
# Nonmetric
trasformationの描画
#
for s in range(n_ind):
plt.plot(ck_del[s],
ck_d_hat[s], label=pers_names[s])
plt.xlabel(r'$\delta$', fontsize=14)
plt.ylabel(r'$\hat{d}$', fontsize=14)
plt.legend()
plt.show()
fout.write('\nCoordinates\n')
fout.write(f'{" ":>10s}')
for i in range(r):
sv = f'dim-{i+1}'
fout.write(f'{sv:>10s}')
fout.write('\n')
for i in range(n):
fout.write(f'{obj_names[i]:>10s}')
for j in range(r):
fout.write(f'{x_1[i][j]:>10.5f}')
fout.write('\n')
fout.write('\n\nWeights\n')
fout.write(f'{" ":>10s}')
for i in range(r):
sv = f'dim-{i+1}'
fout.write(f'{sv:>10s}')
fout.write('\n')
for i in range(n_ind):
fout.write(f'{pers_names[i]:>10s}')
for j in range(r):
fout.write(f'{w_1[i][j]:>10.5f}')
fout.write('\n')
#
# 解の布置の描画
#
xdim = 1
ydim = 2 if r >= 2 else xdim
while True:
#
# 項目の座標xの布置
#
xdim -= 1
ydim -= 1
plt.plot(x_1.T[xdim],
x_1.T[ydim], 'o')
for i in range(n):
plt.text(x_1.T[xdim][i], x_1.T[ydim][i], obj_names[i])
plt.xlabel(f'Dim-{xdim+1}')
plt.ylabel(f'Dim-{ydim+1}')
plt.title('Coordinates')
plt.show()
#
# 重みの座標wの布置
#
plt.plot(w_1.T[xdim],
w_1.T[ydim], 'o')
for i in range(n_ind):
plt.text(w_1.T[xdim][i], w_1.T[ydim][i], pers_names[i])
plt.xlabel(f'Dim-{xdim+1}')
plt.ylabel(f'Dim-{ydim+1}')
plt.title('Weights')
plt.show()
print('Choose dimensions for
plot. Set a negative integer when you want to quit.')
xdim =
int(input(f'xdim(<={r}) = '))
xdim = xdim if xdim <= r
else 1
if xdim < 1:
break
ydim =
int(input(f'ydim(<={r}) = '))
ydim = ydim if ydim <= r
else r
if ydim < 1:
break
#
# Stressの計算
#
sum_istress = 0.0
for i in range(n_ind):
dijk = np.empty((n,n))
xi = x_1 * (w_1[i]**0.5)
for j in range(n):
for
k in range(n):
dijk[j][k] = np.sum((xi[j] - xi[k])**2)**0.5
m_dijk = np.mean(dijk)
sum_den = np.sum((dijk -
m_dijk)**2)
sum_num = np.sum((dijk -
d_ijk[i])**2)
sum_istress += sum_num /
sum_den
stress = (sum_istress / n_ind) ** 0.5
print(f'STRESS = {stress:.5f}')
fout.write(f'\nSTRESS =
{stress:.5f}\n')
fout.close()
リスト２ 非計量単調変換（BlockAverage_v1.py）
import copy
def nm_transf(dist):
"""
The up-and-down-blocks algorithm
I.
Borg & P. J. F. Groenen(2005). Modern multidimensional scaling, 2nd ed.,
(p. 206).
J.
B. Kruskal (1964). Nonmetric multidimensional scaling.
Psychometrika, 115-129 (pp. 126-128).
"""
a_dist = copy.copy(dist)
n = len(a_dist)
#
# Check
a_dist[0],...,a_dist[n-1] for monotonicity
#
Lpos = 0 # Start from a_dist[0]
Rpos = 0
transf_done = False
while transf_done == False:
#
# Check
a_dist[Lpos],...,a_dist[Rpos] for monotonicity
#
Lok
= False
Rok
= False
#
# Check the
right value a_dist[Rpos]
# If necessary,
adjust the values
#
if
Rpos >= n - 1:
Rok
= True
elif
a_dist[Rpos] <= a_dist[Rpos+1]:
Rok = True
else:
v = a_dist[Rpos] * (Rpos - Lpos + 1) + a_dist[Rpos+1]
v /= Rpos - Lpos + 2
Rpos += 1
for i in range(Lpos, Rpos+1):
a_dist[i] = v
#
# Check
the left value a_dist[Lpos]
# If
necessary, adjust the values
#
if
Lpos <= 0:
Lok = True
elif
a_dist[Lpos-1] <= a_dist[Lpos]:
Lok = True
else:
v = a_dist[Lpos] * (Rpos - Lpos + 1) + a_dist[Lpos-1]
v /= (Rpos - Lpos + 2)
Lpos -= 1
for i in range(Lpos, Rpos+1):
a_dist[i] = v
if
Lok and Rok: # Monotonicity is satisfied. No adjustment
is done
Rpos += 1
Lpos = Rpos
if
Rpos > n-1:
transf_done = True
return a_dist
リスト３ 極値探索法を用いるためのクラス型（myopt.py）
import numpy as np
class for_w:
"""
wの極値探索のためのクラス型
"""
def __init__(self, dist2, n,
p, r, x):
self.dist2 = dist2
self.n = n
self.p = p
self.r = r
self.x = x
# xは固定
def f(self, z):
"""
１次元配列（ベクトル）zを2次元配列wに並べ換える
"""
w =
np.resize(z, (self.p, self.r))
w =
w**2
d2_hat = []
for
s in range(self.p):
d2_temp = np.zeros((self.n,self.n))
v = self.x * (w[s]**0.5)
for i in range(self.n):
for j in range(self.n):
d2_temp[i][j] = np.sum((v[i] - v[j])**2)
d2_hat.append(d2_temp)
d2_hat = np.array(d2_hat)
err
= np.sum((self.dist2 - d2_hat)**2)
return err
class for_x:
"""
xの極値探索のためのクラス型
"""
def __init__(self, dist2, n,
p, r, w):
self.dist2 = dist2
self.n = n
self.p = p
self.r = r
self.w = w
# wは固定
def f(self, z):
"""
１次元配列（ベクトル）zを２次元配列xに並べ換える
"""
x =
np.resize(z, (self.n, self.r))
d2_hat = []
for
s in range(self.p):
d2_temp = np.zeros((self.n,self.n))
v = x * (self.w[s]**0.5)
for i in range(self.n):
for j in range(self.n):
d2_temp[i][j] = np.sum((v[i] - v[j])**2)
d2_hat.append(d2_temp)
d2_hat = np.array(d2_hat)
err
= np.sum((self.dist2 - d2_hat)**2)
return err