Amazon SageMakerのRandom Cut Forestをやってみた（異常検知）

Amazon SageMaker 組み込みアルゴリズム入門

#Python

#統計

#機械学習

#ビッグデータ

#Amazon SageMaker

yoshim

2018.05.14

この記事は公開されてから1年以上経過しています。情報が古い可能性がありますので、ご注意ください。

概要

　こんにちは、データインテグレーション部のyoshimです。
　ブログでrcfについてご紹介したのですが、実際にやってみた部分のご紹介ができていなかったので、本エントリーも追加で記述しようと思います。

なお、中身としてはほとんどAWSのブログのそのままです。ちょっとやってみよう、という方のご参考になれば光栄です。

1.Amazon SageMakerの準備

まずは、Amazon SageMakerの準備をしましょう。
やることとしては、「ノートブックインスタンスの作成」と「IAMロールの作成」の2つです。
手順としてはこちらをご参照ください。

なお、2018年5月14日時点では、Amazon SageMakerが利用できるリージョンは下記の4つのみです。
・米国東部 (バージニア北部)
・米国東部 (オハイオ)
・米国西部 (オレゴン)
・EU(アイルランド)

2.チュートリアルの説明

　 AWSのブログでは、NYCのタクシーの利用客数データを利用したチュートリアルを紹介しています。　今回使うデータは「タクシー乗客数を30分ごとに集計したデータ」とのことですので、ちょっと中身を見てみましょう。
　

timestamp  value
2014-11-02 01:00:00  39197
2014-11-02 01:30:00  35212
2014-09-06 23:00:00  30373
2014-09-06 22:30:00  30313
2015-01-01 01:00:00  30236

これを可視化すると下記のような形になります。

横軸が時間軸で、縦軸が乗客数ですね。
横軸が「6000」となっているポイントはやたらと乗客数が多いようですね。チュートリアルではこの日は「ニューヨークシティマラソン」があったらしく乗客数が多かったみたいです。
また、この図だけだとわかりづらいですが、「乗客数が少なかった日」も当然あり、横軸が「8833」の時は大晦日、横軸が「10090」の時は暴風雪だったそうです。

上記のような日を「異常」と見なして分類できるようにすることが今回のチュートリアルの目的です。(こんな感じ)

3.やってみよう

前置きが長くなりましたが、早速実際にやってみようと思います。
とはいえ、ほとんどがチュートリアルに記載されている内容をコピペすれば実行自体はできます。
参考までに、私が実行したコードを残しておきます。

'''1.データを取得する
'''
%matplotlib inline

import pandas
import urllib.request

data_filename = 'nyc_taxi.csv'
data_source = 'https://raw.githubusercontent.com/numenta/NAB/master/data/realKnownCause/nyc_taxi.csv'

urllib.request.urlretrieve(data_source, data_filename)
taxi_data = pandas.read_csv(data_filename, delimiter=',')
taxi_data.plot(title='Taxi Ridership in NYC')

まずは、データを取得しています。

'''1-2.データの中身を見る
'''

print('次元数:{}'.format(taxi_data.shape)) # 次元数
print('\n=======実際のデータ=======\n{}'.format(taxi_data.head(5))) # 最初の10行を目視で確認
print('\n=======実際のデータ=======\n{}'.format(taxi_data.tail(5))) # 最後の10行を目視で確認

# value値（乗客数）
print('\n=======乗客数の基本統計=======')
print('median{}'.format(taxi_data['value'].median()))
print(taxi_data['value'].describe())

# timestamp（時間帯）
print('\n=======時間帯の基本統計=======')
print(taxi_data['timestamp'].describe())

# 時間帯ごとの乗客数
print('\n=======乗客が多い時間帯トップ30=======')
print(taxi_data.sort_values(by='value', ascending=False).head(30))
print('\n=======乗客が少ない時間帯トップ30=======')
print(taxi_data.sort_values(by='value', ascending=True).head(30))

