Zenn
📑

JavaScriptで時系列データのクラスタリングがしたい

2021/09/20に公開
JavaScript
TypeScript
OpenLayers
#
clustering
tech

Pythonで時系列解析をしたい(時系列クラスタリング)
DTW(Dynamic Time Warping)で台風軌道をクラスタリングする

こんなふうに時系列クラスタリングをJavaScript（TypeScript）でしてみたい。
でも調べるとほとんどみんなPythonのtslean、sklearnやRを使っています。
調べればJavaScriptライブラリがあるのでしょうが、理屈の勉強がてら自分で実装（TypeScript）してみます。

分類法

分類法においてメジャーなのはK-meansとK-medoids。
2つを比べたとき外れ値やノイズに強いとされている、また時系列データ同士の比較は基本的にはデータ間の距離で行うことから、K-medoidsを今回採用します[1]。

  1. データセットの中からk個のデータをクラスターの代表データ（medoid）としてランダムに選ぶ。
  2. 残ったデータを、“距離”が最も近い代表データを持つクラスターに振り分ける。
  3. 各クラスターで、クラスター内の各データの代表データとの距離の合計が最小になるようmedoidを選び直す。
  4. medoidが変わらなくなるまで2~3を繰り替えす。

ただ、K-meansやK-medoidsには、初期値をランダムに選ぶため毎回結果が異なる、運が悪いと分類がうまくいかないことがある、という課題があります。

それをどうにかしようと先人たちが発明したのがk-means++で、1の初期値を選ぶ工程で、クラスター中心がなるべく固まらないように選択します[2]。
収束しやすくするだけで必ず毎回同じ結果となるわけではありません。

K-meansといっていますが、ようは「選択済みの中心データ（代表データ）から距離が遠いデータが選ばれやすくしよう」という考え方なのでK-medoidsでも使えます。
sklearnでもK-medoidsの初期化にこのロジックを使っています。[3]

function randomIndex(v: number | number[], options?:{
  prng?: ()=>number
} ): number {
  // v is length or array of weight.
  const prng = options?.prng ? options.prng : Math.random;

  if (isNumber(v)) {
    const length = v as number;
    if (length < 1 || !Number.isInteger(length)) {
        throw new Error("length must be positive integer.");
    }
    return Math.floor(prng() * length);
  }

  const weights = v as number[];
  const totalWeight = weights.reduce((sum, weight) => {
    if (weight < 0) {
        throw new Error("weight must be positive value.");
    }
    return sum + weight;
  }, 0);
  const rand = prng() * totalWeight;
  let threshold = 0;
  for (let i = 0; i < weights.length - 1; i++){
    threshold += weights[i];
    if (rand < threshold) return i;
  }
  return weights.length - 1;
};

function initKpp(k: number, distanceMatrix: number[][], prng: ()=>number) {
    const size = distanceMatrix.length;
    const medoidIndexes: number[] = [randomIndex(size, { prng })];

    for (let i = 1; i < k; i++) {
        const candidateIndexes = range(size).filter((index) => {
            return !medoidIndexes.includes(index)
        });
        const dissimilarities = candidateIndexes.map((candidateIndex) => {
            return Math.min(...medoidIndexes.map((medoidIndex) => {
                return distanceMatrix[candidateIndex][medoidIndex] ** 2;
            }))
        });
        const totalSquareDissimilarity = sum(dissimilarities);
        const candidatePotentials = dissimilarities.map((dissimilarity) => {
            return dissimilarity / totalSquareDissimilarity;
        });
        medoidIndexes.push(candidateIndexes[randomIndex(candidatePotentials, { prng })]);
    }
    medoidIndexes.sort(function (a, b) {
        if (a < b) return -1;
        if (a > b) return 1;
        return 0;
    });
    return medoidIndexes;
}

※prngはランダムな数を返す関数

距離

次に、クラスターを分けるのに利用する距離について。
点同士の距離だったらシンプルにユークリッド距離を使うのですが、時系列データ同士の距離となるとそうもいきません。

ここではDynamic Time Warpingを使います[4][5][6]。
この方法だと長さが揃っていない時系列データ同士の比較も可能です。

type DtwValue = { cost: number, path: [number, number] };

function isNumber(n: any): boolean {
  return typeof n === "number" && isFinite(n);
};

