2010-04-10
多次元尺度法で遊んでみる(オレ流 R入門)
多次元データをクラスタリングする際に、それらのデータを2次元データに落とし込んで可視化させたいことがあります。そんな時に便利なのが「多次元尺度法」という手法です。
個々のデータ間の距離/類似度が分かっている場合に、それらのデータの座標を求めて、データ構造を復元するようなものです。
詳しい説明は割愛します。知りたい人はwikipediaと金先生の連載を読んで下さい。
体で覚えるタイプなので、とにかく何かデータを処理してみます。
「山手線」の地図を再現
さっそく試してみます。
山手線の各駅同士の直線距離を測っておいて、そのデータから実際の位置関係を復元できるか実験してみます。
山手線全駅の距離を測るのはめんどいので、適当に抜粋してしらべました。
以下のような表になりました。単位はメートルです。
さてさて、この距離表からどのようなデータ構造が再現されるでしょうか?
このデータを統計解析ソフトRで読み込めるようなデータにします。ファイル名はyamanote.txtとしておきます。
0 3609 7880 9026 9217 6385 4742 4092
3609 0 4328 5911 7410 6046 7679 7477
7880 4328 0 2666 6275 7305 11326 11469
9026 5911 2666 0 4069 6284 11504 12003
9217 7410 6275 4069 0 3744 9886 10871
6385 6046 7305 6284 3744 0 6157 7230
4742 7679 11326 11504 9886 6157 0 1592
4092 7477 11469 12003 10871 7230 1592 0
一行目がラベルで、2行目以降はデータが並んでいます。区切り文字はスペースだと思ってみて下さい。
それではこのファイルをRで読み込んで多次元尺度法で処理してみます。
多次元尺度を求めるにはcmdscaleという関数を叩きます。
X <- read.table("c:\\yamanote.txt", header=T) loc <- cmdscale(X)
続いてこれを描画してみます。
plot(loc) text(loc, names(X), col="red")
そうするとこのようなグラフが描画されます。
ん??角度と向きが変ですね。。。
というのも、多次元尺度法ではあくまでも距離関係だけでデータ構造を再現するので、
今回のような実際の地理情報を扱う場合にはちょっと手直ししてあげる必要があります。
手直しといっても簡単です。裏表が逆なのでひっくり返して、時計回りに90度動かしてみます。
ぬぬ!これは・・・東京在住の人はすでに興奮して鳥肌ものですね!?
関東以外の方は山手線の駅名に馴染みがないかもしれないので、参考までに実際の山手線の駅画像を貼付けてみましょう。
この駅画像の上に、先ほど描画したファイルをちょっと角度をずらして、幅を調整しながら重ねてみると・・・。
おおおおおお!ぴったりだぁぁぁ!
感動したぁ!!
「ガンダム」でモビルスーツを分析
お次は実際の距離ではなく、適当な特徴ベクトルを作ってマップしてみようと思います。
世の中にはガンダムおたく、通称ガンヲタの方が多いようで、インターネット上にはモビルスーツのスペック表を公開しているサイトもあったりします。
http://www.interq.or.jp/jupiter/mcmurd/FG.htm
今回はこのスペック表から特徴ベクトルを作って、個々のモビルスーツ/モビルアーマーの類似度を測定してみます。それでもって多次元尺度法で2次元に「見える化」してみようという実験です。
なお素性として使うデータは「頭頂高」「本体重量」「ジェネレータ出力」「スラスター推力」の4項目とします。
データを作る部分はperlでかきました。
たぶんRだけでも出来るんだと思いますが、Rは初心者なので、細かい処理をどう書けばいいのかわからないのでperlに逃げました。
use strict; use warnings; #-- データ読み込み my $vectors; my @names; while(<DATA>){ chomp $_; my @f = split("\t", $_); my $id = $f[0]; my $name = $f[1]; my $height = $f[2]; $height =~ s/m//g; my $weight = $f[3]; $weight =~ s/t//g; my $kwat = $f[4]; $kwat =~ s/kw//g; my $gain = $f[5]; $gain =~ s/kg//g; my $weapon = $f[6]; $vectors->{$name} = { height => $height, weight => $weight, kwat => $kwat, gain => $gain, }; } #-- 単位ベクトルに正規化 while ( my ( $label, $vec ) = each %$vectors ) { unit_length($vec); $vectors->{$label} = $vec; } #-- 個々の距離(1−類似度)を量ってマトリックス化 my $matrix; my @labels; for my $label (keys %$vectors){ my $vec = $vectors->{$label}; for my $target_label( keys %$vectors){ my $target_vec = $vectors->{$target_label}; my $sim = cosine_similarity($vec, $target_vec); my $dis = 1 - $sim; $matrix->{$label}->{$target_label} = $dis; } push(@labels, $label); } #-- 出力 print join(" ", @labels),"\n"; for my $row(@labels){ my @array; for my $col(@labels){ my $item = $matrix->{$row}->{$col}; push @array, $item; } print join(" ", @array), "\n"; } #---- 以下サブルーチン sub unit_length { my $vec = shift; my $norm = norm($vec); while ( my ( $key, $value ) = each %$vec ) { $vec->{$key} = $value / $norm; } } sub norm { my $vec = shift; my $norm; for ( values %$vec ) { $norm += $_**2; } sqrt($norm); } sub cosine_similarity { my ( $vector_1, $vector_2 ) = @_; my $inner_product = 0.0; map { if ( $vector_2->{$_} ) { $inner_product += $vector_1->{$_} * $vector_2->{$_}; } } keys %{$vector_1}; my $norm_1 = 0.0; map { $norm_1 += $_**2 } values %{$vector_1}; $norm_1 = sqrt($norm_1); my $norm_2 = 0.0; map { $norm_2 += $_**2 } values %{$vector_2}; $norm_2 = sqrt($norm_2); return ( $norm_1 && $norm_2 ) ? $inner_product / ( $norm_1 * $norm_2 ) : 0.0; } __DATA__ MS-05 ザクI 17.5m 50.3t 899kw 40700kg 105mm マシンガン、240mmバズーカ、ヒートホーク MS-06F ザクII 17.5m 58.1t 976kw 43300kg 120mm マシンガン、240mmバズーカ、ヒートホーク、他 MS-06S ザクII 17.5m 56.2t 976kw 51600kg 120mm マシンガン、240mmバズーカ、ヒートホーク、他 MS-07B グフ 18.2m 58.5t 1034kw 40700kg 5連装75mmマシンガン、ヒートロッド、ヒートサーベル MS-09 ドム 18.6m 62.6t 1269kw 58200kg 360mm バズーカ、拡散ビーム砲、ヒートサーベル MS-09R リックドム 18.6m 43.8t 1199kw 53000kg ジャイアントバズーカ、拡散ビーム砲、ヒートサーベル MS-14A ゲルググ 19.2m 42.1t 1440kw 61500kg ビームライフル、ビームナギナタ MSM-03 ゴッグ 18.3m 82.4t 1740kw 121000kg メガ粒子砲×2、魚雷発射管×2 MSM-04 アッガイ 19.2m 91.6t 1870kw 109600kg 105mm バルカン×4、ロケット弾ランチャー×6 MSM-07 ズゴック 18.4m 65.1t 2480kw 83000kg 240mm ロケット弾×8、メガ粒子砲×8 MSM-10 ゾック 23.9m 167.6t 3849kw 253000kg フォノンメーザー砲、メガ粒子砲×8 MA-05 ビグロ 45.5m 125.5t 17800kw 136100kg 大型メガ粒子砲、4連装ミサイルランチャー×2 MA-08 ビグ・ザム 59.6m 1021.2t 140000kw 580000kg 大型メガ粒子砲、対空メガ粒子砲×28、105mmバルカン×2 MAN-02 ジオング 17.3m 151.2t 9400kw 187000kg 有線制御式メガ粒子砲×2、メガ粒子砲×3 MAN-03 ブラウ・ブロ 60.2m 1735.3t 74000kw 1760000kg 有線制御式メガ粒子砲×4 MAN-08 エルメス 85.4m 163.7t 14200kw 645200kg メガ粒子砲×2、ビット×10 RX-75 ガンタンク 15.0m 56.0t 878kw 88000kg 120mm 低反動キャノン砲×2、4連装ガンランチャー×2 RX-77-2 ガンキャノン 17.5m 51.0t 1380kw 51800kg 240mm キャノン、ビームライフル、60mmバルカン砲×2、他 RX-78-2 ガンダム 18.0m 43.4t 1380kw 55500kg ビームライフル、ビームサーベル、60mmバルカン砲×2、ハイパーバズーカ、他 RGM-79 ジム 18.0m 41.2t 1250kw 55500kg ビームスプレーガン、ビームサーベル、60mmバルカン砲×2、ハイパーバズーカ
こうして出力した結果をテキストファイルに保存してRで処理します。
X <- read.table("c:\\gandum.txt", header=T) loc <- cmdscale(X) plot(loc, type="n") text(loc, names(X), col="red")
結果はこう。
扱ったデータが重量とか大きさ、馬力のような類いのものなので、そのことを念頭に入れてみてみると、なんとなくこの距離感はうなずけるかな、と思います。
ビグザム、ビグロ、エルメス、ジオング、ブラウブロあたりの重量級が程よく散らばってますね。
またさりげなくガンタンクも個性的なポジショニングなのが笑えます。
本当は個々の素性を標準化しないといけないような気がしますが、まぁ雰囲気でてるんで良しとします。
それにしても、これだとなんだかイマイチなんで、画像をプロットしてみましょう。
おお!なんだか猛烈にガンプラ作りたくなってきた!
というわけで、Rと多次元尺度法、結構あそべますyo。
- 83 http://www.google.co.jp/search?sourceid=navclient&hl=ja&ie=UTF-8&rlz=1T4GPEA_jaJP294JP295&q=ベイズ
- 39 http://www.google.co.jp/search?hl=ja&source=hp&q=ナップサック問題&lr=&aq=0r&aqi=g-r10&aql=&oq=nappusakku&gs_rfai=
- 37 http://reader.livedoor.com/reader/
- 34 http://d.hatena.ne.jp/
- 34 http://www.google.co.jp/search?q=ナップサック問題&hl=ja&lr=&sa=2
- 33 http://www.google.co.jp/search?hl=ja&client=firefox-a&rls=org.mozilla:ja:official&q=df+idf+ツール&lr=&aq=f&aqi=&aql=&oq=&gs_rfai=
- 33 http://www.google.co.jp/search?hl=ja&q=ナップサック問題&sourceid=navclient-ff&rlz=1B3GGGL_jaJP323JP323&ie=UTF-8
- 32 http://www.google.co.jp/search?hl=ja&safe=off&client=firefox-a&hs=CKq&rls=org.mozilla:ja:official&q=合わせ鏡+実験&lr=&aq=3&aqi=g10&aql=&oq=合わせ鏡&gs_rfai=
- 30 http://www.google.co.jp/search?q=Initialized+empty+Git+repository+in&ie=utf-8&oe=utf-8&aq=t&rls=com.ubuntu:ja:unofficial&client=firefox-a
- 29 http://www.google.co.jp/search?q=動的計画法&hl=ja&start=10&sa=N