「Pythonによる異常検知」を寄り道写経 ～ 第３章3.5節「時系列データにおける異常検知」機械学習編①スライド窓ｋ近傍法

第３章「時系列データにおける異常検知」

書籍の著者　曽我部東馬 先生、監修　曽我部完 先生

---

この記事は、テキスト「Pythonによる異常検知」第３章「時系列データにおける異常検知」3.5節「時系列データにおける異常検知」3.5.2項「機械学習による時系列データの異常検知」の通称「寄り道写経」を取り扱います。

７回連続で機械学習を利用した時系列データの異常検知を寄り道写経します！
今回はスライド窓ｋ近傍法です。

ではテキストを開いて時系列データの異常検知の旅に出発です🚀

懐中時計のイラスト：「いらすとや」さんより

はじめに

---

テキスト「Pythonによる異常検知」のご紹介

このシリーズは書籍「Pythonによる異常検知」（オーム社、「テキスト」と呼びます）の寄り道写経です。

テキストは、2021年2月に発売され、機械学習等の誤差関数から異常検知を解明して、異常検知に関する実践的なPythonコードを提供する素晴らしい書籍です。
とにかく「誤差関数」と「異常度」の強い結びつきを堪能できる１冊です。

Pythonによる異常検知

引用表記

この記事は、出典に記載の書籍に掲載された文章及びコードを引用し、適宜、掲載文章とコードを改変して書いています。
【出典】
「Pythonによる異常検知」第１版第６刷、著者 曽我部東馬、監修者 曽我部完、オーム社

記事中のイラストは、「かわいいフリー素材集いらすとや」さんのイラストをお借りしています。
ありがとうございます！

かわいいフリー素材集

3.5 時系列データにおける異常検知

---

主に Jupyter Notebook 形式（拡張子 .ipynb）でPythonコードを書きます。

今回の寄り道ポイントです。

• 分類タスクでおなじみの ｋ近傍法（knn）を「スライド窓」単位のデータに適用して異常検知します。
  スライド窓単位に再構成した学習データとテストデータの「ユークリッド距離」を計算して、k個の距離最小値の平均をテストデータの異常度にします。

記事で用いるデータセットはテキスト提供データを引用いたします。

この記事で用いるライブラリをインポートします。

### インポート

# 数値計算、確率計算
import numpy as np 

# 機械学習
from sklearn.neighbors import NearestNeighbors

# 描画
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Meiryo'

スライド窓ｋ近傍法

テキストで初めて知った「スライド窓」。
「窓」は時系列データの移動平均を取る際の「窓」（window）と同じ意味だと思います。
窓幅が３の場合、連続する３期間のベクトルが１つの時間窓になります。
１０期間、窓幅３だと、８個の時間窓が出来ます。
スライド窓はこれらの時間窓をスライドさせる感じ。

詳しい説明はぜひテキストやWeb情報などでご確認ください。
例えば tetsu さんのブログはテキストと同じ分析データを用いて、時間窓のの図示や分析結果のグラフを掲載されています。

k近傍法による時系列データの異常検知 | CCT-recruit

データの読み込み

テキストのデータを引用いたします。
ECGデータと呼ばれる心電図の公開データです。期間の単位はミリ秒です。
テキストによると「時系列異常検知のベンチマーク解析に頻繁に用いられる」そうです。
データ読み込み後、学習データ（100～1999期）とテキストデータ（2001～5999期）に分割します。

### データの読み込みとデータセットの作成 ※テキストのコードを引用、一部改変

# データの読み込み
data = np.loadtxt('./data/discord.txt', delimiter='\t')

# 学習データとテストデータの分割
train_data = data[100:2000, 2]
test_data = data[2001:6000, 2]

# 結果表示
print('train_data.shape:', train_data.shape)
print(train_data[:10])
print('test_data.shape:', test_data.shape)
print(test_data[:10])

【実行結果】

以降の手順はテキストの「STEP」に準拠して記述します。

STEP 1 サンプルを作成する

学習データとテキストそれぞれについて、時間窓サンプルを作成するステップです。
時間窓サンプル作成関数を定義します。

### 時間窓サンプルの作成関数の定義
def make_samples(data, w):    # w: 時間窓幅
    result = [data[i:i+w] for i in range(len(data) - w + 1)]
    return np.array(result)

時間窓幅とｋ近傍法のｋ（近傍点の数）を設定します。

### 設定 ※テキストの設定と同じ
window_width = 50  # 時間窓幅
n_neighbors = 1    # k近傍法のk

では時間窓サンプルを作成します。