function euclideandistance(a: number | number[], b: number | number[]): number {
    const v1 = isNumber(a) ? [a as number] : a as number[];
    const v2 = isNumber(b) ? [b as number] : b as number[];

    if (v1.length !== v2.length) {
        throw new Error("arguments must be same length.");
    }
    let distance = 0;
    for (let i = 0; i < v1.length; i++){
        distance += (v1[i] - v2[i]) ** 2;
    }

    return Math.sqrt(distance);
}

function minVal (values:DtwValue[]): DtwValue{
    const sorted = [...values].sort(function(a,b){
        if( a.cost < b.cost ) return -1;
        if( a.cost > b.cost ) return 1;
        return 0;
    });
    return sorted[0]
}

function dtw(
  ts1: number[] | number[][],
  ts2: number[] | number[][],
): number {
    if (!Array.isArray(ts1) || ts1.length < 1) {
        throw new Error("The first argument is not array.");
    }
    if (!Array.isArray(ts2) || ts2.length < 1) {
        throw new Error("The second argument is not array.");
    }

    const ret: DtwValue[][] = [[{ cost: euclideandistance(ts1[0], ts2[0]), path: [-1, -1] }]];
    for (let i = 1; i < ts1.length; i++){
        ret.push(
            [{ cost: ret[i - 1][0].cost + euclideandistance(ts1[i], ts2[0]), path: [i - 1, 0]}]
        )
    }
    for (let j = 1; j < ts2.length; j++){
        ret[0][j] = { cost: ret[0][j - 1].cost + euclideandistance(ts1[0], ts2[j]), path: [0, j - 1] };
    }
    for (let i = 1; i < ts1.length; i++) {
        for (let j = 1; j < ts2.length; j++) {
            const min = minVal([
                {
                    cost: ret[i - 1][j].cost,
                    path: [i - 1, j]
                },
                {
                    cost: ret[i][j - 1].cost,
                    path: [i, j - 1]
                },
                {
                    cost: ret[i - 1][j - 1].cost,
                    path: [i - 1, j- 1]
                }]
            )
            ret[i][j] = {
                cost: min.cost + euclideandistance(ts1[i], ts2[j]),
                path: min.path
            }
        }
    }

    return ret[ts1.length - 1][ts2.length - 1].cost
}

文字列の類似度を測るのに使われるレーベンシュタイン距離の考え方をもとにしているそうです。

K-medoids実装

 type FitOptions = {
    maxIter?: number;
    initialize?: number;
    prng?: () => number;
}  

function sum(values: number[]): number {
  return values.reduce((sum, weight) => {
    return sum + weight;
  }, 0);
}

function zip<S, T>(arr1: S[], arr2: T[]): [S, T][] {
    if (!Array.isArray(arr1) || !Array.isArray(arr2) || arr1.length !== arr2.length) {
      throw new Error("arguments must be same length array.");
    }

    const length = arr1.length;
    const ret: [S, T][] = [];
    for (let i = 0; i < length; i++){
        ret.push([arr1[i], arr2[i]])
    }
    return ret;
}

function sse(clusters: number[][], medoidIndexes: number[], distanceMatrix: number[][]):number {
    let ret = 0;
    zip(clusters, medoidIndexes).forEach(([clusterIndexes, medoidIndex]) => {
        clusterIndexes.forEach((clusterIndex) => {
            ret += distanceMatrix[medoidIndex][clusterIndex] ** 2
        })
    })
    return ret;
}

function classificate(medoidIndexes: number[], distanceMatrix: number[][]) {
    const clusters: number[][] = medoidIndexes.map(() => []);
    distanceMatrix.forEach((distances, index) => {
        let minDistance = Infinity;
        let minIndex = -1;
        medoidIndexes.forEach((medoidIndex) => {
            const d = distances[medoidIndex];
            if (d < minDistance) {
                minDistance = d;
                minIndex = medoidIndex;
            }
        })
        clusters[medoidIndexes.indexOf(minIndex)].push(index)
    });
    return clusters;
}

class Kmedoids {
    size: number;
    distances: number[][];

    constructor(values: number[][] | number[][][]) {
        const size = values.length;
        const distances: number[][] = values.map(() => {
            return Array.from({ length: size }, () => 0)
        });
        for (let i = 0; i < size; i++) {
            for (let j = i + 1; j < size; j++) {
                const d = dtw(values[i], values[j]);
                distances[i][j] = d;
                distances[j][i] = d;
            }
        }

        this.size = size;
        this.distances = distances;
    }