ここはチュートリアルに記載のない部分なのですが、処理を始める前に一旦データの概要を把握するべく、ちょっと確認をしています。
下記のような結果が出て来ます。

2014年7月1日〜2015年1月末までの、各30分ごとのデータで、だいたい各30分ごとに1.7万人弱の乗客数がいるみたいですね。

こちらは時間帯のカラムに重複がないかを確認したものです。

乗客数が多い時は30分で３万人もさばいているんですね...。
1人10分乗客するとしても、1万人くらいタクシードライバーが必要なんですかね...?
それと、時間帯を考えると終電を逃した人がタクシーを利用している、といったイメージでしょうか？
そこはNYCも日本も変わらないですね。

一方、乗客数が少ない時は30分で100人もいないんですね。まあ、時間帯が時間帯ですしね。

'''2.CSV 形式のデータを変換し、そのデータを Amazon S3 バケットへプッシュしています。
'''

def convert_and_upload_training_data(
    ndarray, bucket, prefix, filename='data.pbr'):
    import boto3
    import os
    from sagemaker.amazon.common import numpy_to_record_serializer

    # convert Numpy array to Protobuf RecordIO format
    serializer = numpy_to_record_serializer()
    buffer = serializer(ndarray)

    # upload to S3
    s3_object = os.path.join(prefix, 'train', filename)
    boto3.Session().resource('s3').Bucket(bucket).Object(s3_object).upload_fileobj(buffer)
    s3_path = 's3://{}/{}'.format(bucket, s3_object)
    return s3_path

bucket = 'your-backet' # <-- 利用するバケット名に変更してください
prefix = 'randum_cut_forest' # <-- S3のパス名になります。ここも変更してください。
s3_train_data = convert_and_upload_training_data(
    taxi_data.value.as_matrix().reshape(-1,1),
    bucket,
    prefix)

上記で取得したデータをS3にCSVファイルとしてアップロードします。こんな感じにファイルが生成されるかと思います。

'''3.モデル生成！
'''

import boto3
import sagemaker

# コンテナの一覧
containers = {
    'us-west-2': '174872318107.dkr.ecr.us-west-2.amazonaws.com/randomcutforest:latest',
    'us-east-1': '382416733822.dkr.ecr.us-east-1.amazonaws.com/randomcutforest:latest',
    'us-east-2': '404615174143.dkr.ecr.us-east-2.amazonaws.com/randomcutforest:latest',
    'eu-west-1': '438346466558.dkr.ecr.eu-west-1.amazonaws.com/randomcutforest:latest'}
region_name = boto3.Session().region_name
container = containers[region_name]

session = sagemaker.Session()

# 予測に用いるインスタンスの設定
rcf = sagemaker.estimator.Estimator(
    container,
    sagemaker.get_execution_role(),
    output_path='s3://{}/{}/output'.format(bucket, prefix),
    train_instance_count=1,
    train_instance_type='ml.c5.xlarge',
    sagemaker_session=session)

# robust_random_cut_forestのハイパーパラメータを設定(チュートリアルとは「feature_dim」パラメータの名前が異なる。)
rcf.set_hyperparameters(
    num_samples_per_tree=200, # 各決定木に渡すデータサイズ。1 / num_samples_per_tree = （予想される異常値の割合）となるようにすると良い
    num_trees=50, # 生成する決定木の数。多くてもせいぜい100個くらいまでで良さそう（iForestベースが同じなので）
    feature_dim=1) # 利用する特徴量の割合?ランダムフォレストみたいに各決定木生成に利用する特徴量の割合を制御しているのかと思います

# inputデータの指定
s3_train_input = sagemaker.session.s3_input(
    s3_train_data,
    distribution='ShardedByS3Key',
    content_type='application/x-recordio-protobuf')

# モデル生成
rcf.fit({'train': s3_train_input})

続いてモデルの生成に移ります。
AWSブログ中に記述しているパラメータ名だとエラーが出たので修正しています。
（feature_dimというパラメータ名に修正してください。）気になるパラメータですが、下記の3点を指定しています。

