はじめに

研究で時系列データの異常検知に関する研究を行っています。そのため、機械学習による異常検知のアルゴリズムについて普段から文献調査などしているわけですが、機械学習を使わずとも異常検知できるという記事を見つけました。

大変興味深かったので、記事の内容に少しだけアレンジを加えて実装してみたいと思います。

対象データ

自分で異常検知用のデータを作成しても良いのですが、それだとどうしても作業的になってしまい、面白くないのでこちらの心電図のデータを使わせていただきました。
ダウンロードしたデータをそのままグラフのさせると以下のようになります。ただし、そのまま出力させるとデータが多過ぎて潰れたグラフになって見にくいので、データの頭から32秒付近までの間のデータを抽出しています。

f:id:Noleff:20200707152820p:plain

異常検知

見つけたい異常

上のグラフからわかるように、17秒付近で他とは振動をしていることがわかると思います。今回はこの異常を検知したいと思います。

特徴量

1. 単純移動平均

一つ目の特徴量は移動平均です。移動平均とはある区間を決めてその区間をずらしていき、各区間内で平均値を求めていく計算になります。
平均は一般的にあまり有効な特徴量となり得ませんが、今回は二つ目の特徴量を際立たせるためにあえて平均を採用しています。（といっても今回の場合この特徴量を使っていないとうまくいってないところもありますが）
以下に移動平均のグラフを載せます。移動平均なため、データの頭は0で埋めていることに注意です。

f:id:Noleff:20200707143036p:plain

2. ローパスフィルタとの差分

フーリエ変換

ローパスフィルタを説明するために、まずフーリエ変換について簡単に触れます。

フーリエ変換とは実変数の複素または実数値関数を別の同種の関数に写す変換である。変換後の関数はもとの関数に含まれる周波数を記述し、しばしばもとの関数の周波数領域表現 (frequency domain representation) と呼ばれる。実質的に、フーリエ変換は関数を振動関数に分解する。（Wikipediaより） https://ja.wikipedia.org/wiki/%E3%83%95%E3%83%BC%E3%83%AA%E3%82%A8%E5%A4%89%E6%8F%9B

言葉にするとわかりにくいですが、同じくWikipediaの載ってあるGIF画像を見ると非常にわかりやすいです。
複雑な振動波形でも、複数の振動波形が合成されてできています。この一つ一つの振動波形を取り出し、その振動波形からどの周期にどのぐらいの振れ幅を持っているかわかります。この振動波形から周波数成分に分解する変換のことをフーリエ変換といいます。
詳しくはこちらのサイトが大変わかりやすかったです。

ローパスフィルタ

続いてローパスフィルタの説明をします。フーリエ変換をすることで、振動波形を周波数成分に変換しました。今回の心電図をデータをフーリエ変換すると以下のようになります。ただし、0から10Hzの範囲のみ出力させています。なお、縦軸は各周期数成分の大きさを表す振幅スペクトルになります。

f:id:Noleff:20200707140545p:plain

大きな振幅スペクトルが二つ、中くらいの振幅スペクトルが一つ、小さな振幅スペクトルが三つできているのがわかると思います。
ローパスフィルタはこれらの周期数の内、Lowな振幅スペクトルの周期数のみ抽出します。つまり、振幅スペクトルが小さい高周波を除去しているイメージです。ローパスフィルタ後の心電図とそれをフーリエ変換したしたものを以下に載せます。

f:id:Noleff:20200707140640p:plain

f:id:Noleff:20200707140709p:plain

今回ローパスフィルタは2.5Hzよりも大きな値を除去したので、二つ目の振幅スペクトルがやや鋭角になってますが、今回はこれでいきたいと思います。

差分の計算

長くなりましたが、ここまでが前置きです。
2つ目の特徴量はローパスフィルタとの差分でした。つまり、ローパスフィルタと元のデータの差分をとり、それを異常度として算出します。上のグラフでいうと、橙がローパスフィルタ後の心電図、青が元のデータの心電図なため、橙から青を引いた絶対値を異常度とします。なお、以下のグラフは0から1の範囲に値が収まるように正規化してます。

f:id:Noleff:20200707135208p:plain