### 時間窓サンプルの作成 ※テキストのコードを引用、一部改変

# 学習データとテストデータの時間窓サンプルを作成
train = make_samples(train_data, window_width)
test = make_samples(test_data, window_width)

# 結果表示
print('train.shape:', train.shape)
print(train[:2, :10], '...')
print('test.shape:', test.shape)
print(test[:2, :10], '...')

【実行結果】
時間窓幅＝５０なので、列数が５０になりました。
学習データ、テストデータそれぞれについて、最初の２行・１０列を表示しています。

STEP 2, 3 サンプルとの距離を計算して、異常度を算出する

テストデータの異常度を算出するステップです。

■ ユークリッド距離の計算とｋ個の平均値の計算を理解する

テストデータの時間窓サンプルの１つ１つについて、学習データの時間窓サンプルの$${1851}$$個の全データ点とのユークリッド距離を算出します。
算出したユークリッド距離の中から値の小さい順にｋ個を取り出して平均値を取ったものを異常度にします。
この計算手順に沿って異常度算出関数を２種類作成しました。

１つ目はユークリッド距離の計算式そのものを使いました。

### 異常度算出関数の定義１
#   学習データとテストデータのユークリッド距離を算出してk個の最小値の平均を算出
def calc_dist(train, test, k):
    distances = [sorted(np.sqrt(((train - test[i])**2).sum(axis=1)))
                 for i in range(len(test))]
    dists_mean = np.array(distances)[:, :k].mean(axis=1)
    return dists_mean

２つ目はノルム計算（np.linalg.norm）でユークリッド距離を算出します。

### 異常度算出関数の定義２ ※ノルム計算を利用：np.linalg.norm()
#   学習データとテストデータのユークリッド距離を算出してk個の最小値の平均を算出
def calc_dist2(train, test, k):
    distances = [sorted(np.linalg.norm(train - test[i], axis=1))
                 for i in range(len(test))]
    dists_mean = np.array(distances)[:, :k].mean(axis=1)
    return dists_mean

ではテストデータの異常度を計算しましょう。
１つ目の関数を使ってみます。

### テストデータの異常度の算出１
test_anomaly = calc_dist(train, test, n_neighbors)
print('test_anomaly.shape:', test_anomaly.shape)
print(test_anomaly[:10], '...')

【実行結果】
２行目は異常度の最初の１０個です。

続いて２つ目の関数を使います。

# テストデータの異常度の算出２ ノルム計算を利用
test_anomaly2 = calc_dist2(train, test, n_neighbors)
print('\ntest_anomaly2.shape:', test_anomaly2.shape)
print(test_anomaly2[:10], '...')

【実行結果】
２行目は異常度の最初の１０個は、１つ目の関数で計算した結果と同じになりました。

■ scikit-learn の knn を利用する

ここで、テキストが用いた異常度計算に移ります。
scikit-learn のｋ近傍法：knn で距離を計算する方法です。

まず knn の学習です。学習データの時間窓サンプルを使います。

### NearestNeighborsの学習 ※デフォルトの距離計算方法はminkowski
# ※テキストのコードを引用、一部改変

# k近傍法インスタンス生成
neigh = NearestNeighbors(n_neighbors=n_neighbors)
# 学習データによるフィッティング
neigh.fit(train)

【実行計算】

続いて、距離の計算と平均値＝異常度の算出です。
「.kneighbors(test)」でテストデータの距離を算出しています。

### テストデータの異常度の算出 ※テキストのコードを引用、一部改変（max値で割らない）

# テストデータの距離の算出
d = neigh.kneighbors(test)[0]
# 異常度の算出：最小距離の平均
d = d.mean(axis=1)
# 結果表示
print('テストデータの異常度')
print(d[:10], '...')

【実行結果】
先程の自作関数１・関数２で計算した異常度と同じ値になりました。

STEP 4 異常部位を判断する

テキストを読む限り、異常度の閾値を算出する考えは無さそうです。
おそらくテキストが「異常部位」と呼ぶ「異常度の最大値のサンプル」を探索する感じです。
テキストと同様に可視化で異常部位を目視確認します。

描画関数を定義します。

### 時間窓を利用したスライド窓k近傍法のテスト結果の描画関数の定義

