気象庁の気象データを機械学習して異常検知モデルを構築できるか試してみた
こんにちは、遠藤です。
こちらは、IoTのセンサーデータをクラウドに貯めてなにかするようなシステムを構築するシリーズの記事になります。
前回は、AWS Sage Makerを利用した異常検知についてチュートリアルを試してみました。
今回も、引き続き「異常検知を備えたIoTシステム」を実装することを想定した検証のレポートになります。
過去30年の気象庁の気象データをRandom Cut Forestに学習させてみた
既存の気象データを利用して、「温度センサーなど挙がってきたデータの異常を推論できないか？」という仮説を検証するために、気象庁のサイトからデータをダウンロードし、Random Cut Forestアルゴリズムを使用して学習させました。
データの取得と加工
気象データの取得とCSVに加工
気象データは気象庁の「過去の気象データ・ダウンロード」のページからダウンロードしました。
ダウンロードサイトでは、「気温」、「湿度」、「降水量」といったデータの種類と場所、期間とデータの間隔を選択してダウンロードすることができます。
今回は、日ごとの「平均気温」、「最高気温」、「最低気温」を使います。ただ、一度にダウンロードできるファイルサイズに制限があり10年分だとギリギリサイズオーバーでした。
1994年から2023年まで、5年分ずつに分けてダウンロードして、以下のようなファイル名で保存しました。
CSVを開いて中身を確認してみると、以下のようにヘッダが5行と、データの種類ごとに「品質情報」、「均質番号」列があります。
また、ファイルのエンコードは「Shift-JIS」です。
エンコードの変換と不要なヘッダ行と列の削除はCSVの読み込み時に削除処理をしますので、そのままS3のバケットの適当なパスにアップロードしました。
CSVを1つのデータフレームにまとめる
指定したS3バケットのパスに含まれるCSVファイルをダウンロードして、エンコードShift-JS、ヘッダー行5行を削除して1つのPandsのデータフレームにまとめるコードは以下になります。
まとめた後に、ヘッダーをつけ直して、「品質番号」と「均質番号」を削除しています。
import os
import boto3
import pandas as pd
from io import StringIO
import urllib.request
#バケットへの接続
s3 = boto3.resource('s3')
bucket_name = 'sage-maker-hello-rcf'
bucket = s3.Bucket(bucket_name)
s3_path = 'temperature_tokyo'
# 指定したパス内のファイルを取得
csv_files = [obj.key for obj in bucket.objects.filter(Prefix=s3_path) if obj.key.endswith('.csv')]
# 各CSVファイルをデータフレームとして読み込み、リストに追加
dataframes = []
for path in csv_files:
filename = path.split(os.sep)[1]
bucket.download_file(path, filename)
df2 = pd.read_csv(filename, encoding="shift-jis", skiprows=5) # 5行のヘッダー行を削除して読み込む
dataframes.append(df2)
# 全てのデータフレームを1つに結合
df = pd.concat(dataframes, ignore_index=True)
# 新しいヘッダーラベルを設定
new_headers = ['日時', '平均気温', '平均気温/品質番号', '平均気温/均質番号', '最高気温', '最高気温/品質番号', '最高気温/均質番号', '最低気温', '最低気温/品質番号', '最低気温/均質番号']
df.columns = new_headers
# 品質番号と均質番号列をドロップ
columns_to_drop = [new_headers[2], new_headers[3], new_headers[5], new_headers[6], new_headers[8], new_headers[9]]
df = df.drop(columns=columns_to_drop)
# 結果を表示
df実行結果
欠損値を埋める
欠損値を確認すると、1行だけ最低気温が欠けています。
欠損値の補完方法はいくつかありますが、今回は時系列データなので有効な値を持つ前の行または後ろの行の値で補完をしておきます。
前の行で補完するにはffillを使います。
参考: Pandas: データフレームについて--02:欠損値の取り扱い #Python - Qiita
日付を数値に変換
今回は、日付のデータも特徴量として使用しました。
年間の周期性が考慮されてくれるかもしれない？ので「月日」のみを使います。
また、Randam Cut Forestの特徴量として使えるのは数値のみです。日付は「2024/05/11」のような「yyyy/MM/dd」形式の文字列なので、MMddを数値に変換します。
今回は、月日も特徴量に含めているので「shingling (重複する可能性の検出)」の処理は行いません。
学習させてみる
データの加工が終わったので、Random Cut Forest で学習して、モデルを推論インスタンスにデプロイします。
以下で解説しているコードは、前回解説したチュートリアルのコードをほぼ流用しています。実際にはもっとシンプルに書けるらしいので、どこかで整理したコードを紹介したいところです。
データフレームをエンコードして、S3にアップロード
まずは、加工したデータフレームをエンコードして、S3にアップロードします。
def convert_and_upload_training_data(
ndarray, bucket, prefix, filename='data.pbr'):
import boto3
import os
from sagemaker.amazon.common import RecordSerializer
# convert Numpy array to Protobuf RecordIO format
serializer = RecordSerializer()
buffer = serializer.serialize(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
prefix = 'sagemaker/temperature_tokyo'
s3_train_data = convert_and_upload_training_data(
train_df.to_numpy(),
bucket_name,
prefix)学習の実行
以下のコードを実行してモデルの学習をします。
特徴量は4つあるので、ハイパーパラメータの「feature_dim」は「4」になります。
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_name, prefix),
instance_count=1,
instance_type='ml.c5.xlarge',
sagemaker_session=session)
# ハイパーパラメータの設定
feature_dim = train_df.shape[1]
rcf.set_hyperparameters(
num_samples_per_tree=200,
num_trees=50,
feature_dim=feature_dim)
s3_train_input = sagemaker.inputs.TrainingInput(
s3_train_data,
distribution='ShardedByS3Key',
content_type='application/x-recordio-protobuf')
rcf.fit({'train': s3_train_input})学習は、2分半くらいで終わりました。
参考: RCF Hyperparameters - Amazon SageMaker
推論インスタンスにデプロイ
学習したモデルを推論インスタンスにデプロイします。
from sagemaker.serializers import CSVSerializer
from sagemaker.deserializers import JSONDeserializer
rcf_inference = rcf.deploy(
initial_instance_count=1,
instance_type='ml.c5.xlarge',
)
rcf_inference.content_type = 'text/csv'
rcf_inference.serializer = CSVSerializer()
rcf_inference.deserializer = JSONDeserializer()推論してみる
学習データで推論した結果を確認してみます。
推論結果は、元のデータフレームに戻しています。日付列が残っているので、異常な日付がわかるようにするためです。
results = rcf_inference.predict(train_df.to_numpy())
scores = [datum['score'] for datum in results['scores']]
df['score'] = pd.Series(scores, index=df.index)
score_mean = df.score.mean()
score_std = df.score.std()
score_cutoff = score_mean + 3*score_std
anomalies = df[df['score'] > score_cutoff]
anomaliesこんな感じで、30年のなかで134日異常な気温の日があったと判定しています。割合的に妥当そうですが、実際どうなんでしょうか？
推論結果の検証
2023年11月7日は、11月に最高気温27.5℃は確かに暑すぎです。30年間の11月7日の気温を抽出して箱ひげ図を作ってみると、たしかに2023年の11月7日は突出していると言えそうです。
日時 平均気温 最高気温 最低気温 score DateTime
310 1994/11/7 16.7 20.5 13.6 0.874211 1994-11-07
675 1995/11/7 16.3 19.8 10.9 0.830533 1995-11-07
1041 1996/11/7 12.7 16.0 8.4 0.865156 1996-11-07
1406 1997/11/7 14.2 17.3 11.4 0.826204 1997-11-07
1771 1998/11/7 13.4 15.0 11.7 0.877888 1998-11-07
2137 1999/11/7 15.5 17.6 13.5 0.845668 1999-11-07
2503 2000/11/7 16.8 19.6 13.7 0.870985 2000-11-07
2868 2001/11/7 13.0 15.4 10.9 0.857479 2001-11-07
3233 2002/11/7 9.6 12.3 7.0 1.027629 2002-11-07
3598 2003/11/7 18.5 21.4 16.0 0.982798 2003-11-07
3965 2004/11/7 17.5 21.7 14.5 0.925778 2004-11-07
4330 2005/11/7 18.7 24.1 13.9 1.008837 2005-11-07
4695 2006/11/7 18.5 24.9 12.6 1.007248 2006-11-07
5060 2007/11/7 16.3 19.0 14.8 0.909147 2007-11-07
5426 2008/11/7 17.3 22.3 13.4 0.924352 2008-11-07
5792 2009/11/7 16.5 20.5 13.0 0.867634 2009-11-07
6157 2010/11/7 14.9 18.6 12.6 0.836000 2010-11-07
6522 2011/11/7 18.6 22.7 15.9 1.003222 2011-11-07
6888 2012/11/7 16.5 20.3 12.9 0.866876 2012-11-07
7253 2013/11/7 15.0 17.8 12.8 0.843843 2013-11-07
7619 2014/11/7 17.7 21.1 14.7 0.921245 2014-11-07
7984 2015/11/7 16.5 20.9 13.0 0.875354 2015-11-07
8350 2016/11/7 11.0 14.4 8.2 0.912685 2016-11-07
8715 2017/11/7 15.5 21.2 9.3 0.872713 2017-11-07
9080 2018/11/7 16.6 20.0 13.8 0.881325 2018-11-07
9446 2019/11/7 16.0 22.6 9.8 0.905404 2019-11-07
9812 2020/11/7 15.5 20.6 11.5 0.831879 2020-11-07
10177 2021/11/7 15.7 19.2 12.2 0.832936 2021-11-07
10542 2022/11/7 13.7 18.2 9.6 0.832457 2022-11-07
10907 2023/11/7 22.3 27.5 17.0 1.231373 2023-11-07
逆に、2002年の11月7日が最高気温12.3℃とこちらも外れてそうです。しかし、2002年の11月7日のscoreは「1.027629」でカットオフの「1.1682454…」を下回っているため正常と判定されています。
こちらは、平均気温と、スコア、異常と判定されたスコアをドットとしてプロットしたグラフになりますが、異常と判定される点が冬に集中していて、春夏秋はほとんど検出されません。
どうもこのモデルはイマイチであることが推測されます。こういった場合は、特徴量を変更して学習し直すのがセオリーらしいですが、今回はここまで。そのあたりはまた次回以降に検証してみたいと思います。
おわりに
前回の記事では、「作成したモデルをデプロイしたモデルを利用してLambda上で推論するあたりを試してみたい」と書いていましたが、まだそこまにたどり着けていません。
気象庁のデータを利用したモデルが、今回の目的である「温度センサーなど挙がってきたデータの異常を推論できないか？」に使えるかはまだわかりませんが、もう少し深堀りしてみたいと思います。
ChatGPTに季節性を考慮した特徴量の改善方法をいくつか教えてもらったので、引き続きそれらを試してみたいと思います。ちょっと脱線している気もしますが、次回のレポートも楽しみにしていてください。