Zenn
👌

telearnを使って気象データをクラスタリングしてみる

2021/12/13に公開
Python
#
クラスタリング
tech

この記事は建築環境/設備 Advent Calendar 2021の13日目の記事です。

tslearnというライブラリを使って時系列クラスタリングで気象データをグループ分けしてみます。気象データ以外でも負荷データの分析などでも使えると思いますが、とりあえずやってみたかったのでやることにしました。今回はtslearnで扱うk-means法やそれに類似した手法を使い、その他のクラスタリングには触れません。

手法の説明はふわっとしているので、より深く正確な内容を学びたい方はクラスタリングについて書かれている本やtslearnやscikit-learnのドキュメントを読んだり、そこで紹介されている論文に当たるとよいと思います。Colabで実行したPythonコードをここにおいています。

クラスタリングとは

クラスタリングとは与えられたデータを類似度に基づいて自動的に分類する手法です。機械学習においてはグループの分類に正解ラベルがあってそれを学習する事を「分類(classification)」といい教師あり学習になりますが、「クラスタリング(clustering)」は教師なし学習です。一般にクラスタリングの結果に正解はありません。分類の意味付けなどの検証作業は各々が別途行います。

時系列ではないデータのクラスタリング（k-means法）

時系列クラスタリングをする前に、まずは一般的な非階層型のクラスタリングとしてk-means法について確認します。k-means法はシンプルなアルゴリズムで以下の手順でk個のグループに分けます。

  • 最初に各データをランダムにk個のグループに割り振ります。
  • あとはグループの割り当てが収束するまで以下を繰り返します。
  • 1.各グループの重心(centroid)を計算します。
  • 2.各点と重心との距離(類似度)を計算し、最も距離が近いグループに再割り当てします。
  • 収束しない(各ステップの計算でグループが変わる数が一定値以下にならない)場合は設定した計算ステップを終えたところで計算を打ち切ります。

グループ数kは結果を見ながら自分で適当な値を決めます。エルボー法やシルエット分析などの結果を参考にしてもよいですが、自動的に決めるのは難しい事が多いようです。英語版のwikipediaの記事に目を通しておくと良いかもしれません。

k-means法では最初に重心を計算したときにたまたま重心同士が近いとうまくグループ分けできません。この問題に対処するため、重心の初期値を決める際にお互いが遠いところに配置されるような確率が高くなるようにしたものをk-means++法といい、通常はこちらを使います。scikit-learnのk-meansのデフォルトもこれです。

試しに適当なデータをk-means++法でクラスタリングしてみます。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans

# データ生成
np.random.seed(seed=42)
x = np.random.rand(30)
y = np.random.rand(30)

# 5つのグループにクラスタリング
n = 5
cust_array = np.array([x.tolist(), y.tolist()]).T
pred = KMeans(n_clusters=n).fit_predict(cust_array)
for i in range(n):
    plt.scatter(x[pred==i], y[pred==i])

5つのグループに色分けできました。年データや日データなどの特徴量を用いて分類するにはこれを使ってもよいと思います。他にもscikit-learnにいろいろなクラスタリング手法が実装されているので試してみても良いと思います。

データの次元数はいくつに増えても良いのですが、前処理にせよ結果の確認にせよPCA(主成分分析)を一緒に使って次元削減する事も多いようです。

時系列データ用のライブラリtslearnとクラスタリング用気象データの用意

今回、時系列クラスタリングに使用するライブラリとしてはtslearnを使用します。通常のクラスタリングとの違いとしては時系列データ同士の類似度をどう定義し計算するかというところが大きいと思うのですが、tsが（おそらく）time-seriesの略であることもあり、時系列解析用の機能がまとめられています。tslearnはscikit-learn、numpy、scipyなどを使って作られています。インストールはpipで行えます。

!pip install tslearn

気象データにはEnergyplusのデータを使います。とりあえず大阪の気象データをPandasのDataFrameに読み込みます。使わないデータも一応ラベルを付けておきます。もっと短い名前の方が使いやすいと思いますが、今回は一度変数に取り込んだら触らないので内容が分かる長い名前にしておきます。

!wget https://energyplus-weather.s3.amazonaws.com/asia_wmo_region_2/JPN/JPN_Osaka.477710_IWEC/JPN_Osaka.477710_IWEC.epw

import pandas as pd
pd.set_option('display.max_columns', 35)

