はじめに

先日、Prophetによる時系列データ分析：異常検知編 の記事を書きました。
「時系列変数と目的変数だけ」、「教師なし」という手軽さは他にはないと思っていましたが、Pycaretでもできるようです。

　
以下は、
異常検知にPycaretが適用された記事です。

適用されたデータは、NYCにおける2014年7月から2015年1月までの30分間隔のタクシー利用者数です。
利用者数の推移グラフの挙動だけでは、どのデータが通常外か？はわからないのですが、学習の結果、2014年の終わりから2015年の初めにかけ、膨大な数の異常検知が見られたとあります。（以下参照）

この膨大な異常検知がみられた時期、NYCは激しい吹雪に見舞われたそうです。

Pycaretのモデルは、吹雪により引き起こされたタクシー乗客者行動の珍しいパターンをとてもよく検出したということになりますので、すごいです。

早速、試してみたくなりました。
NYCタクシー利用者のデータは30分間隔であるのに対し、手元にあるデータはある製品の月次間隔の返品データなので、データ間隔は荒く、データ数も少ないですが、一度トライしてみることにしました。

実行条件など

Google colabで実行

手元にある月次返品データ（当記事の末尾にデータ有）で実行

ds y

Dateです返品数です

ds	y
Dateです	返品数です

このライブラリの使い方について

月次の時系列データ（例：「年月」と「出荷数」等）で先の予測や異常検知を行うためのライブラリです。※ 週次や日時、1時間毎のデータもコード変更で扱えます
PyCaret の教師なし異常検知モジュールを利用しています。Isolation Forestによる異常検知と、任意に選択したモデルによる異常検知（※abod / cluster / cof / histogram / knn / lof / svm / pca / mcd / sod / sos をフォームから選択可）を実行します。

前準備

csvデータは以下の形式（例：年月と返品）としてください。

1列目を『年月』、2列目を『目的変数』としてください。 ※データは1,2列のみとしてください。
[注意] csvデータは文字コードを「UTF-8」としてください。

年月返品 →自動変換→ timestamp value

2019/1 4 →自動変換→ 2019-1-1 4

2019/2 12 →自動変換→ 2019-2-1 12

2019/3 6 →自動変換→ 2019-3-1 6

… …　 →自動変換→ … …

2021/12 1 →自動変換→ 2021-12-1 1

年月	返品	→自動変換→	timestamp	value
2019/1	4	→自動変換→	2019-1-1	4
2019/2	12	→自動変換→	2019-2-1	12
2019/3	6	→自動変換→	2019-3-1	6
…	…	→自動変換→	…	…
2021/12	1	→自動変換→	2021-12-1	1

実行手順

「1．初期設定」にて「グラフテーマ」と「時間間隔設定」を選択する（※初期値はMS(月次)としています。週次や日時データの場合、コードを若干変更する必要があります。）
メニューバーの「ランタイム」から「すべてのセルを実行」をクリック。
ライブラリインストール完了後、「3．データ読込み」にて、[ファイル選択]ボタンをクリックし、分析したい文書（csvファイル）を指定する。
他の異常検知モデル（アルゴリズム）を実行する場合は、「7．異常検知のモデル選択と実行」の「Model選択」（Dropdownメニュー）にて任意のモデルを選択する。（※選択変更後は「Model選択」以下のセルを再実行してください。
グラフ上のプロット(赤)が異常検知と認識された点です。

[注意] 一度実行した後にモデルを変更する場合は、すべてのセルを再実行してください。

1．初期設定

#@title グラフテーマ： ggplot2, seaborn, simple_white, plotly, plotly_white, plotly_dark, presentation, xgridoff, ygridoff, gridon, none]

template = 'plotly' #@param ['ggplot2', 'seaborn', 'simple_white', 'plotly', 'plotly_white', 'plotly_dark', 'presentation', 'xgridoff', 'ygridoff', 'gridon', 'none'] {allow-input: true}

