Pythonで時系列クラスタリングをする

時系列データを分類したいときに、時系列クラスタリングという方法がある。Pythonにはtslearnというパッケージがあって、k-means法によるクラスタリングができる。距離（類似度）として使えるのはユークリッド距離や動的時間伸縮法 (Dynamic Time Warping: DTW)、Soft-DTW。今回はDTWを使うが、DTWは2つの時系列間の類似度を求める方法のひとつで、期間の長さが異なる時系列データ間の類似度も求めることができる。

Google Trendsでは検索キーワードの週ごとの検索数を調べられるので、このデータを使ってtslearnでクラスタリングしてみる。

環境

Windows10のWSL(Ubuntu 18.04)。
$ lsb_release -a No LSB modules are available. Distributor ID: Ubuntu Description: Ubuntu 18.04.2 LTS Release: 18.04 Codename: bionic $ python3 -V Python 3.6.7 $ pip3 show pandas Name: pandas Version: 0.24.1 $ pip3 show matplotlib Name: matplotlib Version: 3.0.3

時系列データの取得

Google Trendsで次のキーワードで「織田信長」「豊臣秀吉」「徳川家康」「武田信玄」「上杉謙信」「伊達政宗」「毛利元就」「長宗我部元親」「真田幸村」「石田三成」「前田利家」「加藤清正」「黒田官兵衛」を検索した。期間は2018年の1年間。ただし、Google Trendsでは同時に検索できるのは5ワードまでで、結果は最も多い検索数を100とした相対的な数値となる。対象期間中では「織田信長」のときに最高値があるので、以下の3組のようにすべての組に「織田信長」を入れることによって、異なる組でも相対数を比較できるようにした。検索結果はそれぞれbusho1.csv、busyo2.csv、busyo3.csvという名前でcsvとして保存。

• 「織田信長」「豊臣秀吉」「徳川家康」「武田信玄」「上杉謙信」
• 「織田信長」「伊達政宗」「毛利元就」「長宗我部元親」「真田幸村」
• 「織田信長」「石田三成」「前田利家」「加藤清正」「黒田官兵衛」

csvデータの読み込み

保存した3つのcsvを読み込む。1つ目のcsv以外は、「織田信長」を除外する。

import pandas as pd
import matplotlib.pyplot as plt

# ターミナル対応
import matplotlib
matplotlib.use('Agg')

def load_tsdata():
    # csvからデータの読み込み

    # 読み込むファルのリスト
    busyo_l = ['busyo1.csv', 'busyo2.csv', 'busyo3.csv']

    df = pd.DataFrame()
    for i, busyo_csv in enumerate(busyo_l):
        # はじめの2行はスキップ
        # 1列目（週）をインデックスに指定し、datetime型で読み込み
        tmp_df = pd.read_csv(busyo_csv, skiprows=2, parse_dates=True, index_col='週')

        # カラム名から「: (日本)」を除去
        names = tmp_df.columns.tolist()
        tmp_df.columns = [x.replace(': (日本)', '') for x in names]

        # 最初のファイル以外では「織田信長」を除外
        if i != 0:
            tmp_df.drop(columns=['織田信長'], inplace=True)

        df = pd.concat([df, tmp_df], axis=1)

    print(df.index.dtype)
    print(df.dtypes)

    # 読み込んだデータのプロット
    plt.figure()
    df.plot()
    plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
    plt.tight_layout()
    plt.savefig('busho.png')
    plt.clf()

    return df

def main():
    df = load_tsdata()

if __name__ == '__main__':
    main()

読み込んだデータをプロットした図。なんだかんだ言って、織田信長の検索数が多い。

tslearnのインストール

tslearnはpipでインストールできる。
$ sudo pip3 install tslearn

クラスター数の決定

取得した数値のままだと数値の大小が結果に影響しそうで、検索数の増加・減少などの変化でクラスタリングするために、データを0-1.0に正規化してからクラスタリングした。

また、クラスタリングはk-means法なので、クラススター数を決める必要がある。シルエット値によるクラスター数の選択についてはk-means法の最適なクラスター数を選択するを参照。今回は簡易的にシルエット値のみで判定した。

