Pythonの機械学習ライブラリtslearnを使った時系列データのクラスタリング

こんにちは、ブレインズテクノロジーの柏木です。

今回はPythonで扱える機械学習ライブラリのtslearnを使って、時系列データをクラスタリングしていきたいと思います。

github.com

tslearnとは

時系列分析のための機械学習ツールを提供するPythonパッケージで、scikit-learnをベースとして作られているみたいです。 主な機能として、クラスタリング、教師ありの分類、複数の時系列を重ねた際の重心の計算ができたりします。 今回使用するに至った一番のモチベーションは、波形や振動などの時系列データに対してクラスタリングできるというところです。

tslearnインストール

pipコマンドでインストールできます。

pip install tslearn

Kshapeというクラスタリング手法

今回tslearnで使用するモジュールとして、Kshapeというクラスタリング手法を時系列データに適用していきたいと思います。 Kshapeは2015年に下記の論文で提唱された方法で、以下の流れで実行されるアルゴリズムになります。

1. 相互相関測定に基づいた距離尺度を使う（Shape-based distance: SBD）
2. 1.に基づいた時系列クラスタの重心を計算（時系列の形状を保持する新しい重心ベースのクラスタリングアルゴリズム）
3. 各クラスタに分割する方法は一般的なkmeansと同じだが、距離尺度や重心を計算する時に上記1と2を使う

細かい数式や理論については論文を確認下さい。 個人的に、内容がしっかり記載されていてわかりやすい論文だと思います。

Paparrizos, John and Luis Gravano. “k-Shape: Efficient and Accurate Clustering of Time Series.” SIGMOD Record 45 (2015): 68.

今回使用するデータセットについて

今回使うデータセットはthe UEA & UCR Time Series Classification RepositoryにあるDataset: SonyAIBORobotSurface1を加工したものになります。 このサイトには様々な時系列データがありますので、色々見てみるのも楽しいかもしれないです。

使用したサンプルデータはこちらに置いておきます。 あるワーク単位でファイルが10個あり、1秒単位のセンサーデータを想定したものになります。

時系列データのクラスタリング

それでは、時系列データのクラスタリングをしていきたいと思います。 実際のPythonコードを記載します。 まず、ライブラリのインポートです。

#Python 3.6.4
import pandas as pd
import numpy as np
import glob
from tslearn.clustering import KShape
from tslearn.preprocessing import TimeSeriesScalerMeanVariance

import matplotlib.pyplot as plt
%matplotlib inline

次に、いくつか関数を定義しておきます。

def read_file_to_dataframe(filenames):
    #ファイルデータを読み込み、データフレームを返す
    dfs = []
    for filename in filenames:
        original_df = pd.read_csv(filename, index_col=None, header=0)
        dfs.append(original_df)
    return dfs

def time_series_data_to_array(dataframes, target_col=''):
    #データフレームを読み込み、それらを時系列の配列にする
    tsdata = []
    for i, df in enumerate(dataframes):
        tsdata.append(df[target_col].values.tolist()[:])
        #それぞれの時系列データの最大の長さを確認
        len_max = 0
        for ts in tsdata:
            if len(ts) > len_max:
                len_max = len(ts)
        #時系列データの長さを揃えるために、最後のデータを付け加える
        for i, ts in enumerate(tsdata):
            len_add = len_max - len(ts)
            tsdata[i] = ts + [ts[-1]] * len_add
    
    tsdata = np.array(tsdata)
    return tsdata

def transform_vector(time_series_array):
    #ベクトルに変換
    stack_list = []
    for j in range(len(time_series_array)):
        data = np.array(time_series_array[j])
        data = data.reshape((1, len(data))).T
        stack_list.append(data)
    #一次元配列にする
    stack_data = np.stack(stack_list, axis=0)
    return stack_data


filenames = sorted(glob.glob('sample_data/sample_data*.csv'))
df = read_file_to_dataframe(filenames=filenames)
tsdata = time_series_data_to_array(dataframes=df, target_col='data')
stack_data = transform_vector(time_series_array=tsdata)

ここでは、使用するファイルを読み込んで配列に変換します。 この時に、仮に時系列のデータの長さがファイル毎に異なる場合、長さを揃える必要があります。 そのため、今回は時系列上一番最後の値で埋め合わせを行なっています。 長さを揃えた後は、一次元のベクトルに変換しています。

