ミンコフスキー距離による個人差尺度法
事物間の心理的距離は、各次元での差異が加法的に処理されるcity block距離、最も差異の大きい次元によって事物間の距離が決まるsup-metric、各次元における差異が総合的に評価されて距離が決まるユークリッド距離などがあると説明されることがある。これらの距離により単位円を描くと、図１のようになる。
図１ 単位円（リスト０のスクリプトで描いた）
これらの距離は、Minkowski distance (metric)の特殊な場合である（Borg & Groenen, 2005）。Minkowski 距離による個人差尺度法（Individual Differences Scaling: INDSCAL；Carroll & Wish, 1974)は、次式で与えられる。
のときがcity block距離、
のときがユークリッド距離、
のときがsup-metricである。のとき、式（1）は次式となる。
ユークリッド距離の場合は、式（1）をCANDECOMP (CANonical DECOMPosition)の形
に変形して解の求める方法が説明されている。CANDECOMPは、PARAFACあるいはCP (Canonical Polyadic decomposition)とも呼ばれている。式（3）の解法として、パラメータ群ごとに交互に最小2乗解を求める（ALS: Alternative Least Squaraes）交互最小２乗法が説明されている。INDSCALの解法について、Takane et al.(1977)において詳しく論じられている。
多次元尺度法において、距離に対して直接極値探索法を適用することができる。式（１）あるいは式（２）に対して直接極値探索法の適用を試みる。
また、与えられたデータ
は距離と単調関係にあると考えられるので、から距離の推定値
をKruskal’s up-and-down-block algorithm (Kruskal, 1964, pp. 127-128)と呼ばれている方法で求める。スクリプトをリスト１に示した。これらのファイルは、圧縮ファイルindscalminkowskifiiles.zipにまとめた。
式（１）あるいは式（２）のモデルに対して交互にウェイト
と布置
の値を極値探索法により求める。極値探索法は、Powell法を用いるので、偏導関数などは不要である。
アルゴリズムは、以下のようなものを設定した。
与えられたデータをと置く。
（１）と単調関係にある距離の推定値をKruskal’s up-and-down-block algorithm (Kruskal, 1964, pp. 127-128)と呼ばれている方法で求める。
（２）を固定して、式（１）あるいは式（２）の値
との差の２乗和が最小になるを極値探索法（Powell法）で求める。
（３）を固定して、式（１）あるいは式（２）の値との差の２乗和が最小になるを極値探索法（Powell法）で求める。
は、中心が原点、各座標の絶対値の最大値が１となるようにスケーリングを行う。のスケーリングに合わせてのスケーリングを行う。
（４）値が収束しているか調べる。収束していないときは、ステップ（２）に戻る。
（５）ステップ（４）において収束していると判定されたときは、現在のとの値に基づいての更新を行う。更新値に変化がある（収束していない）と判定されたときはステップ（２）に戻る。更新値に変化がない（収束している）と判定されたときは終了である。
距離の式が式（１）と式（２）の2通りあるので、スクリプトを式（１）用と式（２）のsup-metric用の2通りを用意した。以下に、それぞれについて説明する。
式（１）の場合：city block距離（
）、ユークリッド距離（
）など
式（１）用のスクリプトをリスト２に示す。極値探索法では、1次元の配列と２次元の配列の型の変換を行って極値探索法の関数scipy.optimize.minimizeを用いているが、そのためのクラス型をリスト３に宣言している。
リスト2による分析の入力データファイルは、CSVファイルとして用意する。例cat9city.csvを図2に示す。９ポイントカテゴリーによる仮想評定データである。
図２ cat9city.csv
図２のデータは、事物A、B、C、D、E、F、G、H、Iの距離を、人P1、P2、・・・について表形式で設定したものである。
１行目の第２列から事物の名前を入力する。
第２行目には、１列目に人の名前を設定する。
第３行目から、先頭に事物の名前を記入した後、その事物と各事物との距離を設定していく。
すべての事物についてのデータを設定後、次の行の１列目に次の人の名前を設定する。名前の設定行の次の行から事物の距離を設定していく。
上の手順で、すべての人の事物についての距離のデータを設定する。
リスト２のスクリプトを実行すると、入力データファイル名の設定が求められる。次の例では、図２のファイル名cat9city.csvが設定されている。入力データファイル名を設定すると、出力ファイル名の設定が求められるので、適当なテキストファイル名を設定する。次の例では、results.txtと設定されている。
Input File Name (*.csv) = cat9city.csv
Output File Name (*.txt) = results.txt
出力ファイル名を設定すると、次にMinkowski距離の
の値の設定が求められる。City block距離を用いる場合は、１と設定すればよい。
Minkowski
distance: p(>=1.0) = 1
の値を設定すると計算が始まる。座標の初期値を設定するためにデータから平均距離の配列が計算され、２重中心化（岡本、2014、pp.197-200）により得られた行列の固有値の正の値のもののスクリープロットが表示される（図３）。
図３
図３のWindowを閉じると、事物の布置を求めるときの次元数の設定が求められる。データの平均値に対して算出された固有値の分布は、事物の布置の大まかな空間構造を反映していると考えられる。図３より空間の次元として２が適当であると推測し、次の例では、２を設定している。
Number of
dimensions(<=5) = 2
次元数を設定すると、極値探索の繰り返しが始まり、途中経過が端末に表示される。極値探索が終了すると、横軸に
、縦軸に
をとったShepard diagramが個人別に表示される（図４）。
図４
図４のWindowを閉じると、事物の布置のグラフが表示される（図５）。
図５
図５のWindowを閉じると、重み
のグラフが表示される（図６）。
図６
図６のwindowを閉じると、次に表示する布置と重みのグラフの横軸（X軸）の設定が求められる。次元として負の値を設定するとグラフ表示は終了する。次の例では、グラフ表示を終了するために、―１を設定している。
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
STRESSは、次式で算出されている。
実行開始時に指定した出力ファイルを開くと、以下のような内容である。
Input Data File = cat9city.csv
p = 1.0
Distances =
[[[1.0 9.0 4.0 5.0 9.0 6.0 6.0 5.0
6.0]
[9.0 1.0 6.0 8.0 1.0 4.0 8.0 9.0 6.0]
[4.0 6.0 1.0 6.0 6.0 3.0 8.0 7.0 2.0]
[5.0 8.0 6.0 1.0 7.0 4.0 2.0 1.0 7.0]
[9.0 1.0 6.0 7.0 1.0 3.0 7.0 8.0 7.0]
[6.0 4.0 3.0 4.0 3.0 1.0 5.0 5.0 4.0]
[6.0 8.0 8.0 2.0 7.0 5.0 1.0 2.0 9.0]
[5.0 9.0 7.0 1.0 8.0 5.0 2.0 1.0 8.0]
[6.0 6.0 2.0 7.0 7.0 4.0 9.0 8.0 1.0]]
[[1.0 9.0 3.0 4.0 9.0 5.0 6.0 4.0 4.0]
[9.0 1.0 7.0 7.0 1.0 5.0 6.0 7.0 7.0]
[3.0 7.0 1.0 4.0 7.0 4.0 6.0 4.0 2.0]
[4.0 7.0 4.0 1.0 6.0 2.0 2.0 1.0 4.0]
[9.0 1.0 7.0 6.0 1.0 4.0 5.0 7.0 7.0]
[5.0 5.0 4.0 2.0 4.0 1.0 3.0 3.0 4.0]
[6.0 6.0 6.0 2.0 5.0 3.0 1.0 2.0 6.0]
[4.0 7.0 4.0 1.0 7.0 3.0 2.0 1.0 4.0]
[4.0 7.0 2.0 4.0 7.0 4.0 6.0 4.0 1.0]]
[[1.0 7.0 4.0 3.0 6.0 4.0 5.0 4.0 6.0]
[7.0 1.0 3.0 8.0 2.0 3.0 8.0 9.0 4.0]
[4.0 3.0 1.0 7.0 3.0 2.0 8.0 7.0 2.0]
[3.0 8.0 7.0 1.0 7.0 5.0 2.0 1.0 8.0]
[6.0 2.0 3.0 7.0 1.0 2.0 7.0 8.0 4.0]
[4.0 3.0 2.0 5.0 2.0 1.0 6.0 6.0 4.0]
[5.0 8.0 8.0 2.0 7.0 6.0 1.0 1.0 9.0]
[4.0 9.0 7.0 1.0 8.0 6.0 1.0 1.0 9.0]
[6.0 4.0 2.0 8.0 4.0 4.0 9.0 9.0 1.0]]]
ck_diff_nm =
8.57876027139738e-22, i_nm_t = 2
Coordinates
dim-1
dim-2
A -0.28357 0.92102
B 0.66163 -1.00000
C 0.56691 0.63181
D -0.74216 0.08574
E 0.42134 -0.93546
F 0.32286 -0.08148
G -0.98889 -0.25875
H -0.95811 0.17240
I 1.00000 0.46472
Weights
dim-1
dim-2
P1 0.34443 0.35109
P2 0.15301 0.44527
P3 0.44338 0.16334
STRESS = 0.00000
読み込まれたデータ、算出された布置と重みなどが出力されている。
基本的には式（１）用のものと同じである。
式（２）用のスクリプトをリスト４に示す。極値探索法では、1次元の配列と２次元の配列の型の変換を行って極値探索法の関数scipy.optimize.minimizeを用いているが、そのためのクラス型をリスト５に宣言している。
リスト４による分析の入力データファイルは、CSVファイルとして用意する。例cat9sup.csvを図S1に示す。９ポイントカテゴリーによる仮想評定データである。
図S1 cat9sup.csv
図S1のデータは、事物A、B、C、D、E、F、G、H、Iの距離を、人P1、P2、・・・について表形式で設定したものである。
１行目の第２列から事物の名前を入力する。
第２行目には、１列目に人の名前を設定する。
第３行目から、先頭に事物の名前を記入した後、その事物と各事物との距離を設定していく。
すべての事物についてのデータを設定後、次の行の１列目に次の人の名前を設定する。名前の設定行の次の行から事物の距離を設定していく。
上の手順で、すべての人の事物についての距離のデータを設定する。
リスト４のスクリプトを実行すると、入力データファイル名の設定が求められる。次の例では、図S1のファイル名cat9sup.csvが設定されている。入力データファイル名を設定すると、出力ファイル名の設定が求められるので、適当なテキストファイル名を設定する。次の例では、results.txtと設定されている。
Input File Name (*.csv) = cat9sup.csv
Output File Name (*.txt) = results.txt
出力ファイル名を設定すると、計算が始まる。座標の初期値を設定するためにデータから平均距離の配列が計算され、２重中心化（岡本、2014、pp.197-200）により得られた行列の固有値の正の値のもののスクリープロットが表示される（図S2）。
図S2
図S2のWindowを閉じると、事物の布置を求めるときの次元数の設定が求められる。データの平均値に対して算出された固有値の分布は、事物の布置の大まかな空間構造を反映していると考えられる。図S2より空間の次元として２が適当であると推測し、次の例では、２を設定している。
Number of
dimensions(<=6) = 2
次元数を設定すると、極値探索における繰り返しが始まり、途中経過が端末に表示される。極値探索が終了すると、横軸に、縦軸にをとったShepard diagramが個人別に表示される（図S3）。
図S3
図S3のWindowを閉じると、事物の布置のグラフが表示される（図S4）。
図S4
図S4のWindowを閉じると、重みのグラフが表示される（図S5）。
図S5
図S5のwindowを閉じると、次に表示する布置と重みのグラフの横軸（X軸）の設定が求められる。次元として負の値を設定するとグラフ表示は終了する。次の例では、グラフ表示を終了するために、―１を設定している。
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.00246
STRESSは、次式で算出されている。
実行開始時に指定した出力ファイルを開くと、以下のような内容である。
Input Data File = cat9sup.csv
Sup-metric
Distances =
[[[1.0 9.0 4.0 4.0 9.0 5.0 5.0 4.0
6.0]
[9.0 1.0 8.0 7.0 2.0 5.0 8.0 8.0 7.0]
[4.0 8.0 1.0 6.0 8.0 4.0 8.0 7.0 2.0]
[4.0 7.0 6.0 1.0 6.0 6.0 2.0 1.0 8.0]
[9.0 2.0 8.0 6.0 1.0 4.0 7.0 7.0 7.0]
[5.0 5.0 4.0 6.0 4.0 1.0 7.0 6.0 3.0]
[5.0 8.0 8.0 2.0 7.0 7.0 1.0 2.0 9.0]
[4.0 8.0 7.0 1.0 7.0 6.0 2.0 1.0 9.0]
[6.0 7.0 2.0 8.0 7.0 3.0 9.0 9.0 1.0]]
[[1.0 9.0 2.0 4.0 9.0 5.0 6.0 4.0 2.0]
[9.0 1.0 8.0 6.0 1.0 5.0 4.0 6.0 8.0]
[2.0 8.0 1.0 3.0 8.0 4.0 4.0 3.0 1.0]
[4.0 6.0 3.0 1.0 5.0 2.0 2.0 1.0 3.0]
[9.0 1.0 8.0 5.0 1.0 5.0 4.0 6.0 7.0]
[5.0 5.0 4.0 2.0 5.0 1.0 2.0 2.0 3.0]
[6.0 4.0 4.0 2.0 4.0 2.0 1.0 2.0 4.0]
[4.0 6.0 3.0 1.0 6.0 2.0 2.0 1.0 3.0]
[2.0 8.0 1.0 3.0 7.0 3.0 4.0 3.0 1.0]]
[[1.0 5.0 4.0 3.0 4.0 4.0 4.0 3.0 6.0]
[5.0 1.0 3.0 7.0 2.0 2.0 8.0 8.0 3.0]
[4.0 3.0 1.0 6.0 3.0 1.0 8.0 7.0 2.0]
[3.0 7.0 6.0 1.0 6.0 6.0 2.0 1.0 8.0]
[4.0 2.0 3.0 6.0 1.0 2.0 7.0 7.0 3.0]
[4.0 2.0 1.0 6.0 2.0 1.0 7.0 6.0 3.0]
[4.0 8.0 8.0 2.0 7.0 7.0 1.0 1.0 9.0]
[3.0 8.0 7.0 1.0 7.0 6.0 1.0 1.0 9.0]
[6.0 3.0 2.0 8.0 3.0 3.0 9.0 9.0 1.0]]]
ck_xw = 5.030582809972423e-06, i_nm_t = 3
Coordinates
dim-1
dim-2
A -0.28409 -0.87977
B 0.70291 1.00000
C 0.52365 -0.63627
D -0.73294 -0.08939
E 0.41428 0.95053
F 0.38731 0.08043
G -1.00000 0.23414
H -0.90302 -0.17093
I 0.89190 -0.48875
Weights
dim-1
dim-2
P1 0.52880 0.52420
P2 0.18475 0.53201
P3 0.52883 0.21219
STRESS = 0.00246
読み込まれたデータ、算出された布置と重みなどが出力されている。
Borg, I. & Groenen, P. J. F. (2005). Modern multidimensional scaling: Theory and applications, 2nd ed. Springer.
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.
import numpy as np
import matplotlib.pyplot as plt
plt.figure(figsize=(7,7))
sup_coord =
np.array([[-1,-1],[1,-1],[1,1],[-1,1],[-1,-1]])
plt.plot(sup_coord.T[0],
sup_coord.T[1], 'b', label='Sup')
city_coord = np.array([[-1,0],[0,-1],
[1,0],[0,1],[-1,0]])
plt.plot(city_coord.T[0],
city_coord.T[1], 'g', label='City')
euc_step = np.linspace(0.0, 2*np.pi,
10000)
euc_x = np.cos(euc_step)
euc_y = np.sin(euc_step)
plt.plot(euc_x, euc_y, 'r',
label='Euclid')
plt.plot([-1.1, 1.1], [0, 0], 'k',
lw=1)
plt.plot([0, 0],[-1.1, 1.1], 'k',
lw=1)
plt.legend(loc='center', fontsize=16)
plt.show()
リスト１ 非計量単調変換（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
リスト２ 式（１）用のスクリプト（indscal_minkowski.py）
"""
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_minkowski 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')
p = float(input('Minkowski distance:
p(>=1.0) = '))
fout.write(f'\np = {p}\n')
obj_names = [v for v in
raw_data.keys()[1:]]
print(obj_names)
v_data = raw_data.values
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_n_dim = 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_n_dim+1),
eig_val[:max_n_dim])
plt.xticks(np.arange(1, max_n_dim+1))
plt.ylabel('Eigen value')
plt.show()
n_dim = int(input(f'Number of
dimensions(<={max_n_dim}) = '))
if n_dim > max_n_dim:
n_dim = max_n_dim
print(f'Number of dimensions
replaced by {n_dim}')
#
# 座標値xと重みw0の初期値
#
x = eig_vec[:,:n_dim] *
(eig_val[:n_dim]**0.5)
w0 = np.ones((n_ind, n_dim))
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
ip in range((n * (n - 1)) // 2):
d_ijk1v[i_pair[ip]][j_pair[ip]] = transf_d[ip]
d_ijk1v[j_pair[ip]][i_pair[ip]] = transf_d[ip]
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, n, n_ind, n_dim, p, x_0) # wの極値探索のためのクラス型
#
# wの極値探索
#
results = so.minimize(opt_w.f, np.resize(w_0, (n_ind * n_dim)),
method='Powell')
print(results)
w =
(np.reshape(results.x, (n_ind, n_dim))**2) # wの設定
opt_x = for_x(d_ijk, n, n_ind, n_dim, p, w) # xの極値探索のためのクラス型
#
# xの極値探索
#
results = so.minimize(opt_x.f, np.resize(x_0, (n * n_dim)),
method='Powell')
print(results)
#print('check_w =\n', w)
x =
np.reshape(results.x, (n, n_dim))
# 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**p)
# 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]**(1/p))
for
j in range(n):
for k in range(n):
d_ijk1[s][j][k] = (np.sum((np.abs(v[j] - v[k]))**p))**(1/p)
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(n_dim):
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(n_dim):
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(n_dim):
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(n_dim):
fout.write(f'{w_1[i][j]:>10.5f}')
fout.write('\n')
#
# 解の布置の描画
#
xdim = 1
ydim = 2 if n_dim >= 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(<={n_dim})
= '))
xdim = xdim if xdim <=
n_dim else 1
if xdim < 1:
break
ydim =
int(input(f'ydim(<={n_dim}) = '))
ydim = ydim if ydim <=
n_dim else n_dim
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]**(1/p))
for j in range(n):
for
k in range(n):
dijk[j][k] = np.sum((np.abs(xi[j] - xi[k]))**p)**(1/p)
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()
リスト３ 配列の型変換を行って極値探索法のための関数を用意するクラス型の宣言（myopt_minkowski.py）
import numpy as np
class for_w:
"""
wの極値探索のためのクラス型
"""
def __init__(self, dist, n,
n_pers, n_dim, p, x):
self.dist = dist
self.n = n
self.n_pers = n_pers
self.n_dim = n_dim
self.p = p
self.x = x
# xは固定
def f(self, z):
"""
１次元配列（ベクトル）zを2次元配列wに並べ換える
"""
w =
np.resize(z, (self.n_pers, self.n_dim))
w =
w**2
d_hat = []
for
s in range(self.n_pers):
d_temp = np.zeros((self.n,self.n))
v = self.x * (w[s]**(1/self.p))
for i in range(self.n):
for j in range(self.n):
d_temp[i][j] = (np.sum((np.abs(v[i] - v[j]))**self.p)) ** \
(1.0/self.p)
d_hat.append(d_temp)
d_hat = np.array(d_hat)
err
= np.sum((self.dist - d_hat)**2)
return err
class for_x:
"""
xの極値探索のためのクラス型
"""
def __init__(self, dist, n,
n_pers, n_dim, p, w):
self.dist = dist
self.n = n
self.n_pers = n_pers
self.n_dim = n_dim
self.p = p
self.w = w
# wは固定
#print('self.p =', self.p)
def f(self, z):
"""
１次元配列（ベクトル）zを２次元配列xに並べ換える
"""
x =
np.resize(z, (self.n, self.n_dim))
d_hat = []
for
s in range(self.n_pers):
d_temp = np.zeros((self.n,self.n))
v = x * (self.w[s]**(1/self.p))
for i in range(self.n):
for j in range(self.n):
d_temp[i][j] = (np.sum((np.abs(v[i] - v[j]))**self.p)) ** \
(1.0/self.p)
d_hat.append(d_temp)
d_hat = np.array(d_hat)
err
= np.sum((self.dist - d_hat)**2)
return err
リスト４ Sup-metric用メインスクリプト（indscal_sup.py）
"""
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_sup 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')
fout.write('\nSup-metric\n')
obj_names = [v for v in
raw_data.keys()[1:]]
print(obj_names)
v_data = raw_data.values
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_n_dim = 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_n_dim+1),
eig_val[:max_n_dim])
plt.xticks(np.arange(1, max_n_dim+1))
plt.ylabel('Eigen value')
plt.show()
n_dim = int(input(f'Number of
dimensions(<={max_n_dim}) = '))
if n_dim > max_n_dim:
n_dim = max_n_dim
print(f'Number of dimensions
replaced by {n_dim}')
#
# 座標値xと重みw0の初期値
#
x = eig_vec[:,:n_dim] *
(eig_val[:n_dim]**0.5)
w0 = np.ones((n_ind, n_dim))
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
ip in range((n * (n - 1)) // 2):
d_ijk1v[i_pair[ip]][j_pair[ip]] = transf_d[ip]
d_ijk1v[j_pair[ip]][i_pair[ip]] = transf_d[ip]
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
# 極値探索の初期値
print('x_0 =\n', x_0)
opt_w = for_w(d_ijk, n, n_ind, n_dim, x_0) # wの極値探索のためのクラス型
#
# wの極値探索
#
results = so.minimize(opt_w.f, np.resize(w_0, (n_ind * n_dim)),
method='Powell')
print(results)
w =
(np.reshape(results.x, (n_ind, n_dim))**2) # wの設定
opt_x = for_x(d_ijk, n, n_ind, n_dim, w) # xの極値探索のためのクラス型
#
# xの極値探索
#
results = so.minimize(opt_x.f, np.resize(x_0, (n * n_dim)),
method='Powell')
print(results)
x =
np.reshape(results.x, (n, n_dim))
# 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
# 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]
for
j in range(n):
for k in range(n):
d_ijk1[s][j][k] = np.max(np.abs(v[j] - v[k]))
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(n_dim):
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(n_dim):
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(n_dim):
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(n_dim):
fout.write(f'{w_1[i][j]:>10.5f}')
fout.write('\n')
#
# 解の布置の描画
#
xdim = 1
ydim = 2 if n_dim >= 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(<={n_dim}) = '))
xdim = xdim if xdim <=
n_dim else 1
if xdim < 1:
break
ydim =
int(input(f'ydim(<={n_dim}) = '))
ydim = ydim if ydim <=
n_dim else n_dim
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]
for j in range(n):
for
k in range(n):
dijk[j][k] = np.max(np.abs(xi[j] - xi[k]))
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()
リスト５ 極値探索法ための配列の型変換用クラス型の宣言（myopt_sup.py）
import numpy as np
class for_w:
"""
wの極値探索のためのクラス型
"""
def __init__(self, dist, n,
n_pers, n_dim, x):
self.dist = dist
self.n = n
self.n_pers = n_pers
self.n_dim = n_dim
self.x = x
# xは固定
def f(self, z):
"""
１次元配列（ベクトル）zを2次元配列wに並べ換える
"""
w =
np.resize(z, (self.n_pers, self.n_dim))
w =
w**2
d_hat = []
for
s in range(self.n_pers):
d_temp = np.zeros((self.n,self.n))
v = self.x * w[s]
for i in range(self.n):
for j in range(self.n):
d_temp[i][j] = np.max(np.abs(v[i] - v[j]))
d_hat.append(d_temp)
d_hat = np.array(d_hat)
err
= np.sum((self.dist - d_hat)**2)
return err
class for_x:
"""
xの極値探索のためのクラス型
"""
def __init__(self, dist, n,
n_pers, n_dim, w):
self.dist = dist
self.n = n
self.n_pers = n_pers
self.n_dim = n_dim
self.w = w
# wは固定
def f(self, z):
"""
１次元配列（ベクトル）zを２次元配列xに並べ換える
"""
x =
np.resize(z, (self.n, self.n_dim))
d_hat = []
for
s in range(self.n_pers):
d_temp =
np.zeros((self.n,self.n))
v = x * self.w[s]
for i in range(self.n):
for j in range(self.n):
d_temp[i][j] = np.max(np.abs(v[i] - v[j]))
d_hat.append(d_temp)
d_hat = np.array(d_hat)
err
= np.sum((self.dist - d_hat)**2)
return err