#@title 時間間隔設定： AS:年(年初) / A:年(年末) / MS:月(月初) / M:月(月末) / W:週(日曜) / D:日(0時) / H:時(0分) / T:分(0秒) / S:秒

Time_interval = 'MS' #@param ["AS","A","MS","M","W","D","H","T","S"] {allow-input: true}

2．ライブラリインストール&インポート

pip install pycaret

#ライブラリインポート
import plotly.express as px
import plotly.io as pio
pio.renderers.default = "colab"
import pandas as pd

3．データ読込み

#@title csvファイル（UTF-8）を指定してください
from google.colab import files
#print('csvファイル(UTF-8）を指定してください')
uploaded = files.upload()

※ライブラリのインストールとインポートが完了すると、[ ファイル選択 ]がうながされるので、クリックし、ファイルを指定。（キャプチャ画像はファイル読込み後の画面イメージ）

前処理

#@title データフレーム格納とDatetime設定
if len(uploaded.keys()) != 1:
    print("アップロードは１ファイルにのみ限ります")
else:
    target = list(uploaded.keys())[0]

data = pd.read_csv(target)
data.columns = ['timestamp','value']

data['timestamp'] = pd.to_datetime(data['timestamp'])
data.head()

timestamp value

0 2015-01-01 14

1 2015-02-01 11

2 2015-03-01 27

3 2015-04-01 13

4 2015-05-01 24

	timestamp	value
0	2015-01-01	14
1	2015-02-01	11
2	2015-03-01	27
3	2015-04-01	13
4	2015-05-01	24

4．データのグラフ表示

※移動平均は 6が月(MA6)、12か月(MA12)としています。変更する場合はコード変更が必要です。

#@title グラフ表示（移動平均 + rangeslider設定含む）
# create moving-averages
data['MA6'] = data['value'].rolling(6).mean()
data['MA12'] = data['value'].rolling(12).mean()
# plot 
import plotly.express as px
fig = px.line(data, x="timestamp", y=['value', 'MA6', 'MA12'], title='Time Series Anomaly Detection', 
              #template = 'plotly_dark'
              template = template
              )
#fig.update_xaxes(rangeslider_visible=True)
fig.update_xaxes(
    rangeslider_visible=True,
    rangeselector=dict(
        buttons=list([
          dict(count=1, label="1m", step="month", stepmode="backward"),
          dict(count=3, label="3m", step="month", stepmode="backward"),
          dict(count=6, label="6m", step="month", stepmode="backward"),
          dict(count=1, label="YTD", step="year", stepmode="todate"),
          dict(count=1, label="1y", step="year", stepmode="backward"),
          dict(step="all")
        ])
    )
)
fig.show()

グラフ上の確認したい点にカーソルを近づけると時期や値が表示されます。
グラフ下のスライドによって表示範囲を任意に変更できます。
グラフ左上のボタンで、[1m]直近1ヵ月、[3m]直近3ヵ月、[6m]直近6か月、[YTD]今年、[1y]直近1年、[all]全表示の切替が可能です。

5．Pycaret実装

#@title Date要素生成（※月次データ以外の場合はコード内の [注記] 確認のこと）
# drop moving-average columns
data.drop(['MA6', 'MA12'], axis=1, inplace=True)
# set timestamp to index
data.set_index('timestamp', drop=True, inplace=True)

# resample timeseries to (hourly,monthly,yearly…etc) 
data = data.resample(Time_interval).sum()

'''
 [注記]　
 Datetime(timestampカラム)から
 day, day_name(曜日) ,day_of_year ,quarter, week_of_year ,hour ,is_weekday
 を生成する。
 月次データを前提としているため、コメントアウト(#)している項目が多いが、
 週次、日時データを扱う場合は、適宜、コメントアウトを外してください。

'''
# creature features from date
#data['day'] = [i.day for i in data.index]
#data['day_name'] = [i.day_name() for i in data.index]
data['day_of_year'] = [i.dayofyear for i in data.index]
data['quarter'] = [i.quarter for i in data.index]
data['week_of_year'] = [i.weekofyear for i in data.index]
#data['hour'] = [i.hour for i in data.index]
#data['is_weekday'] = [i.isoweekday() for i in data.index]
data.head()

