Rand index によるクラスタリング間の類似度の算出
データのクラスタリング間の類似度を測る尺度として、ランド指数 (Rand Index) というものがあるようです。プログラムを書きながら試してみたので、実例とともに計算方法をまとめてみたいと思います*1。なお Wikipedia の英語ページは以下にあります。日本語ページは未作成のようです。
Rand index - Wikipedia, the free encyclopedia
クラスタリングの対象になるような手ごろなデータがないかと探してみたところ、国土地理院のウェブページで都道府県庁所在地の緯度経度座標が公開されているのを見つけました。今回はこのデータを使ってみることにします。
都道府県庁の経度緯度
まずは HTML ファイルを取得して、この後の処理で使いやすいように csv 形式に変換しておきます。元データは緯度経度が「度分秒」の形式になっていますが、これも十進法に変換しておきます。シンプルな HTML だったので、以下のようなコードで一連の処理を実行できました。なお、awk に渡す直前の sed は全角空白を半角空白に置換しているものです。
$ wget http://www.gsi.go.jp/KOKUJYOHO/kencho/kenchobl.html $ sed 's/\r//g' kenchobl.html | grep '<TD>' | tail -n +4 | sed 's/<[^>]*>//g' | sed 's/ / /g' | awk ' (NR % 3 == 1) { gsub(" ", ""); printf("%s,", $0); } (NR % 3 == 2) { printf("%f," , $1 + $2 / 60 + $3 / 3600); } (NR % 3 == 0) { printf("%f\n", $1 + $2 / 60 + $3 / 3600); } ' >kenchobl.csv
出来上がりはこんなファイルになります。
$ head -n 5 kenchobl.csv 北海道,141.346944,43.064167 青森,140.740000,40.824444 岩手,141.152500,39.703611 宮城,140.871944,38.268889 秋田,140.102500,39.718611
これを k-means 法でクラスタリングします。今回は R の kmeans 関数を利用してみました。R は使い慣れていないのですが、どうやら下記のように実行できるようです。クラスタリング間の類似度を算出したいので、何度も繰り返して複数回の結果を出しておきます。ここでは 100 回実行しました。また、クラスタ数は 5 としました。このあたりの設定には特に根拠はありません。
data <- read.csv('kenchobl.csv', header=F) results = matrix(nrow=nrow(data), ncol=100) for (i in 1:ncol(results)) { results[,i] <- matrix(kmeans(data.frame(lng=data[,2], lat=data[,3]), 5)$cluster, nrow=nrow(results), ncol=1) } write.table(results, 'results.csv', sep=',', row.names=F, col.names=F)
出力される results.csv はこのようなファイルになります。全体は 47 行 100 列の行列で、各行が都道府県、各列が 1 回ごとの k-means の結果になり、行列の値がクラスタ番号を表します。
$ head -n 5 results.csv | cut -f1-20 -d, 3,1,4,1,2,1,2,3,4,5,4,5,2,1,5,3,1,1,1,3 3,1,4,1,2,1,2,3,4,5,4,5,2,1,5,3,1,1,1,3 1,1,3,1,2,1,2,3,4,5,4,5,2,3,5,3,1,1,1,3 1,1,3,1,2,1,2,3,2,5,4,5,2,3,5,3,1,1,1,3 1,1,3,1,2,1,2,3,4,5,4,5,2,3,5,3,1,1,1,3 $ awk -F, '{ print NF }' results.csv | uniq -c 47 100
確認のため、いくつか結果を見てみます。1 回目の k-means の結果はこのようになっていました。近い都道府県が同一クラスタに割り当てられています。
$ paste -d, results.csv kenchobl.csv | cut -f1,101 -d, | awk -F, '{ prefs[$1] = prefs[$1] " " $2 } END { for (i in prefs) print i prefs[i] }' | sort 1 岩手 宮城 秋田 山形 福島 新潟 2 広島 山口 愛媛 福岡 佐賀 長崎 熊本 大分 宮崎 鹿児島 沖縄 3 北海道 青森 4 石川 福井 岐阜 愛知 三重 滋賀 京都 大阪 兵庫 奈良 和歌山 鳥取 島根 岡山 徳島 香川 高知 5 茨城 栃木 群馬 埼玉 千葉 東京 神奈川 富山 山梨 長野 静岡
それに対して 2 回目の k-means の結果はこのようになっており、実行するたびに結果が異なっていることがわかります。
$ paste -d, results.csv kenchobl.csv | cut -f2,101 -d, | awk -F, '{ prefs[$1] = prefs[$1] " " $2 } END { for (i in prefs) print i prefs[i] }' | sort 1 北海道 青森 岩手 宮城 秋田 山形 2 鳥取 島根 岡山 広島 山口 徳島 香川 愛媛 高知 3 福島 茨城 栃木 群馬 埼玉 千葉 東京 神奈川 新潟 山梨 長野 静岡 4 富山 石川 福井 岐阜 愛知 三重 滋賀 京都 大阪 兵庫 奈良 和歌山 5 福岡 佐賀 長崎 熊本 大分 宮崎 鹿児島 沖縄
さて、ようやく本題に入ります。ランド指数は、このような二つの分け方の間の類似度を表す尺度です。すべてのデータ対のうち、同一クラスタに割り当てられたかどうかが一致するデータ対の比率として表されます。上記の 2 回の k-means の結果を使って考えてみます。下図のように行と列に 47 都道府県を並べ、すべての都道府県ペアに対して、k-means で同一クラスタに割り当てられたかどうかを○×で記入します。
たとえば、北海道と青森は、1 回目の k-means でも 2 回目の k-means でも同じクラスタに割り当てられましたので、セルの値を「○ ○」とします。一方、北海道と岩手は、1 回目の k-means では異なるクラスタに割り当てられ、2 回目の k-means では同じクラスタに割り当てられましたので、「× ○」とします。北海道と茨城の場合は、1 回目も 2 回目も異なるクラスタに割り当てられましたので「× ×」とします。このようにして各セルに○×を埋めた後、すべてのセルのうち「○ ○」と「× ×」の占める割合を計算したものがランド指数になります。
これは簡単な計算なので PHP の関数として書いてみました。安直に実装すると、このようなコードになります (これは効率の悪い実装なので、すぐ後でより効率的な実装に書き直します)。
<?php function calculateRI($c1, $c2) { $n = count($c1); $numer = 0; $denom = 0; for ($i = 0; $i < $n; ++$i) { for ($j = $i + 1; $j < $n; ++$j) { if (($c1[$i] == $c1[$j]) == ($c2[$i] == $c2[$j])) { ++$numer; } ++$denum; } } return $numer / $denom; }
続いて、ファイルからデータを読み込んでランド指数を計算するコードを記述します。
<?php $c1 = array(); $c2 = array(); while ($line = fgets(STDIN)) { list($i, $j) = explode(',', trim($line)); $c1[] = $i; $c2[] = $j; } print calculateRI($c1, $c2) . "\n";
このプログラムを実行すると、このようにランド指数を計算できました。
$ cut -d, -f1,2 results.csv | php calculateRI.php 0.84736355226642
ここで実装した calculateRI 関数の計算量はデータ数を N としたとき O(N^2) になりますが、少し工夫すると計算量を削減できます。まず、2 回の k-means の結果を分割表 (contingency table) の形にまとめます。行方向を 1 回目の k-means のクラスタ番号、列方向を 2 回目の k-means のクラスタ番号として、各セルには該当するデータの個数を記入します。たとえば、1 行 1 列のセルには、岩手、宮城、秋田、山形の 4 県が該当するので、セルの値は 4 となります。すべてのセルに値を埋めると以下のようになります。
次に、この分割表の各セルから二つの都道府県を選ぶときの組み合わせ数を計算します。セルの値を n とすれば n * (n - 1) / 2 通りになります。これをまとめたものが次の表です。
この表を用いると、最初の表での「○ ○」「× ×」の個数を簡単に計算できます。○×の組み合わせを次の表のようにまとめると、色を付けた各セルの値は上表での同色セルの合計値に一致します。その他のセルの値は引き算で計算できます。各セルの値からランド指数を計算すると RI = (154 + 762) / 1081 = 0.84736... となり、先ほどのプログラムで計算した値と一致します。
これを PHP のコードで書くと、次のようになります。
<?php function calculateRI($c1, $c2) { $n = count($c1); $k1 = array_unique($c1); $k2 = array_unique($c2); $cell = array_fill_keys($k1, array_fill_keys($k2, 0)); for ($i = 0; $i < $n; ++$i) { ++$cell[$c1[$i]][$c2[$i]]; } $total = comb2($n); $oo = array_sum(array_map(function ($r) { return array_sum(array_map('comb2', $r)); }, $cell)); $ox = array_sum(array_map('comb2', array_count_values($c1))) - $oo; $xo = array_sum(array_map('comb2', array_count_values($c2))) - $oo; $xx = $total - $oo - $ox - $xo; return ($oo + $xx) / $total; } function comb2($n) { return $n * ($n - 1) / 2; }
さて、ここまで説明してきたランド指数ですが、一つ問題があります。以下のようにして、ランダムなクラスタリングでランド指数を計算してみます。
$ yes | head -n 47 | awk 'BEGIN { srand(); } { printf("%d,%d\n", rand() * 5, rand() * 5); }' | php calculateRI.php 0.6845513413506
実行結果が示すように、ランダムなクラスタリング同士でも 0.68 という値が出てしまいました。ランド指数は、二つのクラスタリングの間に相関がないときの値が一定にならないので、そのまま扱うには不便なことがあります。調整ランド指数 (Adjusted Rand Index) は、相関がない場合に値が 0 となるように計算式を調整したものです。
調整ランド指数では、2 回のクラスタリングが互いに独立だとした場合の「○ ○」「× ×」の期待値を計算します。「○ ○」について考えると、2 回のクラスタリングが独立なら期待値は下式のように計算できます。同様に「× ×」について計算すると 659.14 と求まり、「○ ○」と「× ×」の合計の期待値は 51.14 + 659.14 = 710.28 になります。
調整ランド指数は、ランド指数の分子と分母のそれぞれから、この期待値を引いたものになります。今回の値で計算すると ARI = (916 - 710.28) / (1081 - 710.28) = 0.554921 になりました。
最後に、調整ランド指数を計算する関数を PHP で書いてみます。コメントアウトしている部分は途中の式になります。途中ではたくさんの項が出現しますが、あちこちキャンセルして整理すると最終的に return 文に書かれた式になり、これが Wikipedia の ARI の計算式に一致します。
<?php function calculateARI($c1, $c2) { $n = count($c1); $k1 = array_unique($c1); $k2 = array_unique($c2); $cell = array_fill_keys($k1, array_fill_keys($k2, 0)); for ($i = 0; $i < $n; ++$i) { ++$cell[$c1[$i]][$c2[$i]]; } $oo = array_sum(array_map(function ($r) { return array_sum(array_map('comb2', $r)); }, $cell)); $c1o = array_sum(array_map('comb2', array_count_values($c1))); $c2o = array_sum(array_map('comb2', array_count_values($c2))); $total = comb2($n); // $ox = $c1o - $oo; // $xo = $c2o - $oo; // $xx = $total - $oo - $ox - $xo // = $total + $oo - $c1o - $c2o; // // $expected = ($c1o * $c2o + ($total - $c1o) * ($total - $c2o)) / $total // = 2 * $c1o * $c2o / $total - $c1o - $c2o + $total; // // $numer = $oo + $xx - $expected // = 2 * $oo - 2 * $c1o * $c2o / $total; // $denom = $total - $expected // = $c1o + $c2o - 2 * $c1o * $c2o / $total; return ($oo - $c1o * $c2o / $total) / (($c1o + $c2o) / 2 - $c1o * $c2o / $total); }
[2014-02-08] 追記: Git の勉強の一環として、https://github.com/y-uti/ari にプログラムを公開しました。
