産業機械の異常状態をローカルレベルモデルを使って検知してみた

#やったこと
以前紹介した[ローカルレベルモデルを用いた異常検知](https://qiita.com/hrkz_szk/items/ea082ca07460ab8b8813)で、実際のデータを用いて解析してみました。
今回用いたデータは異常値を含むデータを多く公開してくれているこちらの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.

#ライブラリ
```python
#基本ライブラリ
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
```

#データ
```python
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

```python
#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')
```

![download1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/246473/cb5fb939-f553-98cd-3028-5b976629c142.png)

#モデリング
全開の記事と違って、今回はPythonのstatsmodelにあるUnobservedComponentsを用いて推定しました。

```python
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|0.066|0.070|
|sigma2.level|0.2173|0.001|233.527|0.000|0.215|0.219|


|||||
|:--:|:--:|:--:|:--:|
|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|

```python
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)
```

![download2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/246473/dbfa3e03-4674-9fbf-bd4e-36b0c23e6b81.png)

#異常検知
```python
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')
```
![download3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/246473/44979f2d-7465-c96f-c00d-d1f4584206c7.png)
3回あるという異常値をうまく推定できています。ただ閾値をどうするかというところについては、主観的判断になってしまうのでしょうか…。統計的アプローチだとしかたないのかも…。