epwfile = "./JPN_Osaka.477710_IWEC.epw"
labels = ["Year(-)","Month(-)","Day(-)","Hour(-)","Minute(-)","Data Source and Uncertainty Flags(-)",
          "Dry Bulb Temperature(degC)","Dew Point Temperature(degC)","Relative Humidity(%)",
          "Atmospheric Station Pressure(Pa)","Extraterrestrial Horizontal Radiation(Wh/m2)",
          "Extraterrestrial Direct Normal Radiation(Wh/m2)","Horizontal Infrared Radiation from Sky(Wh/m2)",
          "Global Horizontal Radiation(Wh/m2)","Direct Normal Radiation(Wh/m2)",
          "Diffuse Horizontal Radiation(Wh/m2)","Global Horizontal Illuminance(lx)",
          "Direct Normal Illuminance(lx)","Diffuse Horizontal Illiminance(lx)","Zenith Luminance(cd/m2)",
          "Wind Direction(deg)","Wind Speed(m/s)","Total Sky Cover(0-10)","Opaque Sky Cover(-)",
          "Visibility(km)","Ceiling Height(m)","Present Weather Observation(-)","Present Weather Code(-)",
          "Precipitable Water","Aerosol Optical Depth(-)","Snow Depth",
          "Days Since Last Snowfall day(day)","Albedo(-)","liquid precipitation depth(mm)",
          "liquid Precipitation quantity(h)"]
df_epw = pd.read_csv(epwfile, header=7, names=labels)

乾球温度のトレンドグラフを書いてデータをチェックしてみましょう。横軸は時間(hour)で一週間分を一つの折れ線にしており、あまり重なると見づらいので月毎にグラフを作成しています。

import datetime

fig, axes = plt.subplots(12, figsize=(8.0, 18.0))
plt.subplots_adjust(hspace=0.5)
for i in range(12):
    data = df_epw[df_epw['Month(-)']==i+1]['Dry Bulb Temperature(degC)'].values
    begin = 0
    end = 24*7 - 1
    ax = axes[i]
    while begin < data.size - 1:
        ax.plot(data[begin:end])
        ax.set_ylim(-5, 40)
        ax.text(-5, 30, datetime.datetime(2021, i+1, 1).strftime('%b'))
        begin = end + 1
        end = min(end+24*7, data.size-1)

気温の動きは似たような形をしていることが多いのが分かります。正直もうちょっと変な動きがたくさんあるかと思っていましたがこれでやってみましょう。

データをtslearnで使える形式に整理します。今回は1時間データ24個（1日分）を1つの時系列データとして扱うことにし、年間データをまとめてやると見づらそうなので12月だけを抽出します。なお、tslearnで使うデータはnumpy.ndarrayである必要があります。

import calendar

month = 12
daynum = calendar.monthrange(2021, month)[1]
db = np.reshape(df_epw[df_epw['Month(-)']==month]['Dry Bulb Temperature(degC)'].values, (daynum, 24, 1))

tslearnにはデータをテキストから読み込む関数もありますが独特なフォーマットなので、任意の形式のデータを読み込んで自分で変換してやるのが良い気がします。今回は使わなかったのですがtslearn.utils.to_time_seriesやtslearn.utils.to_time_series_datasetなどの関数も用意されているので利用してみてもよいかもしれません。

クラスタリングの前処理

クラスタリングをする前にtslearn.preprocessingにあるツールを使って時系列データの前処理をします。次のようなツールがあります。

  • TimeSeriesScalerMeanVariance:標準化に使う。平均(デフォルト0)と分散(デフォルト1)を指定してスケーリングする
  • TimeSeriesScalerMinMax:正規化に使う。最小値(デフォルト0)と最大値(デフォルト1)を指定してスケーリングする
  • TimeSeriesResampler:元のデータを線形補間して指定した要素数のデータをとりだす

今回はとりあえずTimeSeriesScalerMeanVarianceを使って平均0と分散1にデータをスケーリングします。

from tslearn.preprocessing import TimeSeriesScalerMeanVariance

# クラスタリングするデータ
data = db   
print(np.max(data),np.min(data))

# 前処理として標準化を行う
nm_data = TimeSeriesScalerMeanVariance().fit_transform(data)
print(np.max(nm_data),np.min(nm_data))

クラスタ数をとりあえず3に決めておきます。また、グラフのy軸のレンジも先ほどの出力をみてここで決めておきます。

# クラスタ数
n = 3

# 元データのプロットレンジ
ymax = 20
ymin = -5

# 標準化データのプロットレンジ
nm_ymax = 2.5
nm_ymin = -2.5

時系列データにおけるk-means法

TimeSeriesKMeansが時系列データにおけるk-means法のクラスです。時系列データの類似度として、まずは単純に同じ時刻の点と点の距離の総和を使ってみましょう。それぞれの時刻で重心を求めるので全ての時系列データの長さが同じでなければ計算できません。ここでの類似度はユークリッド距離という定規で測ったような直線距離(一般的な距離)であるのでeuclideanと指定します。

from tslearn.clustering import TimeSeriesKMeans