    fit(
        k: number,
        options: FitOptions,
    ): Promise<[number[][], number[]]> {
        return new Promise((resolve, reject) => {
            if (this.size < k) {
                reject(new Error("k is too large."));
            }

            const maxIter = isNumber(options?.maxIter) ? options.maxIter as number : 100;
            const initialize = isNumber(options?.initialize) ? options.initialize as number : 5;
            const prng = options?.prng ? options.prng : random();
        
            let medoidIndexes = initKpp(k, this.distances, prng);
            let clusters: number[][] = classificate(medoidIndexes, this.distances);
            let initScore = sse(clusters, medoidIndexes, this.distances);
            let inited = 1;
            while (inited < initialize) {
                const candidateMedoidIndexes = initKpp(k, this.distances, prng);
                const candidateClusters: number[][] = classificate(candidateMedoidIndexes, this.distances);
                const candidateInitScore = sse(candidateClusters, candidateMedoidIndexes, this.distances);

                if (candidateInitScore < initScore) {
                    medoidIndexes = candidateMedoidIndexes;
                    clusters = candidateClusters;
                    initScore = candidateInitScore;
                }
                inited++;
            }

            let cnt = 0;
            while (cnt++ < maxIter) {
                const newMedoidIndexes: number[] = [];
                clusters.forEach((clusterIndexes) => {
                    let newMedoid = -1;
                    let minVariance = Infinity;
                    clusterIndexes.forEach((index) => {
                        const intraDistances = clusterIndexes.map((otherIndex) => {
                            return this.distances[index][otherIndex]
                        })
                        const variance = sum(intraDistances);
                        if (variance < minVariance) {
                            minVariance = variance;
                            newMedoid = index;
                        }
                    })
                    newMedoidIndexes.push(newMedoid);
                })
                newMedoidIndexes.sort(function (a, b) {
                    if (a < b) return -1;
                    if (a > b) return 1;
                    return 0;
                });
                if (medoidIndexes.toString() === newMedoidIndexes.toString()) break;
                medoidIndexes = newMedoidIndexes;
                clusters = classificate(newMedoidIndexes, this.distances);
            }
            if (cnt++ === maxIter) reject(new Error());

            resolve([clusters, medoidIndexes])
        });
    }
}

kの決定

K-meansにしろK-medoidsにしろ、K-＊＊＊法はクラスター数kを実行者が決めないといけません。
そんな都合よく事前にKがわかってることなんて普通は稀で、k=2~テキトウな数で試して最も当てはまりの良いものを選ぶ必要があります。

選び方としてはエルボー法というグラフの傾きを見ながら最適なクラスター数を推定する方法がよく紹介されていますが、自分でやってみても正直よくわからない。差が微妙なときに目で見て判断できる気がしない。

CVI

一般に、

  • Intra Cluster Distance （クラスター内距離）が最小（ぎゅっとまとまっている）
  • Inter Cluster Distance （クラスター間距離）が最大（クラスターが十分に離れている）

ほど良いクラスターであると言われます[7]。

その"良さ"を測る指標(Cluster Validity Index)は様々な人によって開発されており[8]、その中でもデータ間の距離だけで計算ができるシルエット法（Silhouette index）とXie-Beni index、Davies–Bouldin indexを今回は試してみます[9]。

Silhouette index[10]

a_i = 自身が所属するクラスター内の他のデータとの距離の平均
d_{ij} = 自身が所属するクラスター以外のクラスターjの各データとの距離の平均
b_i = d_{ij}の最小値（=最も近所のクラスターの各データとの距離の平均）
i = 1, ..., n(データ数), j = 1, ... , nc(クラスター数)

s_i = (b_i - a_i) / max(a_i, b_i)
ただしクラスター内のデータがa_iのみの場合はs_i=0

クラスターごとにs_iの平均ave(s_j)を求める。
ave(s_j)を全クラスターで合計し、クラスター数で割る。

この値が最大となるクラスター数をkとして選択する。

