y_uti のブログ

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

k-means 法と DBA で時系列データをクラスタリングする

前回の記事で紹介した DBA を利用することで、DTW を距離尺度として、時系列データに k-means 法を適用できます。そこで今回は、これまでの記事でも利用してきた台風の経路情報を題材に、k-means 法による時系列データのクラスタリングを試してみました。下図は 2001 年から 2015 年までの 355 個の台風を 5 クラスタに分類した結果です。今回の記事の図も、気象庁のウェブページ (http://www.data.jma.go.jp/fcd/yoho/typhoon/index.html) で公開されているデータを加工して作成したものです。

K-means 法はデータを k 個に分類するアルゴリズムです。Wikipedia に説明があります。今回は、クラスタの中心を求める計算に DBA を利用し、各データを最も近い中心のクラスタに割り当て直す処理に DTW を利用します*1
k平均法 - Wikipedia

K-means 法を用いて時系列データをクラスタリングするプログラムを以下のように実装しました。require されている各プログラムは、これまでの記事で作成してきたものです*2。kmeans 関数が wikipediaアルゴリズム全体に対応します。initialize_indices 関数でランダムにクラスタを割り振った後*3クラスタ割り当てが変化しなくなるまで update_centroids と assign_to_nearest を交互に繰り返します。update_centroids から DBA を実行する際には、そのクラスタに割り当てられた時系列のなかで系列長が最大のものを初期値としています。また、DBA 自体が反復処理を行いますが、これを収束するまで繰り返すと大きな時間がかかるので、DBA の反復回数は最大 10 回に制限しました。

<?php
require_once __DIR__ . '/dtw.php';
require_once __DIR__ . '/dba.php';

function kmeans($sequences, $k, $maxiter = 100, $maxdbaiter = 10) {
    $indices = initialize_indices($sequences, $k);
    $centroids = update_centroids($sequences, $indices, $k, $maxdbaiter);
    for ($iter = 0; $iter < $maxiter; ++$iter) {
        $next = assign_to_nearest($sequences, $centroids);
        if ($next === $indices) {
            break;
        }
        $indices = $next;
        $centroids = update_centroids($sequences, $indices, $k, $maxdbaiter);
    }
    return array($indices, $centroids);
}

function initialize_indices($sequences, $k) {
    $indices = range(0, $k - 1);
    for ($i = $k; $i < count($sequences); ++$i) {
        $indices[] = mt_rand(0, $k - 1);
    }
    shuffle($indices);
    return $indices;
}

function assign_to_nearest($sequences, $centroids) {
    $indices = array();
    foreach ($sequences as $s) {
        $mindist = INF;
        $nearest = 0;
        for ($i = 0; $i < count($centroids); ++$i) {
            if (($dtw = dtw($centroids[$i], $s)) < $mindist) {
                $mindist = $dtw;
                $nearest = $i;
            }
        }
        $indices[] = $nearest;
    }
    return $indices;
}

function update_centroids($sequences, $indices, $k, $maxiter) {
    $centroids = array();
    for ($i = 0; $i < $k; ++$i) {
        $s = filter_sequence($sequences, $indices, $i);
        $c = get_longest($s);
        for ($j = 0; $j < $maxiter; ++$j) {
            if (($nc = dba($s, $c)) === $c) {
                break;
            }
            $c = $nc;
        }
        $centroids[] = $c;
    }
    return $centroids;
}

function filter_sequence($sequences, $indices, $target) {
    $filtered = array();
    for ($i = 0; $i < count($indices); ++$i) {
        if ($indices[$i] == $target) {
            $filtered[] = $sequences[$i];
        }
    }
    return $filtered;
}

function get_longest($sequences) {
    $longest = array();
    foreach ($sequences as $s) {
        if (count($s) > count($longest)) {
            $longest = $s;
        }
    }
    return $longest;
}

これを実行するプログラムは次のとおりです。第一引数にクラスタ数を指定し、第二引数以降に対象のファイルを列挙します。結果は、各データが割り当てられたクラスタ番号を標準出力に書き出し、DBA によって計算された各クラスタのセントロイドを centroid_%d.csv という名前のファイルに出力します。

<?php
require_once __DIR__ . '/kmeans.php';

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

function writeSequence($sequence, $filename) {
    $contents = implode("\n", array_map(
        function ($latlng) { return implode(',', $latlng); },
        $sequence));
    file_put_contents($filename, $contents . "\n");
}

$k = $argv[1];
$sequences = array();
for ($i = 2; $i < $argc; ++$i) {
    $sequence = readCsv($argv[$i]);
    $sequences[] = $sequence;
}

list ($indices, $centroids) = kmeans($sequences, $k);

echo implode(' ', $indices) . "\n";
for ($i = 0; $i < count($centroids); ++$i) {
    writeSequence($centroids[$i], sprintf('centroid_%d.csv', $i));
}

以下は、作成したプログラムを実行して結果を確認したものです。5 クラスタに分類しているので、centroid_4.csv までのファイルが出力されます。クラスタ番号ごとに色分けして描画したものが冒頭の図です。

$ php kmeans_main.php 5 T*.csv >indices.txt
$ cat indices.txt
0 0 2 2 0 3 2 0 3 2 3 4 4 2 3 0 3 4 0 3 1 4 2 2 4 2 2 1 4 0 1 1 1 0 3 4 1 2 3 2 3 3 4 2 3 2 3 4 4 4 3 4 4 0 4 0 0 0 2 2 0 1 0 2 1 1 1 3 1 4 0 2 1 1 1 1 0 2 1 0 3 1 3 3 4 0 4 1 3 0 3 3 0 1 1 1 1 2 2 0 1 4 4 2 1 1 0 4 3 2 0 2 3 4 0 1 1 2 3 2 0 1 2 2 2 2 2 1 0 0 2 3 1 0 3 1 4 0 3 2 1 4 3 2 2 2 2 2 4 1 2 1 3 0 0 0 3 4 1 0 2 2 0 4 4 4 3 1 2 4 2 2 2 1 1 1 1 2 0 0 2 4 1 2 0 2 0 2 2 4 2 4 2 1 1 2 0 2 2 2 2 0 3 4 3 1 2 3 2 2 2 1 3 1 2 1 1 2 2 0 2 2 1 0 1 0 0 3 2 1 0 0 2 2 0 3 1 2 1 4 0 3 4 1 3 4 2 2 2 2 2 2 3 1 1 2 2 3 2 0 3 0 4 2 0 1 1 0 3 3 2 1 3 2 2 2 2 2 1 0 2 2 0 2 2 2 2 0 4 4 0 4 0 3 0 3 2 3 1 1 2 3 1 4 2 2 2 2 2 4 2 3 0 0 1 2 0 1 1 4 1 2 0 3 3 1 1 2 2 2 2 4 1 2 4 0 3 2 1 2 3 3 0 3 0 4 4 3 2 3 0 2 4 2 3
$ head -n 5 centroid_0.csv
14.542635658915,133.41240310078
14.790774907749,140.68560885609
17.001976284585,130.03320158103
17.407272727273,125.70545454545
19.70752688172,124.77419354839

クラスタの centroid は以下のとおりです。前回の記事でも確認したように、centroid の軌跡はでこぼこになるようです。


前回の記事で試した、時系列ごとの平均を取ってから centroid を更新する方法で、k-means を実行してみます。クラスタリングの結果は次のようになりました。

この方法で得られた centroid は以下のとおりです。前回と同様に、滑らかな軌跡が得られる一方で、元の時系列に比べて短くなってしまうという特徴が見られました。


*1:各データとクラスタの中心との距離を DTW で計算し、その結果によってデータをクラスタに割り当て直します。

*2:k-medoids 法と DTW による時系列データのクラスタリング - y_uti のブログ, DTW Barycenter Averaging で時系列データの平均を求める - y_uti のブログ を参照してください。

*3:ただし、空のクラスタができないように、最初の $k 個を 0 から $k - 1 に割り当ててから残りをランダムに割り当て、最後に全体をシャッフルして戻しています。