km_euclidean = TimeSeriesKMeans(n_clusters=n, random_state=42, metric="euclidean")
labels_euclidean = km_euclidean.fit_predict(nm_data)

# 標準化したデータのプロット
fig, axes = plt.subplots(n, figsize=(8.0, 18.0))
plt.subplots_adjust(hspace=0.5)
for i in range(n):
    ax = axes[i]
    # データのプロット
    for xx in nm_data[labels_euclidean == i]:
        ax.plot(xx.ravel(), "k-", alpha=.2)
    # 重心のプロット
    ax.plot(km_euclidean.cluster_centers_[i].ravel(), "r-")
    # 軸の設定とテキストの表示
    ax.set_xlim(0, 24)
    ax.set_ylim(nm_ymin, nm_ymax)
    datanum = np.count_nonzero(labels_euclidean == i)
    ax.text(0.5, (nm_ymax*0.9+nm_ymin*0.1), f'Cluster {(i + 1)} : n = {datanum}')
    if i == 0:
        ax.set_title("Euclidean k-means (z-score normalization)")


赤線が重心です。1日の終わりの気温が低いものと高いもの、気温が下がっていっているものの3つのグループに分けられました。思ったよりきれいに分けられている印象です。

せっかくなので元のデータをプロットしてみましょう。元のデータでは日平均値にばらつきがあるように見えますが形が似たものがまとめられているように思います。

今後は両方プロットすると冗長になるので以降は標準化したデータのみプロットしますが、元のデータでプロットした方を見ると想像と違っていることがあるかもしれません。こちらでは両方プロットしているので確認してみてください。

Dynamic Time Warping（DTW、動的時間伸縮法）

TimeSeriesKMeansでは類似度にdtwやsoftdtwを使用する事もできます。DTWは時系列の類似度の評価方法としてよく出てくるので先ほどのユークリッド距離の総和での評価との違いを確認しておきます。ユークリッド距離の総和での評価は、波形が似ていても片方が少し遅れた場合(位相がずれた場合)に類似度が低く判断されることがあり、DTWはその問題に対応する手法になっているようです。さらに、時系列データ同士の長さが異なっていても計算することができるという利点もあります。

DTWの計算は動的計画法を使います。二つの時系列x_1,x_2,x_3とy_1,y_2,y_3の間のDTWを例に考えます。x_1とy_2の距離をd_{12}などとして埋めた表が左の図です。次に、下端と左端は初期値として自身の表位置の距離を埋めた表をつくり、残りの位置を左、左下、下の3つのうちから最も小さい値と自身の表位置の距離の和で埋めた表を作ります。そうやって埋めていったときの右表の右上の値がDTWになります。

このブログの動画が分かりやすいと思います。DTWが長さの異なるデータの類似度も評価できることも分かると思います。

DTWを使ったクラスタリングを見てみましょう。重心はそのクラスタ内のそれぞれとのDTWの和が最も小さくなるものになっていると思います。

km_dtw = TimeSeriesKMeans(n_clusters=n, random_state=42, metric="dtw")
labels_dtw = km_dtw.fit_predict(nm_data)

fig, axes = plt.subplots(n, figsize=(8.0, 18.0))
plt.subplots_adjust(hspace=0.5)
for i in range(n):
    ax = axes[i]
    # データのプロット
    for xx in nm_data[labels_dtw == i]:
        ax.plot(xx.ravel(), "k-", alpha=.2)
    # 重心のプロット
    ax.plot(km_dtw.cluster_centers_[i].ravel(), "r-")
    # 軸の設定とテキストの表示
    ax.set_xlim(0, 24)
    ax.set_ylim(nm_ymin, nm_ymax)
    datanum = np.count_nonzero(labels_dtw == i)
    ax.text(0.5, (nm_ymax*0.9+nm_ymin*0.1), f'Cluster {(i + 1)} : n = {datanum}')
    if i == 0:
        ax.set_title("DTW k-means (z-score normalization) (z-score normalization)")

多少クラスタの分け方は変わりましたが、今回のデータは位相のずれのようなものがあまりなさそうなので特に大きな違いはなさそうです。

soft-DTWを使ったクラスタリングを見てみましょう。DTWが微分不可能であるためにDTWを目的関数に使って最適化を行うと計算が不安定になる事があり、そのような時にsoft-DTWを使います。DTWは添え字をそれぞれ選択しているので重心もガタガタになるのですが、soft-DTWは滑らかになります。滑らかさのパラメータgammaは自分で適当に決める必要があります。

km_softdtw = TimeSeriesKMeans(n_clusters=n, random_state=42, metric="softdtw",
                              metric_params={"gamma": .8})
labels_softdtw = km_softdtw.fit_predict(nm_data)