function silhouette(args: {
    clusters: number[][];
    distanceMatrix: number[][];
}) {
    const {clusters, distanceMatrix } = args;
    
    function ai(i: number, clusterIndexes: number[], distanceMatrix: number[][]): number {
        let ret = 0;
        clusterIndexes.forEach((clusterIndex) => {
            ret += distanceMatrix[i][clusterIndex];
        });
        ret /= (clusterIndexes.length - 1);
        return ret;
    }

    function bi(i: number, clusters: number[][], distanceMatrix: number[][]): number {
        let ret = Infinity;

        clusters.forEach((clusterIndexes) => {
            if (!clusterIndexes.includes(i)) {
                let value = 0;
                if (clusterIndexes.length > 0)  {
                    clusterIndexes.forEach((clusterIndex) => {
                        value += distanceMatrix[i][clusterIndex];
                    })
                    value /= clusterIndexes.length;
                }
                if (value < ret)
                    ret = value;
            }
        })
        return ret;
    }
    const values: number[] = [];
    clusters.forEach((clusterIndexes) => {
        let value = 0;
        if (clusterIndexes.length > 1) {
            clusterIndexes.forEach((clusterIndex) => {
                const a = ai(clusterIndex, clusterIndexes, distanceMatrix);
                const b = bi(clusterIndex, clusters, distanceMatrix);
                value += (b - a) / Math.max(a, b);
            })
            values.push(value / clusterIndexes.length);
        }
    })
    return sum(values) / clusters.length;
}

Xie-Beni index[11]

クラスター内誤差平方和(SSE)とクラスター間の距離（medoidでは代表データ間の距離）から計算します。

SSE = \sum^{nc}_{i=1} \sum^{}_{x\in{C_i}} dist(x - m_i)^2
D_{ij} = dist(m_i - m_j)^2
XB(nc) = \frac{(SSE / n)}{min_{i,j\neq i}(D_{ij})}

XBが最小となるncをkとして選択する。

XB*

クラスター内誤差平方の平均（SSE / n）でなく、各クラスターのクラスター内誤差をクラスター内データ数で割った値を評価に使う改良型XBが提案されています[12]。

S_i = \frac{\sum^{}_{x\in{C_i}} dist(x - m_i)^2}{len(C_i)}
D_{ij} = dist(m_i - m_j)^2
XB^*(nc) = \frac{max_{i=1,...,nc}(S_i)}{min_{i,j\neq i}(D_{ij})} 
function xb(args: {
    clusters: number[][];
    medoidIndexes: number[];
    distanceMatrix: number[][];
}) {
    const { clusters, medoidIndexes, distanceMatrix } = args;
    let separation = Infinity;
    for (let i = 0; i < medoidIndexes.length; i++) {
        for (let j = i + 1; j < medoidIndexes.length; j++) {
            const candidate = distanceMatrix[medoidIndexes[i]][medoidIndexes[j]] ** 2;
            if (candidate < separation)
                separation = candidate;
        }
    }
    let compactness = 0;
    zip(clusters, medoidIndexes).forEach(([clusterIndexes, medoidIndex]) => {
        let candidate = 0;
        clusterIndexes.forEach((clusterIndex) => {
            candidate += distanceMatrix[medoidIndex][clusterIndex] ** 2
        });
        candidate /= clusterIndexes.length;
        if (candidate > compactness)
            compactness = candidate;
    });
    return compactness / separation;
}

Davies–Bouldin index[11:1]

この指標も見積もり方が違うだけでクラスター内距離とクラスター間距離から求めます。

R_{ij} = \frac{(S_i + S_j)}{D_{ij}}
DB(nc) = \frac{\sum^{cn}_{i=1} max_{j=1,...,nc, j\neq i}(R_{ij})}{nc}

DBが最小となるncをkとして選択する。

DB*

こちらも改良型が提案されています[13]。

DB^*(nc) = \frac{\sum^{cn}_{i=1} \frac{max_{k=1,...,nc, k\neq i}(S_i + S_k)}{min_{l=1,...,nc, l\neq i}(D_{il})}}{nc} 
function db(args: {
    clusters: number[][];
    medoidIndexes: number[];
    distanceMatrix: number[][];
}) {
    const { clusters, medoidIndexes, distanceMatrix } = args;
    const averageIntraDistances = zip(clusters, medoidIndexes).map(([clusterIndexes, medoidIndex]) => {
        let distance = 0;
        clusterIndexes.forEach((clusterIndex) => {
            distance += distanceMatrix[medoidIndex][clusterIndex]
        });
        return distance / clusterIndexes.length;
    });
    const k = averageIntraDistances.length;
    let ret = 0;
    for (let i = 0; i < k; i++){
        let otherAverageIntraDistance = 0, minInterClusterDistance = Infinity;
        for (let j = 0; j < k; j++) {
            if (i !== j) {
                if(averageIntraDistances[j] > otherAverageIntraDistance)
                    otherAverageIntraDistance = averageIntraDistances[j];
                if (distanceMatrix[medoidIndexes[i]][medoidIndexes[j]] < minInterClusterDistance)
                    minInterClusterDistance = distanceMatrix[medoidIndexes[i]][medoidIndexes[j]];
            }
        }
        ret += (averageIntraDistances[i] + otherAverageIntraDistance) / minInterClusterDistance;
    }
    return ret / k;
}