timestamp value day_of_year quarter week_of_year

2015-01-01 14 1 1 1

2015-02-01 11 32 1 5

2015-03-01 27 60 1 9

2015-04-01 13 91 2 14

2015-05-01 24 121 2 18

timestamp	value	day_of_year	quarter	week_of_year
2015-01-01	14	1	1	1
2015-02-01	11	32	1	5
2015-03-01	27	60	1	9
2015-04-01	13	91	2	14
2015-05-01	24	121	2	18

#@title Pycaretセットアップ - セットアップが完了したらEnterキーを押してください
# init setup
from pycaret.anomaly import *
s = setup(data, session_id = 123)

#@title 利用可能モデルの表示
# check list of available models
models()

ID Name Reference

abod Angle-base Outlier Detection pyod.models.abod.ABOD

cluster Clustering-Based Local Outlier pyod.models.cblof.CBLOF

cof Connectivity-Based Local Outlier pyod.models.cof.COF

iforest Isolation Forest pyod.models.iforest.IForest

histogram Histogram-based Outlier Detection pyod.models.hbos.HBOS

knn K-Nearest Neighbors Detector pyod.models.knn.KNN

lof Local Outlier Factor pyod.models.lof.LOF

svm One-class SVM detector pyod.models.ocsvm.OCSVM

pca Principal Component Analysis pyod.models.pca.PCA

mcd Minimum Covariance Determinant pyod.models.mcd.MCD

sod Subspace Outlier Detection pyod.models.sod.SOD

sos Stochastic Outlier Selection pyod.models.sos.SOS

ID	Name	Reference
abod	Angle-base Outlier Detection	pyod.models.abod.ABOD
cluster	Clustering-Based Local Outlier	pyod.models.cblof.CBLOF
cof	Connectivity-Based Local Outlier	pyod.models.cof.COF
iforest	Isolation Forest	pyod.models.iforest.IForest
histogram	Histogram-based Outlier Detection	pyod.models.hbos.HBOS
knn	K-Nearest Neighbors Detector	pyod.models.knn.KNN
lof	Local Outlier Factor	pyod.models.lof.LOF
svm	One-class SVM detector	pyod.models.ocsvm.OCSVM
pca	Principal Component Analysis	pyod.models.pca.PCA
mcd	Minimum Covariance Determinant	pyod.models.mcd.MCD
sod	Subspace Outlier Detection	pyod.models.sod.SOD
sos	Stochastic Outlier Selection	pyod.models.sos.SOS

6．Pycaret実行 - Isoration Forest

#@title Isoration Forestによる学習
# train model
iforest = create_model('iforest', fraction = 0.1)
iforest_results = assign_model(iforest)
iforest_results.head()

timestamp value day_of_year quarter week_of_year Anomaly Anomaly_Score

2015-01-01 14 1 1 1 0 -0.066863

2015-02-01 11 32 1 5 0 -0.125465

2015-03-01 27 60 1 9 0 -0.086232

2015-04-01 13 91 2 14 0 -0.081770

2015-05-01 24 121 2 18 0 -0.107713

timestamp	value	day_of_year	quarter	week_of_year	Anomaly_Score
2015-01-01	14	1	1	1	-0.066863
2015-02-01	11	32	1	5	-0.125465
2015-03-01	27	60	1	9	-0.086232
2015-04-01	13	91	2	14	-0.081770
2015-05-01	24	121	2	18	-0.107713

#@title Isoration Forestによる異常検知データ表示：先頭5データのみ
# check anomalies
iforest_results[iforest_results['Anomaly'] == 1].head()

timestamp value day_of_year quarter week_of_year Anomaly Anomaly_Score

2016-01-01 48 1 1 53 1 0.039575