fig, axes = plt.subplots(n, figsize=(8.0, 18.0))
plt.subplots_adjust(hspace=0.5)
for i in range(n):
    ax = axes[i]
    # データのプロット
    for xx in nm_data[labels_softdtw == i]:
        ax.plot(xx.ravel(), "k-", alpha=.2)
    # 重心のプロット
    ax.plot(km_softdtw.cluster_centers_[i].ravel(), "r-")
    # 軸の設定とテキストの表示
    ax.set_xlim(0, 24)
    ax.set_ylim(nm_ymin, nm_ymax)
    datanum = np.count_nonzero(labels_softdtw == i)
    ax.text(0.5, (nm_ymax*0.9+nm_ymin*0.1), f'Cluster {(i + 1)} : n = {datanum}')
    if i == 0:
        ax.set_title("soft-DTW k-means (z-score normalization) (z-score normalization)")

重心が滑らかになっていることが分かります。

カーネル法

KernelKMeansはk-means法で類似度を判定する際にカーネル法を用いるクラスです。カーネル法ではデータを別の高次元空間へ写像したときの、その空間での距離を類似度として使用します。ただし、実際にはデータをいちいち写像せずにカーネル関数を使用して高次元空間での距離(内積)だけを計算しています。これをカーネルトリックと呼びます。

基本的には時系列データ用のカーネル関数（Global Alignment Kernel：GAK）を使用するのだと思うのですが、scikit-learnのpairwise_kernelも使用できるようです。sigmaはsoft-DTWでの滑らかさのようなパラメータですが、autoにしておけば自動で推奨値に設定してくれるようです。

from tslearn.clustering import KernelKMeans

gak_km = KernelKMeans(n_clusters=n, random_state=42, kernel="gak",
                      kernel_params={"sigma": "auto"})
labels_gak = gak_km.fit_predict(nm_data)

fig, axes = plt.subplots(n, figsize=(8.0, 18.0))
plt.subplots_adjust(hspace=0.5)
for i in range(n):
    ax = axes[i]
    # データのプロット
    for xx in nm_data[labels_gak == i]:
        ax.plot(xx.ravel(), "k-", alpha=.2)
    # 軸の設定とテキストの表示
    ax.set_xlim(0, 24)
    ax.set_ylim(nm_ymin, nm_ymax)
    datanum = np.count_nonzero(labels_gak == i)
    ax.text(0.5, (nm_ymax*0.9+nm_ymin*0.1), f'Cluster {(i + 1)} : n = {datanum}')
    if i == 0:
        ax.set_title("gak k-means (z-score normalization)")

クラスタリングの結果はこれまでのものとそれほど変わり映えがありませんが、KernelKMeansは重心を計算してくれないため赤線がありません。

k-Shape法

k-Shape法もk-means法に類似したアルゴリズムで類似度として規格化した相互相関を用いたShape-based distance（SBD）というものを使用しており、SBDの特性に基づいた重心の計算を行っているようです。k-Shape法についてはtslearnのドキュメントにはあまり情報が書いてないので内容を知るには元の論文を読む必要がありますが、解説している日本語ブログもいくつかあるので、まずはそちらに軽く目をとおしても良いかもしれません。

from tslearn.clustering import KShape

ks = KShape(n_clusters=n, random_state=42)
labels_ks = ks.fit_predict(nm_data)

fig, axes = plt.subplots(n, figsize=(8.0, 18.0))
plt.subplots_adjust(hspace=0.5)
for i in range(n):
    ax = axes[i]
    # データのプロット
    for xx in nm_data[labels_ks == i]:
        ax.plot(xx.ravel(), "k-", alpha=.2)
    # 重心のプロット
    ax.plot(ks.cluster_centers_[i].ravel(), "r-")
    # 軸の設定とテキストの表示
    ax.set_xlim(0, 24)
    ax.set_ylim(nm_ymin, nm_ymax)
    datanum = np.count_nonzero(labels_ks == i)
    ax.text(0.5, (nm_ymax*0.9+nm_ymin*0.1), f'Cluster {(i + 1)} : n = {datanum}')
    if i == 0:
        ax.set_title("k-Shape (z-score normalization)")


k-Shape法はこれまでの手法と少し結果が違います。気温が下がっていっているもの分け方に特に違いが出ており、重心の形がこれまでとかなり異なっています。

まとめ

tslearnで使える時系列クラスタリングの使い方を紹介しました。元々のモチベはk-shapeを使ってみたいというものだったのですが、実際やってみると普通のk-meansでもそれなりにクラスタリングができているように思えました。いずれにせよ特徴を把握しつついろいろ試しながら、自分の目的にあった方法を選択したいと思います。
私自身まだまだ勉強中なのでもう少し使ってみてノウハウ的なものをつかみたいです。

Discussion