テスト

type OptimizeOptions = FitOptions & {
    maxCluster?: number;
    cvi?: "xb" | "db" | "silhouette";
};

async function optimize(
    kmedoids: Kmedoids,
    options: OptimizeOptions,
): Promise<{
    k: number;
    clusters: number[][];
    medoids: number[];
    scores: { k: number, score: number }[];
}>{
    const max = (typeof options?.maxCluster !== "undefined" && isNumber(options.maxCluster) && options.maxCluster > 2) ?
        (options.maxCluster < kmedoids.size ? options.maxCluster : kmedoids.size) : Math.floor(kmedoids.size / 2);
    const ks = range(2, max + 1);

    let cvi: (args: cviOptions) => number, minimize: boolean;
    switch (options.cvi) {
        case ("xb"):{
            cvi = xb;
            minimize = true;
            break;
        }
        case ("db"):{
            cvi = db;
            minimize = true;
            break;
        }
        default: {
            cvi = silhouette;
            minimize = false;
        }
    }

    let bestScore = minimize ? Infinity : -Infinity, bestMedoids: number[] = [], bestClusters: number[][] = [], bestK = -1;
        
    const scores: { k: number, score: number }[] = [];
    for (let i = 0; i < ks.length; i++) {
        const [clusters, medoidIndexes] = await kmedoids.fit(ks[i], options);
        const score = cvi({
            clusters,
            medoidIndexes,
            distanceMatrix: kmedoids.distances
        });
        scores.push({
            k: ks[i],
            score
        });
        if (minimize ? score < bestScore : score > bestScore) {
            bestScore = score;
            bestMedoids = medoidIndexes;
            bestClusters = clusters;
            bestK = ks[i];
        }
    }
    return {
        scores,
        k: bestK,
        medoids: bestMedoids,
        clusters: bestClusters
    }
}

1次元

ここを参考に、ランダム、線形、sin、cosの４種類のデータを作って分類してみます。

function range(a, b) {
  let start = Number.isInteger(b) ? a : 0;
  let end = (Number.isInteger(b) ? b : a) as number;

  if (!Number.isInteger(start) || !Number.isInteger(end)) {
    throw new Error("arguments must be integer.");
  }
  if (start > end) [start, end] = [end, start];
  return Array.from({ length : end - start }, (_, i) => i + start);
}

function rand(length) {
    const ret = [];
    range(length).forEach(() => {
	ret.push(Math.random());
    });
    return ret;
}
function linear(length) {
    const ret = [];
    range(length).forEach((i) => {
	ret.push( i  * 0.01 + Math.random() * 0.1);
    });
    return ret;
}
function sin(length) {
    const ret = [];
    range(length).forEach((i) => {
	ret.push(Math.sin(i / 100 * 2 * Math.PI * 2) + 0.3 * Math.random());
    });
    return ret;
}
function cos(length) {
    const ret = [];
    range(length).forEach((i) => {
	ret.push(Math.cos(i / 100 * 2 * Math.PI) + 0.3 * Math.random());
    });
    return ret;
}

長さ100のデータを5つ*4種用意します。
0 - 4 : ランダム
5 - 9 : 線形
10 - 14 : sin
15　- 19 : cos

最大10クラスで最適化してみます。

シルエット法


シルエット法によるとベストのkは3でした。

[
  [0, 1, 2, 3, 4, 5, 6, 7, 8, 9],
  [10, 11, 12, 13, 14],
  [15, 16, 17, 18, 19]
]

どうやらランダムも線形も同じクラスタということらしいです。

XB*


XBではなんと8クラスがベストという結果に。

[
  [0],
  [1],
  [2],
  [3],
  [4],
  [5, 6, 7, 8, 9],
  [10, 11, 12, 13, 14],
  [15, 16, 17, 18, 19],
]