2016-02-01 57 32 1 5 1 0.003787

2016-04-01 39 92 2 13 1 0.003659

2016-07-01 43 183 3 26 1 0.004296

2017-01-01 39 1 1 52 1 0.022689

timestamp	value	day_of_year	quarter	week_of_year	Anomaly	Anomaly_Score
2016-01-01	48	1	1	53	1	0.039575
2016-02-01	57	32	1	5	1	0.003787
2016-04-01	39	92	2	13	1	0.003659
2016-07-01	43	183	3	26	1	0.004296
2017-01-01	39	1	1	52	1	0.022689

#@title Isoration Forestによる異常検知のグラフ表示
import plotly.graph_objects as go
# plot value on y-axis and date on x-axis
fig = px.line(iforest_results, x=iforest_results.index, y=['value'], title='Unsurpervised Anomaly Detection - Isolation Forest', template = template)
fig.update_xaxes(
    rangeslider_visible=True,
    rangeselector=dict(
        buttons=list([
          dict(count=1, label="1m", step="month", stepmode="backward"),
          dict(count=3, label="3m", step="month", stepmode="backward"),
          dict(count=6, label="6m", step="month", stepmode="backward"),
          dict(count=1, label="YTD", step="year", stepmode="todate"),
          dict(count=1, label="1y", step="year", stepmode="backward"),
          dict(step="all")
        ])
    )
)
# create list of outlier_dates
outlier_dates = iforest_results[iforest_results['Anomaly'] == 1].index
# obtain y value of anomalies to plot
y_values = [iforest_results.loc[i]['value'] for i in outlier_dates]
fig.add_trace(go.Scatter(x=outlier_dates, y=y_values, mode = 'markers', 
                name = 'Anomaly', 
                marker=dict(color='red',size=8)))

fig.show()

グラフ上の確認したい点にカーソルを近づけると時期や値が表示されます。
グラフ下のスライドによって表示範囲を任意に変更できます。
グラフ左上のボタンで、[1m]直近1ヵ月、[3m]直近3ヵ月、[6m]直近6か月、[YTD]今年、[1y]直近1年、[all]全表示の切替が可能です。

#@title TSNEプロット
plot_model(iforest)

7．異常検知のモデル選択と実行

#@title Model選択：abod /cluster /cof /histogram /knn /lof /svm /pca /mcd /sod /sos
Select_Model = 'svm' #@param ["abod","cluster","cof","histogram","knn","lof","svm","pca","mcd","sod","sos"] {allow-input: true}

#@title 選択したModelによる学習（表示する場合は#外す）
# train model
AD_Pycaret = create_model(Select_Model, fraction = 0.1)
AD_Pycaret_results = assign_model(iforest)
#AD_Pycaret_results.head()

#@title 選択モデルによる異常検知データ表示：先頭5データのみ（※表示する場合は#外す）
# check anomalies
#AD_Pycaret_results[iforest_results['Anomaly'] == 1].head()

#@title 選択モデルによる異常検知のグラフ表示
import plotly.graph_objects as go
# plot value on y-axis and date on x-axis
fig = px.line(AD_Pycaret_results, x=AD_Pycaret_results.index, y=['value'], title='Unsurpervised Anomaly Detection - Selected Model', template = template)
fig.update_xaxes(
    rangeslider_visible=True,
    rangeselector=dict(
        buttons=list([
          dict(count=1, label="1m", step="month", stepmode="backward"),
          dict(count=3, label="3m", step="month", stepmode="backward"),
          dict(count=6, label="6m", step="month", stepmode="backward"),
          dict(count=1, label="YTD", step="year", stepmode="todate"),
          dict(count=1, label="1y", step="year", stepmode="backward"),
          dict(step="all")
        ])
    )
)
# create list of outlier_dates
outlier_dates = AD_Pycaret_results[AD_Pycaret_results['Anomaly'] == 1].index
# obtain y value of anomalies to plot
y_values = [AD_Pycaret_results.loc[i]['value'] for i in outlier_dates]
fig.add_trace(go.Scatter(x=outlier_dates, y=y_values, mode = 'markers', 
                name = 'Anomaly', 
                marker=dict(color='red',size=8)))

