#やったこと
以前紹介したローカルレベルモデルを用いた異常検知で、実際のデータを用いて解析してみました。
今回用いたデータは異常値を含むデータを多く公開してくれているこちらのNAB(https://github.com/numenta/NAB)というサイトから持ってきた'machine_temperature_system_failure'というもので、機械に備え付けられているセンサーが読み取った機械の内部温度の推移を示したものです。異常値は3ヶ所あって、最初のものは予定されたシャットダウン、2つ目は検知が難しいもので、それが3つ目の異常値に繋がり、機械止まることになったというものです。
（原文）Temperature sensor data of an internal component of a large, industrial mahcine. The first anomaly is a planned shutdown of the machine. The second anomaly is difficult to detect and directly led to the third anomaly, a catastrophic failure of the machine.
#ライブラリ
#基本ライブラリ
import numpy as np
import pandas as pd
#図形描画ライブラリ
import matplotlib.pyplot as plt
plt.style.use('seaborn-darkgrid')
from matplotlib.pylab import rcParams
rcParams['figure.figsize'] = 15,4
# 統計モデル
import statsmodels.api as sm
#データ
data = pd.read_csv('https://raw.githubusercontent.com/numenta/NAB/master/data/realKnownCause/machine_temperature_system_failure.csv')
data.info()
RangeIndex: 22695 entries, 0 to 22694
Data columns (total 2 columns):
timestamp 22695 non-null object
value 22695 non-null float64
dtypes: float64(1), object(1)
memory usage: 354.7+ KB
#change the type of timestamp for plot
data['timestamp'] = pd.to_datetime(data['timestamp'])
#change fahrenheit to c
data['value'] = (data['value'] - 32) * 5/9
#plot
rcParams['figure.figsize'] = 16,4
data.plot(x='timestamp', y='value')
#モデリング
全開の記事と違って、今回はPythonのstatsmodelにあるUnobservedComponentsを用いて推定しました。
mod_local_level = sm.tsa.UnobservedComponents(data.value, 'local level')
res_local_level = mod_local_level.fit()
print(res_local_level.summary())
| sales_ad | dummy_ad | |
|---|---|---|
| 0 | 22.834287 | 0.0 |
| 1 | 26.212849 | 0.0 |
| 2 | 25.638334 | 0.0 |
| 3 | 26.333011 | 0.0 |
| 4 | 25.669880 | 0.0 |
Unobserved Components Results
| Dep. Variable: | value | No. Observations: | 22695 |
|---|---|---|---|
| Model: | local level | Log Likelihood | -19945.532 |
| Date: | Sat, 11 May 2019 | AIC | 39895.065 |
| Time: | 13:39:42 | BIC | 39911.124 |
| Sample: | 0 | HQIC | 39900.287 |
| Covariance Type: | opg |
| coef | std err | z | P>z | [0.025 0.975] |
|---|---|---|---|---|
| sigma2.irregular | 0.0679 | 0.001 | 55.691 | 0.000 |
| sigma2.level | 0.2173 | 0.001 | 233.527 | 0.000 |
| Ljung-Box (Q): | 1673.00 | Jarque-Bera (JB): | 282754.21 |
| Prob(Q): | 0.00 | Prob(JB): | 0.00 |
| Heteroskedasticity (H): | 0.94 | Skew: | 1.38 |
| Prob(H) (two-sided): | 0.01 | Kurtosis: | 20.07 |
rcParams['figure.figsize'] = 15,5
plt.plot(data.value.iloc[1:], label='Observed')
plt.plot(res_local_level.fittedvalues[1:], label='Predicted')
plt.legend(loc='lower right', borderaxespad=0, fontsize=15)
#異常検知
sigma = res_local_level.params[0] + res_local_level.params[1]
anomaly_detection = (data.value - res_local_level.fittedvalues[1:])*sigma
plt.plot(anomaly_detection)
plt.title('anomaly detection')
3回あるという異常値をうまく推定できています。ただ閾値をどうするかというところについては、主観的判断になってしまうのでしょうか…。統計的アプローチだとしかたないのかも…。