import pandas as pd

from tslearn.clustering import TimeSeriesKMeans, silhouette_score
from tslearn.preprocessing import TimeSeriesScalerMinMax
from tslearn.utils import to_time_series_dataset

def tsclusteringN(ts_data, names):
    # クラスタリング

    # 正規化
    ts_dataset = TimeSeriesScalerMinMax().fit_transform(ts_data)

    metric = 'dtw'
    n_clusters = [n for n in range(2, 6)]
    for n in n_clusters:
        print('クラスター数 =', n)

        # metricが「DTW」か「softdtw」なら異なるデータ数の時系列データでもOK
        km = TimeSeriesKMeans(n_clusters=n, metric=metric, verbose=False, random_state=1).fit(ts_dataset)

        # クラスタリングの結果
        print('クラスタリング結果 =', km.labels_)

        # -1から1の範囲の値。シルエット値が1に近く、かつシルエット値をプロットしたシルエット図でクラスター間の幅の差が最も少ないクラスター数が最適
        # 今回はシルエット値のみを確認
        print('シルエット値 =', silhouette_score(ts_dataset, km.labels_, metric=metric))
        print()

def main():
    df = load_tsdata()
    tsclusteringN(df.values.transpose(), df.columns)

if __name__ == '__main__':
    main()

結果は以下の通り。データ数が少ないためか、クラスター数2のときのシルエット値が最も1に近い。
クラスター数 = 2 クラスタリング結果 = [0 0 0 0 0 0 0 1 0 1 0 1 1] シルエット値 = 0.20001431824275723 クラスター数 = 3 クラスタリング結果 = [0 2 2 0 2 0 2 1 2 1 2 1 1] シルエット値 = 0.13154201699633325 クラスター数 = 4 クラスタリング結果 = [0 2 2 0 2 0 2 3 2 3 2 3 1] シルエット値 = 0.0772278603770549 クラスター数 = 5 クラスタリング結果 = [0 2 4 4 2 0 4 2 2 3 4 3 1] シルエット値 = 0.015061953461705777

時系列クラスタリング

クラスター数を2として時系列クラスタリングして、結果をプロット。

import pandas as pd

from tslearn.clustering import TimeSeriesKMeans, silhouette_score
from tslearn.preprocessing import TimeSeriesScalerMinMax
from tslearn.utils import to_time_series_dataset

import matplotlib.pyplot as plt

# ターミナル対応
import matplotlib
matplotlib.use('Agg')

def plot_clustering(km, ts_dataset, names, n_clusters):
    # クラスタリングの結果をプロット

    # クラスターごとの中心をプロット
    for i, c in enumerate(km.cluster_centers_):
        plt.plot(c.T[0], label=i)

    plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
    plt.show()
    plt.savefig('ts_clust_center.png')
    plt.clf()

    # クラスターごとのプロット
    for i in range(n_clusters):
        for label, d, t in zip(km.labels_, ts_dataset, names):
            if label == i:
                plt.plot(d, label=t)

        plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
        plt.show()
        plt.savefig('ts_labeled{}.png'.format(i))
        plt.clf()

def tsclustering(ts_data, names):
    # 正規化
    ts_dataset = TimeSeriesScalerMinMax().fit_transform(ts_data)

    n_clusters = 2
    metric = 'dtw'

    # metricが「DTW」か「softdtw」なら異なるデータ数の時系列データでもOK
    km = TimeSeriesKMeans(n_clusters=n_clusters, metric=metric, verbose=False, random_state=1).fit(ts_dataset)

    # クラスタリングの結果
    print('クラスタリング結果 =', km.labels_)

    plot_clustering(km, ts_dataset, names, n_clusters)

def main():
    df = load_tsdata()
    #tsclusteringN(df.values.transpose(), df.columns)
    tsclustering(df.values.transpose(), df.columns)

if __name__ == '__main__':
    main()

結果は以下の通り。うまく分かれているのかどうか確信が持てないが・・

クラスター中心のプロット。

クラスター0に分類されたデータ（正規化済み）のプロット。

クラスター1に分類されたデータ（正規化済み）のプロット。