思ったよりきれいに異常度を算出できませんでした（汗）。
ですが、心電図に異常が起きている17秒付近で異常度が大きくなっていることは見てわかると思います。

プログラム

環境

OS：windows10
python：3.7.6
numpy：1.18.1
pandas：1.0.1
scikit-learn：0.22.1

コード

上の描画のコードなど一部省略します。ほとんど参考から流用したので、気になる方はそちらをご覧ください。

生データ

最初に読み込むecg.csvはこんな感じです。timeが200スタートだったり、signal1とsignal2の違いは詳しく調べてないのでわかりません。今回はsignal2の方を使用しました。データはtimeからわかるように0.004秒ごとに記録されています。

f:id:Noleff:20200707141231p:plain

関数

# データ取得
def read_df(path):
    df = pd.read_csv(path)
    return df

# Timeから時間の差を追加
def differ_time(df):
    # Timeから時間の差を割り当てる
    df['dif_sec'] = df['time'].diff().fillna(0)
    df['cum_sec'] = df['dif_sec'].cumsum()
    return df

# 正規化
def minxmax_scaler(df):
    scaler = MinMaxScaler()
    scaler.fit(df)
    # 正規化したデータを新規のデータフレームに
    df_mc = pd.DataFrame(scaler.transform(df), columns=df.columns) 
    # 正規化したデータをリストに
    # mc_list = scaler.fit_transform(df)
    return df_mc

以下べた書き

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import MinMaxScaler

# データ読み込み
df = read_df('./ecg.csv')
dif_df = differ_time(df) # 秒数を計算 
test_df = dif_df.iloc[0:8192, :] # 使うデータだけ抽出

N = len(test_df) # FFTのサンプル数 
dt = 0.004 # サンプリング周波数 
freq = np.fft.fftfreq(N, d=dt) # 周波数

# 高速フーリエ変換(FFT)
F = np.fft.fft(test_df['signal2'].values) 
F_abs = np.abs(F/(N/2))

# ローパスフィルタ
threshold_period = 0.4
threshold_freq = 1 / threshold_period
F_lowpass = np.where(abs(freq) > threshold_freq, 0, F)
lowpass = np.fft.ifftn(F_lowpass).real

# 特徴量
sma = test_df['signal2'].rolling(25).mean().fillna(0) # 0.1秒ごとの移動平均
diff = abs(lowpass - test_df['signal2']) # ローパスフィルタとの差分

# 結果をまとめるデータフレーム
result_df = pd.DataFrame(diff, columns=['sma', 'diff', 'anomaly', 'label'])
result_df['sma'] = sma
result_df['diff'] = diff

 # 異常度を正規化
result_df = minxmax_scaler(result_df)

# 移動平均とローパスフィルタとの差分を乗算
result_df['anomaly'] = result_df['sma'] * result_df['diff'] 

# ラベル振り
result_df.loc[result_df['sma'] == 0, 'label'] = 2
result_df.loc[result_df['anomaly'] > 0.6, 'label'] = 1
result_df['label'] = result_df['label'].fillna(0)

これで、結果をまとめたresult_dfができました。

結果

移動平均とローパスフィルタを乗算した結果を示します。

plt.figure(figsize=(16,4))
plt.plot(test_df['cum_sec'], minxmax_scaler(test_df)['signal2'], color='blue', label='row')
plt.plot(test_df['cum_sec'], result_df['anomaly'], color='red', label='anomaly')
plt.legend()
plt.xlabel('sec')
plt.ylabel('anomaly')
plt.show()

f:id:Noleff:20200707144458p:plain

異常度が0から1の範囲に収まるようにしているので、閾値を設定すれば異常検知できます。今回は異常度が0.6より大きい箇所を異常としています。
グラフを見ると、心電図の17秒付近のみが0.6を越えていることがわかります。

今度は散布図として出力させます。

c_list = ['blue', 'red', 'green']

for i, c in enumerate(c_list):
    plt.scatter(result_df.loc[result_df['label']==i, 'diff'], result_df.loc[result_df['label']==i, 'sma'], color=c, label=i)

plt.xlabel('diff')
plt.ylabel('sma')
plt.legend()
plt.show()