Loading web-font TeX/Math/Italic

January 5, 2012

Pythonで多次元尺度構成法を実装してみる

多次元尺度構成法(Multidimensional scaling; MDS)は多変量解析の手法です.よくデータ間の(非)類似度の情報可視化に用いられます.MDSは基本的には似ているものは近くに,似ていないものは遠くに配置するような座標を求めます.ここでは古典的(計量的; metric)MDSをPythonで実装してみます.

古典的MDSは以下の手順で可視化.実装に必要なところのみ.

  • 要素の値が距離の2乗と見なせる非類似度行列Sを用意する
  • Sにヤング・ハウスホルダー変換を施してPとする
  • Pをスペクトル分解する(固有値・固有ベクトルを求める)
  • 固有値の大きい方から2~3個選び,対応する固有ベクトルを取り出す
  • 各固有ベクトルの要素値をプロットする(2個の固有ベクトルの時は2次元,3個の固有ベクトルの時は3次元)

一応,下部に参考リンクを挙げましたが,詳しい説明は検索すれば山ほど….ヤング・ハウスホルダー変換は距離の2乗の行列に両側から中心化行列をかける演算のことらしいです.中心化行列は距離の基準の原点を重心に移動するための行列らしいです.

ヤング・ハウスホルダー変換:
SがN×N行列のとき
P = -\frac{1}{2}({\bf I}-\frac{1}{n}{\bf 11}')S({\bf I}-\frac{1}{n}{\bf 11}')
と計算します.1は要素がすべて1のN次ベクトル[1,...,1]'です.なので11'は要素がすべて1のN×N行列.

アメリカの空港間の距離が与えられているときMDSで可視化してみます.
atl   chi   den   hou    la    mi    ny    sf   sea    dc
atl     0
chi   587     0
den  1212   920     0
hou   701   940   879     0
 la  1936  1745   831  1374     0
 mi   604  1188  1726   968  2339     0
 ny   748   713  1631  1420  2451  1092     0
 sf  2139  1858   949  1645   347  2594  2571     0
sea  2182  1737  1021  1891   959  2734  2408   678     0
 dc   543   597  1494  1220  2300   923   205  2442  2329     0

以下コード.要numpy, matplotlib.

#!/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()
view raw mds.py hosted with ❤ by GitHub


結果はこんなかんじ.ラベルは手作業で付けました.


ちなみに実際の位置はこんなかんじ.

View US 10 Airport in a larger map

結果を180度回転すればだいたいなんとなくあってるかな?

上の例では非類似度行列を実際の距離の行列としましたが,MDSはデータ間の非類似度を測ることができたらそれっぽく2次元(3次元)にマッピングできるよ,ってことなので色々使えるかも知れません.

ちなみにRだったらcmdscaleでMDSが実装されてて楽

No comments:

Post a Comment