それでは、加工したデータを学習して実際にクラスタリングしていきます。

seed = 0
np.random.seed(seed)
#相互相関を計算するために、正規化する必要があります。
#TimeSeriesScalerMeanVarianceがデータを正規化してくれるクラスになります。
stack_data = TimeSeriesScalerMeanVariance(mu=0.0, std=1.0).fit_transform(stack_data)

#KShapeクラスのインスタンス化
ks = KShape(n_clusters=2, n_init=10, verbose=True, random_state=seed)
y_pred = ks.fit_predict(stack_data)

#クラスタリングして可視化
plt.figure(figsize=(16,9))
for yi in range(2):
    plt.subplot(2, 1, 1 + yi)
    for xx in stack_data[y_pred == yi]:
        plt.plot(xx.ravel(), "k-", alpha=.2)
    #plt.plot(ks.cluster_centers_[yi].ravel(), "r-")
    plt.title("Cluster %d" % (yi + 1))

plt.tight_layout()
plt.show()

以下、結果をプロットしたものになります。

KShapeクラスの引数n_clustersでクラスタ数を指定することができ、今回はn_clusters=2としました。 実際のデータセットのページを見てみると、AIBOがカーペットを歩いている場合とコンクリートを歩いている場合の二種類のデータとなっています。 それを踏まえて、結果を見てみると良い感じに分離できているように見えます。 今回はクラスタ数を2としましたが、クラスタリングする際によくあるクラスタ数をどのように決めたら良いかという問題について、エルボー法で確認した結果を次に示します。

エルボー法によるクラスターの数の計算

エルボー法とは... 各点から割り当てられたクラスタ中心との距離の二乗の合計をクラスタ内誤差平方和（SSE）として計算します。 クラスタ数を変えて、それぞれのSSE値をプロットし、「肘」のように曲がった点を最適なクラスタ数とする方法です。 また、SSEが小さいほど歪みのない良いモデルであり、うまくクラスタリングできていることになります。

エルボー法について、wikipediaに記載されています。
Elbow method (clustering) - Wikipedia

distortions = []

#1~10クラスタまで計算 
for i in range(1,11):
    ks = KShape(n_clusters=i, n_init=10, verbose=True, random_state=seed)
    #クラスタリングの計算を実行
    ks.fit(stack_data)
    #ks.fitするとks.inertia_が得られる
    #inertia_でSSEを取得できる
    distortions.append(ks.inertia_)

plt.plot(range(1,11), distortions, marker='o')
plt.xlabel('Number of clusters')
plt.ylabel('Distortion')
plt.show()

n_clusters=2で肘のように曲がっているので、今回はクラスタ数2で正解のようです。 ただし、一般的にエルボー法ではこのようにはっきりと曲がりがある場合は少なく、あまり精度はよくないようです。 他にもクラスタ数を決める方法は色々ありシルエット分析、X-means、GAP法などがあります。 興味がある方はX-means、GAP法は論文がありますので、参考文献をご覧下さい。

クラスタリングの決め方について、wikipediaに記載されています。
Determining the number of clusters in a data set - Wikipedia

まとめ

今回はtslearnを使って、時系列データのクラスタリングをしてみました。 振動や波形データをうまくクラスタリングできるので、非常に重宝すると思います。 tslearnにはKshape以外にもTimeSeriesKMeansという手法もあり、距離尺度としてDTW（動的時間伸縮法）やSoft-DTWを使ってクラスタリングをすることができます。さらに時系列に対して重心線を引いてくれたりもします。
今回はデータ点数も少なく振動も緩やかであるため、クラスタ数を決めることが容易でしたが、一般的に時系列データのクラスタ数を決めることは難しい問題だったりするので、その辺りを決める良い方法を模索していきたいと思います。

ブレインズテクノロジーでは「共に成長できる仲間」を募集中です。
採用ページはこちら

参考文献および参考ウェブサイト

• tslearn’s documentation — tslearn 0.1.18.3 documentation
• Silhouette (clustering) - Wikipedia
• k-meansの最適なクラスター数を調べる方法
• Pythonでクラスタリング k-meansからk-medoidsを改良する - 見習いデータサイエンティストの隠れ家
• X-means: Extending K-means with Efficient Estimation of the Number of Clusters
• Estimating the number of clusters in a data set via the gap statistic