num_samples_per_tree

各決定木に渡すデータサイズです。多ければ多いほど精度が高くなる、という訳ではありません。
元の論文中でも、「1 / num_samples_per_tree = （予想される異常値の割合）」となるように設定することを推奨しています。
ただ、最終的にはこのパラメータはアプリケーションによって最適値が異なるのでチューニングが必要です。

num_trees

作成する決定木の数です。
多い方が計算結果が安定するものの、処理に時間がかかるため、最初はとりあえず100個くらいで進めると良いかもしれません。
このパラメータもアプリケーションによって最適値が異なるのでチューニングが必須ですね。

feature_dim

利用する特徴量の割合？だと思われます。例えば、0.5とかを指定すると各決定木生成時に利用する特徴量は元サブサンプリングが保有する特徴量の半分のみを使う、とかの挙動になるのかもしれません。（すいません、本パラメータについてはまだよくわかってないです） (2018年8月19日追記:データセット内の特徴量の数でした。なので、正の整数を指定します)

'''4.異常スコアを取得する準備
'''

from sagemaker.predictor import csv_serializer, json_deserializer

# 推論エンドポイントを作成
rcf_inference = rcf.deploy(
    initial_instance_count=1,
    instance_type='ml.c5.xlarge',
)

rcf_inference.content_type = 'text/csv'
rcf_inference.serializer = csv_serializer
rcf_inference.deserializer = json_deserializer


'''5.異常スコアを取得
トレーニングセット全体で推論を実行する。
'''

results = rcf_inference.predict(taxi_data.value.as_matrix().reshape(-1,1))
scores = [datum['score'] for datum in results['scores']]
taxi_data['score'] = pandas.Series(scores, index=taxi_data.index)

score_mean = taxi_data.score.mean()
score_std = taxi_data.score.std()

score_cutoff = score_mean + 3*score_std # 平均スコアよりも3標準偏差離れているポイントを異常値としてみなす
anomalies = taxi_data[taxi_data['score'] > score_cutoff]

続いて、先ほど作成したモデルから異常スコアを取得しています。
大事なのは、「異常とみなす閾値」をどうするかですが、今回は「平均スコアよりも3標準偏差離れている」データポイントを異常とみなすように設定しています。ここの部分についても利用するアプリケーション、データセットによって異なってくる部分かと思います。

'''6.予測結果の可視化
異常値のデータポイントをハイライトして可視化する
'''
%matplotlib inline

import matplotlib.pyplot as plt

fig, ax1 = plt.subplots()
ax2 = ax1.twinx()

ax1.plot(taxi_data['value'], alpha=0.8)
ax1.set_ylabel('Taxi Ridership', color='C0')
ax1.tick_params('y', colors='C0')

ax2.plot(taxi_data['score'], color='C1')
ax2.plot(anomalies.index, anomalies.score, 'ko')
ax2.set_ylabel('Anomaly Score', color='C1')
ax2.tick_params('y', colors='C1')

fig.suptitle('Taxi Ridership in NYC')
plt.show()

最後に、ここまでの結果を可視化しています。

青色が乗客数、オレンジ色が異常スコアです。正直これだけだとなんとも判別つかないのですが、なんとなくニューヨークシティマラソンや、大晦日、暴風雪だった日は異常とみなしてくれていそうですね。

4.おまけ

なんとなく気になったのでもう少し結果を確認してみました。

'''6-2.予測結果内容をもう少し詳細に確認する
「1-2.データの中身を見る」で確認した内容と照らし合わせる。
'''

# 異常度の基本統計
print('\n=======異常度の基本統計=======')
print(taxi_data['score'].describe())

# 異常度が高い時間帯
print('\n=======異常度が高い時間帯トップ30=======')
print(taxi_data.sort_values(by='score', ascending=False).head(30))

# 異常度が低い時間帯
print('\n=======異常度が低い時間帯トップ30=======')
print(taxi_data.sort_values(by='score', ascending=True).head(30))

まずは異常スコアの基本統計を見てみました。