Up

個人差尺度法（INDSCAL: Individual Differences Scaling）

重み付き距離に対する直接交互最小2乗法

INDSCAL (Carroll & Wish, 1974)では、個人の空間における事物との距離は次式で与えられる。

Carroll & Wish (1974)では、式（1）をCANDECOMP (CANonical DECOMPosition)の形

に変形して解の求め方が説明されている。CANDECOMPは、PARAFACとかCP (Canonical Polyadic decomposition)とも呼ばれている。式（２）の解法として、パラメータ群ごとに交互に最小2乗解を求める（ALS: Alternative Least Squaraes）交互最小２乗法が説明されている。

多次元尺度構成法において、直接距離との差を極値探索法により最小にする布置を求める方法が試みられている（計量モデル、非計量モデル）が、直接、式（１）のモデルに対して交互にウェイトと布置の値を極値探索法により求めてみた。距離（１）を直接対象とする方法であるので、ユークリッド距離以外のr-metricの場合に簡単に拡張できる。また、極値探索法は、Powell法を用いたので、偏導関数などは不要である。

以下のようなアルゴリズムになる。

与えられたデータをと置く。

（１）を固定して、式（１）の値の2乗との2乗との差の２乗和が最小になるを極値探索法で求める。

（２）を固定して、式（１）の値の2乗との2乗との差の２乗和が最小になるを極値探索法で求める。

は、中心が原点、各座標の絶対値の最大値が１となるようにスケーリングを行う。のスケーリングに合わせてのスケーリングを行う。

（３）値が収束しているか調べる。収束していないときは、ステップ（１）に戻る。

上のアルゴリズムによるINDSCALのスクリプトをリスト１のように用意した。ファイルは、圧縮ファイルindscalmetricfiles.zipにまとめた。

リスト１では、極値探索法の関数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ファイルとして用意する。例rawdist_rnd.csvを図１に示す。

図１　rawdist_rnd.csv

図１のデータは、事物A、B、C、D、E、F、G、H、Iの距離を、人P1、P2、・・・について表形式で設定したものである。

１行目の第２列から事物の名前を入力する。

第２行目には、１列目に人の名前を設定する。

第３行目から、先頭に事物の名前を記入した後、その事物と各事物との距離を設定していく。

すべての事物についてのデータを設定後、次の行の１列目に次の人の名前を設定する。名前の設定行の次の行から事物の距離を設定していく。

上の手順で、すべての人の事物についての距離のデータを設定する。

設定したデータの分析は、リスト１のスクリプトの実行時に入力データファイルとして指定する。

リスト１のスクリプトを実行すると、入力データファイル名の設定が求められる。次の例は、図１のファイル名を設定したものである。

Input File Name (*.csv) = rawdist_rnd.csv

入力データファイル名を設定すると、ファイルが読み込まれて分析が始まる。

まず、の初期値が設定される。これは、距離データの平均値に対して２重中心化により得られる座標値を初期値としている。の初期値はすべて１としている。

２重中心化された配列の固有分解の固有値が、正の値のみスクリープロットの形式で表示される（図２）。

図２

このスクリープロットから、布置の次元数の検討を付ける。図２のグラフの場合、次元数は２が適当であると思われる。

図２のWindowを閉じると、次元数を聞いてくる。次の例は、次元数を２と設定したものである。

Number of dimensions(<=5) = 2

次元数を設定すると極値探索法の計算が始まる。解が収束すると事物の布置のグラフが表示される（図３）。

図３

図３のWindowを閉じると、重みのグラフが表示される（図４）。

図４

図４のWindowを閉じると、次に描画する次元（座標軸）の設定が求められる。負の整数を設定すると、描画終了である。

次の例では、描画を終了するためにー１を設定している。

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.10153

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.

リスト１　計量INDSCAL

import numpy as np

import matplotlib.pyplot as plt

import pandas as pd

import scipy.optimize as so

from myopt import * # 極値探索用関数のクラス型

fin_nm = input('Input File Name (*.csv) = ')

with open(fin_nm, 'r') as f:

raw_data = pd.read_csv(f)

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)

d_ijk = []

for i in range(n_ind):

d_ijk.append(v_data[1 + (n + 1)*i:(n + 1)*(i+1)])

d_ijk = np.array(d_ijk)

print('Distances =\n', d_ijk)

#

# 初期値設定用の平均距離の配列

#

m_dist = np.mean(d_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 # 初期値

i_trial = 0

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

#

# 解の布置の描画

#

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}')

リスト２　極値探索のためのクラス型

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

Home