シルエット法と比べて、ランダムと線形の区別はできていますが、今回の例ではXBはランダムを全て別のクラスとみなすようです。

DB*


DBはXBと同じ結果となりました。

地理座標

台風をその進路から分類してみます。データは気象庁の台風位置表　令和2年（2020年）より入手しました。

この年は23号まで発生しています。
これを最大11クラスで最適化してみます。

JSONへの加工(Python)

import pandas as pd
import json

df = pd.read_csv("table2020.csv",encoding='shift_jis')
names = df['台風名'].unique()
ret = []
for name in names:
    typhoon = {
        'name': name,
        'coordinates': []
    }
    for index, row in df[df['台風名'] == name].iterrows():
        typhoon['coordinates'].append([row['経度'],row['緯度']])
    ret.append(typhoon)

with open("typhoon.json", 'w') as outfile:
    json.dump(ret, outfile, indent=4)

Openlayersでの表示

<!DOCTYPE HTML>
<html>
<head>
    <meta charset="utf-8">
    <meta http-equiv="X-UA-Compatible" content="IE=edge">
    <meta name="viewport" content="width=device-width, initial-scale=1.0">
    <title>tcjs sample</title>

    <link rel="stylesheet" href="https://cdn.jsdelivr.net/gh/openlayers/openlayers.github.io@master/en/v6.6.1/css/ol.css" type="text/css">
    <style>
      .map {
        height: 400px;
        width: 100%;
      }
    </style>

    <script src="https://cdn.jsdelivr.net/gh/openlayers/openlayers.github.io@master/en/v6.6.1/build/ol.js"></script>
</head>

<body>
    <div id="default" class="map"></div>
    <script type="text/javascript">
      const typhoon = [
	    //略, typhoon.jsonの内容
      ];

      const sourceLayer = new ol.layer.Vector({
        source: new ol.source.Vector({
            features: new ol.format.GeoJSON().readFeatures(
              {
                type: "FeatureCollection",
                features: typhoon.map(item=>{
                  return {
                      type: "Feature",
                      geometry: {
                        type: "LineString",
                        coordinates: item.coordinates,
                      },
                      properties: {},
                  }
                }),
              },
              {
                dataProjection: "EPSG:4326",
                featureProjection: "EPSG:3857",
              }
            ),
            attributions: [
              "気象庁「<a href='http://www.data.jma.go.jp/fcd/yoho/typhoon/position_table/table2020.html'>台風位置表　令和2年（2020年）</a>」を加工して作成"
            ]
        }),
        style: new ol.style.Style({
          stroke: new ol.style.Stroke({
            color: 'green',
            width: 2,
          }),
        }),
      });
      new ol.Map({
        target: 'default',
        layers: [
          new ol.layer.Tile({
            source: new ol.source.OSM()
          }),
          sourceLayer
        ],
        view: new ol.View({
          center: ol.proj.fromLonLat([135, 30]),
          zoom: 3
        }),
      });

      const randomColors = (num)=>{
          const diff = 360 / num;
          const colors =[];
          const base = Math.floor(prng() * 360);
          for (let i = 0; i < num; i++) {
            const color = (base + i * diff) % 360;
            colors.push("hsl(" + color + ", 100%, 50%)");
          }

          return colors;
      }
      const testdata = typhoon.map((item)=>{
        return item.coordinates
      })

      const kmedoids = new Kmedoids(testdata);
      const maxK = Math.floor(testdata.length/2);
      const cvi = "silhouette";
      tcjs.kmedoids.optimize(kmedoids, {
          prng,
          cvi,
          maxCluster: maxK
      }).then((result)=>{  
          const layers = [new ol.layer.Tile({
            source: new ol.source.OSM()
          })];

          const colors = randomColors(result.k);
          result.clusters.forEach((clasterIndexes,i)=>{
              const sourceLayer = new ol.layer.Vector({
                  source: new ol.source.Vector({
                      features: new ol.format.GeoJSON().readFeatures(
                        {
                          type: "FeatureCollection",
                          features: clasterIndexes.map(i=>{
                            return {
                                type: "Feature",
                                geometry: {
                                  type: "LineString",
                                  coordinates: testdata[i],
                                },
                                properties: {
                                  medoid: result.medoids.includes(i)
                                },
                            }
                          }),
                        },
                        {
                          dataProjection: "EPSG:4326",
                          featureProjection: "EPSG:3857",
                        }
                      ),
                      attributions: [
                        "気象庁「<a href='http://www.data.jma.go.jp/fcd/yoho/typhoon/position_table/table2020.html'>台風位置表　令和2年（2020年）</a>」を加工して作成"
                      ]
                  }),
                  style : (feature)=>{
                    return new ol.style.Style({
                      stroke: new ol.style.Stroke({
                        color: colors[i],
                        width: feature.getProperties()["medoid"]? 5:2,
                      }),
                    })
                  },
              });
              layers.push(sourceLayer);
          })

          new ol.Map({
            target: 'result',
            layers,
            view: new ol.View({
              center: ol.proj.fromLonLat([135, 30]),
              zoom: 3
            }),
          });
        });	    
    </script>