def plot_anomaly(data, anomaly_score, window_width, n_neighbors):
    # 描画領域の設定 gs:gridspecで上下のaxesの高さを5:2にする
    fig = plt.figure(figsize=(10, 5))
    gs = fig.add_gridspec(2, 1, height_ratios=(5, 2))
    ax1 = fig.add_subplot(gs[0, 0])
    ax2 = fig.add_subplot(gs[1, 0])
    # テストデータの観測値の描画
    ax1.plot(data, lw=0.9)
    ax1.set(ylabel="ECG's value")
    ax1.grid(lw=0.5)
    # テストデータの異常度の描画
    ax2.plot(anomaly_score, color='tab:red', lw=0.9)
    ax2.sharex(ax1)
    ax2.set(ylabel='異常度')
    ax2.grid(lw=0.5)
    # 全体修飾
    fig.suptitle(f'スライド窓k近傍法のテスト結果  $width=${window_width}, '
                 f'$k=${n_neighbors}\n上段：観測値、下段：異常度')
    fig.supxlabel('時間[ミリ秒]')
    plt.tight_layout()
    plt.show()

ではテストデータとテストデータの異常度を可視化します。
テキストの図3.44に相当します。

### 異常度の可視化 ★図3.44に相当
plot_anomaly(test_data[1+window_width:], test_anomaly, window_width,
             n_neighbors)

【実行結果】
上段の青い線が観測値、下段の赤い線が異常度です。
観測値は2200～2300ミリ秒あたりに他の時間と異なる部位があります。
同じ時間の異常度は値が大きくなっています。大きな山ができています。
うまく異常検知できているようです！

さまざまなパラメータで検証

テキストの図3.45にならって、時間窓幅やｋ近傍法のｋの値を変えて、異常度に対する影響を確認しましょう。

時間窓サンプル作成～異常度の可視化を一気通貫で実行できるよう、関数化します。

### 実行関数の定義
def sliding_window_knn(train_data, test_data, window_width, n_neighbors):
    # 時間窓サンプルの作成
    train = make_samples(train_data, window_width)
    test = make_samples(test_data, window_width)
    # テストデータの異常度の算出
    test_anomaly = calc_dist(train, test, n_neighbors)
    # 異常度の可視化
    plot_anomaly(test_data[1+window_width:], test_anomaly, window_width,
                 n_neighbors)

■ 時間窓幅＝１、ｋ＝１

###  時間窓幅1, k=1 ★図3.44に相当
# 設定
window_width = 1  # 時間窓幅
n_neighbors = 1   # k近傍法のk
# 実行
sliding_window_knn(train_data, test_data, window_width, n_neighbors)

【実行結果】
異常度の２箇所で値が大きくなっています。

■ 時間窓幅＝５０、ｋ＝１００

###  時間窓幅50, k=100 ★図3.45左に相当
# 設定
window_width = 50  # 時間窓幅
n_neighbors = 100  # k近傍法のk
# 実行
sliding_window_knn(train_data, test_data, window_width, n_neighbors)

【実行結果】
異常度の値が大きい部位は１箇所になりました。
ただし、他の部位の値も大きくなっている感触です。

■ 時間窓幅＝５０、ｋ＝２００

###  時間窓幅50, k=200 ★図3.45中央に相当
# 設定
window_width = 50  # 時間窓幅
n_neighbors = 200  # k近傍法のk
# 実行
sliding_window_knn(train_data, test_data, window_width, n_neighbors)

【実行結果】
全ての異常度の山が大きくなってしまい、異常部位が分かりにくくなりました。

移動平均データの活用

テキストによると、データを「移動平均」でスムージングすると距離のバラツキが抑えられ、異常部位の異常度が周辺の正常部位よりも数桁上がる効果が得られるとのことです。
やってみましょう！

テキストオリジナルコードと実装を変えて移動平均算出関数を定義します。

### 移動平均算出関数の定義
def calc_moving_average(data, width):
    window = np.ones(width) / width 
    return np.convolve(data, window, mode='valid') 

移動平均データを用いた異常度の算出を実行しましょう。
なお、可視化に際して、上段の観測値も移動平均結果をプロットしています。

###  時間窓幅50, k=1, 移動平均 ★図3.45右に相当 ※観測値も移動平均値になっています
# ※テキストと異なるコード

# 設定
n = 20             # 移動平均の窓幅
window_width = 50  # 時間窓幅
n_neighbors = 1    # k近傍法のk

# 移動平均の算出
train_mv = calc_moving_average(train_data, n) 
test_mv = calc_moving_average(test_data, n)

# スライド窓k近傍法による異常検知の実行
sliding_window_knn(train_mv, test_mv, window_width, n_neighbors)

【実行結果】
異常部位の異常部位が際立った感じがいたします。

今回の寄り道写経は以上です。

---

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

---

note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！

１．のんびり統計

統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！

３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！

４．楽しい写経 ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

６．データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

７．Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。