y_uti のブログ

統計、機械学習、自然言語処理などに興味を持つエンジニアの技術ブログです

k-medoids 法と DTW による時系列データのクラスタリング

過去の台風の経路情報を題材として、k-medoids 法による時系列データのクラスタリングを試してみました。距離の尺度には、以前の記事*1でも試した Dynamic Time Warping (DTW) を利用しました。K-medoids 法と DTW については、それぞれ Wikipedia に説明があります。
k-medoids - Wikipedia, the free encyclopedia
Dynamic time warping - Wikipedia, the free encyclopedia

結果は次のようになりました。気象庁の台風位置表 (http://www.data.jma.go.jp/fcd/yoho/typhoon/position_table/) から過去の台風の経路情報を取得して、k-medoids 法で 5 クラスタに分類した結果です。今回の記事では、この図を描画するまでの手順を説明します。

データの準備

過去の台風の資料は気象庁のウェブページで公開されています*2。「台風位置表」のページから、2001 年以降の台風の資料を PDF ファイルでダウンロードできます。今回は、これらの PDF ファイルから時刻ごとの緯度経度を抽出して、時系列データとして利用します。
気象庁｜過去の台風資料

上記のウェブページから 2001 年以降のすべての台風の PDF ファイルをダウンロードして、それぞれテキストファイルに変換します。テキストファイルへの変換には pdftotext を利用しました。2001 年の台風 1 号での実行例は次のとおりです。pdftotext の処理結果が T0101.txt に保存されます。

$ wget http://www.data.jma.go.jp/fcd/yoho/data/typhoon/T0101.pdf
$ pdftotext T0101.pdf

次に、このテキストファイルから緯度と経度を抽出して csv ファイルに変換します。以前の記事にも書きましたが、このような処理はデータに合わせてアドホックな工夫が必要になります。今回は以下のようなスクリプトを作成して処理しました。

#!/bin/bash
tr ' -' '\n'                     |\
grep '\(\.\|^[NEWS]$\)'          |\
sed 's/\(.\)\([NEWS]\)$/\1\n\2/' |\
awk '
    BEGIN {
        nextval = "lat";
    }
    $1 ~ /[NS]/ { latdir = $1; next; }
    $1 ~ /[EW]/ { lngdir = $1; next; }
    nextval == "lat" && NR != 1 {
        writeline(lat, latdir, lng, lngdir);
    }
    nextval == "lat" { lat = $1 % 100; nextval = "lng"; next; }
    nextval == "lng" { lng = $1      ; nextval = "lat"; next; }
    END {
        writeline(lat, latdir, lng, lngdir);
    }

    function writeline(lat, latdir, lng, lngdir) {
        printf("%.1f,%.1f\n", latdir == "N" ? lat : -lat, lngdir == "E" ? lng : 360 - lng);
    }
'

スクリプトの前半では、まず tr でデータを値ごとに改行した後、grep で緯度経度の数値と、北緯、東経等を表す文字を含む行を抽出します。データの中には、緯度経度の数値の後に空白無しで北緯、東経などの文字が続く場合もあったので、そのパターンを sed で拾って整えます。ここまでの処理で次の状態になります。元データの pdf ファイル*3と内容を比較してみると、緯度経度をうまく抽出できていることがわかります。

$ cat T0101.txt | tr ' -' '\n' | grep '\(\.\|^[NEWS]$\)' | sed 's/\(.\)\([NEWS]\)/\1\n\2/' | head
11.8
N
120.0
E
12.3
119.6
13.1
119.5
14.1
119.3

後半の awk スクリプトは、ここから緯度と経度を一行にまとめて csv 形式で出力します*4*5。南緯、西経の場合にはそれぞれ北緯、東経の値に変換しています。

さて、このスクリプトを extract_latlng.sh として、次のように実行します。

$ extract_latlng.sh <T0101.txt >T0101.csv

処理結果は以下のとおりです。各行が時刻ごとの台風の位置を表しており、1 列目が北緯、2 列目が東経の値になります。

$ head T0101.csv
11.8,120.0
12.3,119.6
13.1,119.5
14.1,119.3
15.3,119.1
16.2,119.1
17.1,119.1
17.6,119.1
18.0,119.2
18.4,119.4
距離行列の計算

ここまでの処理で、2001 年から 2015 年までの 355 個の台風の経路情報を csv 形式で抽出できました。ここで、これらのすべての組み合わせについて DTW による距離を計算し、355 行 355 列の距離行列を作成します。

時系列データの距離を計算するプログラムは以下のとおりです。これは、以前の記事で作成したプログラムをそのまま利用しています。

<?php
function dtw($a, $b, $distance = 'euclid') {
    $d = array_fill(0, count($a) + 1, array_fill(0, count($b) + 1, INF));
    $d[0][0] = 0;
    for ($i = 1; $i <= count($a); ++$i) {
        for ($j = 1; $j <= count($b); ++$j) {
            $d[$i][$j] =
                min([$d[$i - 1][$j - 1], $d[$i][$j - 1], $d[$i - 1][$j]])
                + $distance($a[$i - 1], $b[$j - 1]);
        }
    }
    return $d[count($a)][count($b)];
}

function euclid($a, $b) {
    return sqrt(array_sum(array_map(function ($x, $y) {
                    return pow($x - $y, 2);
                }, $a, $b)));
}

function readCsv($filename) {
    return array_map(
        function ($line) { return explode(',', $line); },
        file($filename, FILE_IGNORE_NEW_LINES));
}

echo dtw(readCsv($argv[1]), readCsv($argv[2])) . "\n";

このプログラムを dtw.php として、すべての組み合わせについて DTW 距離を計算します。次のように使い捨てのスクリプトを作成して実行しました。結果は 355 行 355 列の対称行列になります。なお、手元の環境ではこの行列の計算に数時間かかりました。結果が対称行列になることはわかっているので、実装を工夫すれば計算時間は半分で済むはずですが、このスクリプトでは何も工夫せずに行列の全ての要素を計算しています。

$ cat build_distance_matrix.sh
#!/bin/bash
for a in T*.csv; do
    for b in T*.csv; do
        echo -n $(php dtw.php $a $b),
    done
    echo
done

$ ./build_distance_matrix.sh | sed -s ',$//' >dtw_distance_matrix.csv

得られた行列の左上 5 行 5 列を表示させてみます。対称行列になっていること、対角成分の値が 0 になっていることが確かめられます。

$ head -n 5 dtw_distance_matrix.csv | cut -d, -f1-5
0,212.8224774021,359.27459542857,354.74975039293,188.25687465072
212.8224774021,0,383.83645213076,201.19359206238,162.92601640849
359.27459542857,383.83645213076,0,270.84594012043,225.63593029795
354.74975039293,201.19359206238,270.84594012043,0,133.08475475109
188.25687465072,162.92601640849,225.63593029795,133.08475475109,0
K-medoids 法によるクラスタリング

前段の処理で得られた DTW の距離行列を利用して、今回の記事の主題である k-medoids 法によるクラスタリングを実行します。

K-medoids 法は、次の手順でデータを k 個のクラスタに分類する手法です*6。

  1. 与えられたデータから k 個を適当に選択する。これを k 個のクラスタそれぞれの medoid とする
  2. 収束するまで次の操作を繰り返す
    • 各データを、k 個のクラスタのなかで medoid との距離が一番近いものに割り当てる
    • 各クラスタについて、クラスタ内の各データとの距離の総和が最小になるデータを、そのクラスタの medoid とする

K-means 法ではクラスタ内の各データとの二乗距離の総和が最小になる点を計算して centroid としますが、k-medoids 法は、与えられたデータの中から medoid を選択するという点が k-means 法と異なります。K-medoids 法はデータ間の距離だけが求まればよいので、あらかじめ全データ間の距離行列を与えれば実行中に距離を計算する必要はありません。

K-medoids 法は以上のような単純なアルゴリズムで、一般的なプログラミング言語で簡単に実装できます。今回は PHP で次のように実装しました。

<?php
function kmedoids($distances, $k, $maxiter = 100) {
    $medoids = initialize_medoids($distances, $k);
    $indices = false;
    for ($iter = 0; $iter < $maxiter; ++$iter) {
        $next = assign_to_nearest($distances, $medoids);
        if ($next === $indices) {
            break;
        }
        $indices = $next;
        $medoids = update_medoids($distances, $indices, $k);
    }
    return array($indices, $medoids);
}

function initialize_medoids($distances, $k) {
    $medoids = range(0, count($distances) - 1);
    shuffle($medoids);
    return array_slice($medoids, 0, $k);
}

function assign_to_nearest($distances, $medoids) {
    $k = count($medoids);
    $indices = array();
    foreach ($distances as $d) {
        $mindist = INF;
        $nearest = 0;
        for ($i = 0; $i < $k; ++$i) {
            if ($d[$medoids[$i]] < $mindist) {
                $mindist = $d[$medoids[$i]];
                $nearest = $i;
            }
        }
        $indices[] = $nearest;
    }
    return $indices;
}

function update_medoids($distances, $indices, $k) {
    $n = count($distances);
    $mindists = array_fill(0, $k, INF);
    $medoids = array_fill(0, $k, false);
    for ($i = 0; $i < $n; ++$i) {
        $m = $indices[$i];
        $dist = 0;
        for ($j = 0; $j < $n; ++$j) {
            if ($indices[$j] == $m) {
                $dist += $distances[$i][$j];
            }
        }
        if ($dist < $mindists[$m]) {
            $mindists[$m] = $dist;
            $medoids[$m] = $i;
        }
    }
    return $medoids;
}

function readCsv($filename) {
    return array_map(
        function ($line) { return explode(',', $line); },
        file($filename, FILE_IGNORE_NEW_LINES));
}

$distances = readCsv($argv[1]);
$k = $argv[2];

list ($indices, $medoids) = kmedoids($distances, $k);
echo implode(' ', $indices) . "\n";
echo implode(' ', $medoids) . "\n";

これを kmedoids.php として以下のように実行します。引数に距離行列の csv ファイルとクラスタ数を指定します。

$ php kmedoids.php dtw_distance_matrix.csv 5 >kmedoids_result.txt

実行結果は以下のようになります。一行目は、各データ (台風) が何番目のクラスタに割り当てられたかを表します。二行目は、各クラスタについて、入力データの何番目がそのクラスタの medoid かを表します。

$ cat kmedoids_result.txt
1 1 1 1 1 4 1 1 4 1 0 4 2 1 4 1 4 2 1 4 4 2 1 1 3 1 1 3 4 1 0 0 4 1 0 2 0 1 4 1 0 0 2 1 4 1 4 2 2 2 4 3 3 3 2 4 4 4 3 1 1 4 1 3 1 0 4 4 4 2 1 1 3 3 4 3 4 1 4 0 4 1 4 4 2 1 2 4 0 1 0 0 4 0 4 4 3 1 1 3 3 2 3 3 4 4 0 2 4 1 1 1 4 4 1 0 1 1 4 1 0 4 1 1 1 1 1 0 1 1 1 4 0 1 4 4 2 4 2 1 3 2 4 1 1 3 1 3 3 4 1 4 0 1 1 1 0 2 0 1 1 1 4 2 4 2 4 4 1 4 1 1 1 4 1 4 4 1 4 1 1 2 4 1 4 1 4 1 1 2 1 2 1 3 4 1 4 1 1 1 1 1 4 2 4 4 1 0 1 1 1 0 4 4 3 3 3 1 1 4 1 1 4 1 0 4 1 4 1 4 4 4 1 1 1 4 3 1 4 2 1 0 4 4 4 4 1 1 1 1 3 1 4 4 4 1 1 0 1 1 0 0 2 1 1 4 4 4 4 4 1 4 4 1 3 1 1 1 4 4 1 1 0 1 1 1 1 1 2 2 1 2 4 4 1 4 1 4 1 0 1 4 0 2 1 3 1 1 3 2 3 3 1 4 4 3 1 0 1 2 4 1 1 4 4 0 4 1 3 1 3 2 3 3 3 4 3 1 0 1 0 0 0 4 0 2 2 4 1 4 0 1 2 3 4
299 199 48 147 167
クラスタリング結果の可視化

あとは、この結果を可視化すれば冒頭の図が得られます。今回も d3.js を利用しました。クラスタ番号ごとに色を変えて台風の軌跡を描画しています。また、同じ結果をクラスタごとに分けて描画したものが以下の 5 枚の図です。各クラスタの medoid を赤色の軌跡で表しています。

*1:Dynamic Time Warping による時系列データの類似度計算 - y_uti のブログ

*2:「気象庁 | 著作権・リンク・個人情報保護について」の利用規約に従って、自由に利用することが認められています。

*3:http://www.data.jma.go.jp/fcd/yoho/data/typhoon/T0101.pdf

*4:北緯や東経を示す記号は数値の後に出現し、また、直前のデータと同じ場合には省略されます。このため、やや複雑なスクリプトになってしまいました。

*5:緯度の数値を 100 で割った余りとしているのは、たとえば 7112113.9 のように日時と緯度が繋がったデータがあるためです。これは 7 月 11 日 21 時に北緯 13.9 度と読みます。緯度が 10 度未満のデータがあると今回のスクリプトでは処理できませんが、幸い、そのようなデータはなさそうでした。

*6:これは Wikipedia の Algorithms の項に Voronoi iteration method として記載されている方法です。具体的な計算方法には、このほかに PAM, CLARA といった方法があるようです。