</body>
</html>

シルエット法

[
  [3, 4, 7, 8, 9, 11, 12, 13],
  [0, 1, 2, 5, 6, 10, 14, 15, 16, 17, 18, 19, 20, 21, 22]
]


太い線はmedoidです。
西へ抜けていくグループと北へ向かうグループの2つという、あまり面白くない結果となりました。

XB*


ぱっとみわかりにくいですが、k=6が最小となります。

[
  [3, 4],
  [7, 8],
  [9],
  [11, 13],
  [12],
  [0, 1, 2, 5, 6, 10, 14, 15, 16, 17, 18, 19, 20, 21, 22]  
]

冒頭で表示していた図です。
データが1つしかないクラスターがいくつかあり、やや分け過ぎな感もあります。
が、medoidを眺めるとどれも確かに特徴的な経路をしており、妥当な気もしてくるので難しいです。
なお2番目にscoreが小さいのはk=2でした。

DB*


k=6が最小となり結果だけ見ればXBと同じですが、scoreが僅差だったXBに対してDBでは明確に差が出ました。

[
  [3, 4],
  [7, 8],
  [9],
  [11, 13],
  [12],
  [0, 1, 2, 5, 6, 10, 14, 15, 16, 17, 18, 19, 20, 21, 22]
]

クラスターの分類結果もXBと同じでした。

まとめ

先人の知恵を借りつつ、JavaScript（TypeScript）で時系列データの分類を行いました。
tcjs

あくまで初歩であり、より精度を高めたければ以下の改善が考えられます。

  • DTW以外の距離を試す（DDTW, DBAなど）
  • ほかのCVIを導入する
  • k-means++以外の初期化法を導入する

残課題は

  • ざっと作っただけなので計算の高速化、省メモリ化
  • irisデータセットみたいな評価用の時系列データって何かないのか

既存のライブラリを使えばサクッとできることですが、ありものを使うだけではk-means++やDTW、各CVIについて詳しく知ることはなかったと思うので(特にCVIは日本語だとエルボー法やX-meansばかり検索にひっかかるので)、たまには車輪の再発明をするのも悪くないです。

脚注

  1. Improved k-medoids clustering based on cluster validity index and object density ↩︎

  2. k-means++: the advantages of careful seeding ↩︎

  3. sklearnなど各ライブラリではさらにクラスター中心候補を複数生成し、クラスター内誤差平方和が最も小さいものを初期値として採用する、という工夫を行っています。 ↩︎

  4. A global averaging method for dynamic time warping, with applications to clustering ↩︎

  5. 日本語のわかりやすい解説：DTW(Dynamic Time Warping)動的時間伸縮法 ↩︎

  6. DTW以外の距離尺度の参考：Local Shapeletを用いた時系列分類に最適な距離尺度の選択 ↩︎

  7. https://web.iitd.ac.in/~bspanda/clustering.pdf ↩︎

  8. Understanding of Internal Clustering Validation Measures ↩︎

  9. X-meansというK-meansのKを自動で決定するアルゴリズムではBICが指標として使われていますが、時系列データにおいてBICがどう見積もるのが正しいのかは勉強不足でよくわからず。 ↩︎

  10. Silhouettes: A graphical aid to the interpretation and validation of cluster analysis ↩︎

  11. New indices for cluster validity assessment ↩︎ ↩︎

  12. 同論文内では、最適なncの前後でIntra Cluster Distanceが急激に減少するという仮定のもとさらに拡張したXB** も提案されています。 ↩︎

  13. XB*同様、さらに拡張したDB** が提案されています。 ↩︎

Discussion