Multidimensional Scaling

mds.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
多次元尺度構成法
Multidimensional Scaling; MDS
"""

import numpy as np
import matplotlib.pyplot as plt

# airline distances between 10 US cities
d = """0,587,1212,701,1936,604,748,2139,2182,543;
587,0,920,940,1745,1188,713,1858,1737,597;
1212,920,0,879,831,1726,1631,949,1021,1494;
701,940,879,0,1374,968,1420,1645,1891,1220;
1936,1745,831,1374,0,2339,2451,347,959,2300;
604,1188,1726,968,2339,0,1092,2594,2734,923;
748,713,1631,1420,2451,1092,0,2571,2408,205;
2139,1858,949,1645,347,2594,2571,0,678,2442;
2182,1737,1021,1891,959,2734,2408,678,0,2329;
543,597,1494,1220,2300,923,205,2442,2329,0"""

# データ読み込み
D = np.array(np.mat(d))
N = len(D)

# 距離の2乗の行列 (arrayだと要素同士の掛け算になる)
S = D * D

# 中心化行列
one = np.eye(N) - np.ones((N,N))/N

# ヤング・ハウスホルダー変換
P = - 1.0/2 * one * S * one

# スペクトル分解
w,v = np.linalg.eig(P)
ind = np.argsort(w)
x1 = ind[-1] # 1番
x2 = ind[-2] # 2番
print w[x1],w[x2]

# 標準されたデータの固有値が求められているので標準偏差を掛けて可視化
s = P.std(axis=0)
w1 = s[x1]
w2 = s[x2]
for i in range(N):
    plt.plot(w1*v[i,x1],w2*v[i,x2],'b.')
plt.draw()
plt.show()