fig.show()
print("※ Select_Modelは "+"["+Select_Model+"]")

グラフ上の確認したい点にカーソルを近づけると時期や値が表示されます。
グラフ下のスライドによって表示範囲を任意に変更できます。
グラフ左上のボタンで、[1m]直近1ヵ月、[3m]直近3ヵ月、[6m]直近6か月、[YTD]今年、[1y]直近1年、[all]全表示の切替が可能です。

最後に

「時系列変数」をdatetime型に変え、これから取り出した月，日，曜日，週，1年の何日目…等を変数として利用し、学習… 時系列変数と目的変数だけあれば適用できるのはありがたい。
数多くのモデル（アルゴリズム）が選択できるのはいい。Pycaretならでは。
今回適用した月次データは、まぁそこそこか。より時間間隔が短いデータほど効果を発揮するはずだから、（とてもお手軽でもあるので）ぜひ利用したい。

参考サイト

使用したデータ

ds y
2015/1/1 0:00 14
2015/2/1 0:00 11
2015/3/1 0:00 27
2015/4/1 0:00 13
2015/5/1 0:00 24
2015/6/1 0:00 17
2015/7/1 0:00 29
2015/8/1 0:00 27
2015/9/1 0:00 36
2015/10/1 0:00 48
2015/11/1 0:00 40
2015/12/1 0:00 40
2016/1/1 0:00 48
2016/2/1 0:00 57
2016/3/1 0:00 54
2016/4/1 0:00 39
2016/5/1 0:00 41
2016/6/1 0:00 21
2016/7/1 0:00 43
2016/8/1 0:00 32
2016/9/1 0:00 40
2016/10/1 0:00 64
2016/11/1 0:00 36
2016/12/1 0:00 48
2017/1/1 0:00 39
2017/2/1 0:00 42
2017/3/1 0:00 43
2017/4/1 0:00 7
2017/5/1 0:00 69
2017/6/1 0:00 10
2017/7/1 0:00 10
2017/8/1 0:00 7
2017/9/1 0:00 47
2017/10/1 0:00 81
2017/11/1 0:00 17
2017/12/1 0:00 46
2018/1/1 0:00 14
2018/2/1 0:00 14
2018/3/1 0:00 18
2018/4/1 0:00 13
2018/5/1 0:00 23
2018/6/1 0:00 17
2018/7/1 0:00 6
2018/8/1 0:00 6
2018/9/1 0:00 28
2018/10/1 0:00 31
2018/11/1 0:00 12
2018/12/1 0:00 14
2019/1/1 0:00 11
2019/2/1 0:00 10
2019/3/1 0:00 9
2019/4/1 0:00 3
2019/5/1 0:00 15
2019/6/1 0:00 5
2019/7/1 0:00 9
2019/8/1 0:00 6
2019/9/1 0:00 13
2019/10/1 0:00 23
2019/11/1 0:00 31
2019/12/1 0:00 17
2020/1/1 0:00 9
2020/2/1 0:00 19
2020/3/1 0:00 10
2020/4/1 0:00 9
2020/5/1 0:00 3
2020/6/1 0:00 6
2020/7/1 0:00 10
2020/8/1 0:00 5
2020/9/1 0:00 21
2020/10/1 0:00 0
2020/11/1 0:00 8
2020/12/1 0:00 2
2021/1/1 0:00 5
2021/2/1 0:00 3
2021/3/1 0:00 3
2021/4/1 0:00 3
2021/5/1 0:00 1
2021/6/1 0:00 8
2021/7/1 0:00 0
2021/8/1 0:00 4
2021/9/1 0:00 7
2021/10/1 0:00 12
2021/11/1 0:00 11
2021/12/1 0:00 2
2022/1/1 0:00 4
2022/2/1 0:00 4