Gunosyデータ分析ブログ

Gunosyで働くデータエンジニアが知見を共有するブログです。

さくっとトレンド抽出: Pythonのstatsmodelsで時系列分析入門

久しぶりの投稿になってしまいましたが、ニュースパス(現在CM放映中!!)開発部の大曽根です。 作業中はGrover Washington Jr のWinelightを聴くと元気が出ます。参加ミュージシャンが素晴らしいですね。

なぜ時系列分析をするのか

数値を非常に重視している弊社では、数値を知るためのツールとしてRedashやChartioおよびSlackへの通知を活用しています。現在の数値を理解する上では、長期のトレンド(指標が下がっているのか、上がっているのか)を知ることが重要です。しかし、日々変化するデータ(特に売上やKPIと言われる指標)は、ばらつきも大きく、変化を適切に捉えることが難しいこともあります。

特にSlackなどへの通知を行っていると、日々の変化に囚われがちです。例えば、弊社ではニュースを扱っていますが、サッカーや競馬のニュースであれば土日に、技術系のニュースであれば平日に偏る傾向があります。日々の数値の浮き沈みに注目していると全体的なトレンドの変化を見逃す可能性もあります。そこで、時系列分析を行うことで変化を適切に捉えたり、予測をある程度正確にすることができます。 今回は周期的なデータとトレンド成分を分割して表現できる季節調整モデルについて紹介します (1週間であれば移動平均などでもある程度のノイズ除去は可能)。

季節調整

季節調整データをざっくり説明すると、時系列のデータを 観測値 = トレンド成分 + 季節成分 + ノイズ成分 で説明するモデルです。

f:id:dr_paradi:20170201104931p:plain

※詳しい説明は北川源四郎 『時系列解析入門』の12章を参考にいただけると幸いです。

弊社の提供するアプリもですが、人間の生活に密着している以上は、人間の生活リズムに影響を受けます。大きく分けて「月要因」、「曜日要因」、「時間要因」がありますが、今回は曜日に注目してサンプルを実装しようかと思います。

実演

本来ならサービスのデータを使いたいのですが、当然使えないので、東京電力さんのデータを使って実装したいと思います。 まずいつものpandas+jupyter+matplotlibで生データを可視化してみます。

data.gunosy.io

import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline

# csvファイルの読み込み
row_data = pd.read_csv('tokyo2016.csv', index_col=0)

df = row_data.reset_index()
df['date'] = pd.to_datetime(df['date'], format='%Y/%m/%d')
df['power'] = df['power'].values.astype(int)
daily_data = df.groupby('date')['power'].sum().reset_index() #元データが1時間ごとなので日毎のデータにする
daily_data = daily_data.set_index(['date'])
daily_data.plot()

すると下記のようなデータが得られます。1年分のデータを可視化する人間の目でもなんとなく傾向は掴めますが、細かい範囲ではデータが波打っており正確なトレンドを把握することは難しいかと思います。

f:id:dr_paradi:20170201013929p:plain

そこで季節調整モデルの出番です。今回は pythonのstatsmodelsというライブラリを用いてデータをトレンド成分 + 季節成分に分けていきます。 StatsModels: Statistics in Python — statsmodels documentation

追記: Pythonのversionは3.5.1、statsmodelsのバージョンは0.8.0rc1を使用しています。

ソースは下記の通りです。

import statsmodels.api as sm
res = sm.tsa.seasonal_decompose(daily_data.values, freq=7)
res.plot()

freq には周期を入力します。今回はdailyのデータということで1週間を周期として設定します。 出力は以下のようになります。グラフの上から 観測データ(生データ)トレンド成分季節成分残差になります。 それぞれの値の関係として

観測データ(生データ) = トレンド成分 + 季節成分 + 残差 が成り立ちます。

f:id:dr_paradi:20170201014033p:plain

季節成分を除いたトレンドを見てみます。

trend = res.trend
trend = pd.DataFrame({'trend': trend, 'date':daily_data.index})
trend['date'] = pd.to_datetime(trend['date'], format='%Y-%m-%d')
trend = trend.set_index(['date'])
trend = trend.plot()

f:id:dr_paradi:20170201010642p:plain

生データと比較すると滑らかになっているのがわかるかと思います。 5月前半のゴールデンウィークなどの大型連休で凹んでいるので、事前に大きな休日などの変化要因を除くことでよりトレンドを正確に把握できるかもしれません。

また、消費電力外気温と適温との差 と考えることもできるので消費電力から気温のトレンドがより分かるようになります。今回は単年でのトレンドを追っていますが、年度ごとに比較することで年ごとの傾向がよりわかりやすくなるでしょう。

次に曜日成分を見てみます。

seasonal_data = res.seasonal
seasonal = pd.DataFrame({'seasonal': seasonal_data, 'date':daily_data.index})
seasonal['date'] = pd.to_datetime(seasonal['date'], format='%Y-%m-%d')
seasonal['weekday'] = seasonal['date'].dt.dayofweek
seasonal[0:7]

f:id:dr_paradi:20170201011320p:plain

weekdayの0は月曜日なので、土日に消費電力は大きく下がることがわかります。おそらく家庭での消費電力よりも工場やオフィスでの消費電力が多いことが影響していそうです。また、月曜日の消費電力が少なく出ていますが、こちらは月曜日に祝日が多いことも影響してそうです。事前にデータを除去したり、平滑化することでトレンド抽出がより正確になる可能性もあります。

以上のように、ただデータを見るだけではなく、季節成分とトレンド成分を分離することで、トレンドの把握が容易になるだけではなく、いろいろな知見を得ることができます。

おまけ: 時間別に見てみる

上記では、1日単位での変動を分析しましたが、時間帯によっても季節変動の差がでるので時間帯別でも同様にトレンド成分を分解してみてみます。

df = row_data.reset_index()
morning_power = df[df['time'] == '8:00'].reset_index()
morning_power['date'] = pd.to_datetime(morning_power['date'], format='%Y/%m/%d')
morning_power['power'] = morning_power['power'].values.astype(int)
morning_data = morning_power.groupby('date')['power'].sum().reset_index()
morning_data = morning_data.set_index(['date'])
res = sm.tsa.seasonal_decompose(morning_data.values, freq=7)
res.plot()

朝8時台のみ抽出したデータ f:id:dr_paradi:20170131013747p:plain

昼の12時台のみ抽出したデータ f:id:dr_paradi:20170131013738p:plain

非常に当たり前の結論になりますが、 トレンド成分を見ると、朝は冬のほうが消費電力が高く、昼は夏の方が消費電力が高いことがわかります。(生データで見てもわかるといえばわかりますが…) 時間帯別に別にトレンドを見ることで当たり前の結論も可視化で明確にすることができます。

まとめ

  • サービスを開発、改善するにあたってのKPIは多くの要因によって変化することが多いので正確に評価するのはなかなか大変です
  • 今回紹介したPythonのstatsmodelを用いて季節変動などの周期的なノイズを除去することでトレンドの把握がしやすくなります
  • 短期的な変動に左右されずにサービスの状態を把握してサービス改善につなげましょう

今後

ノイズ除去後のデータで変化検知とか異常値検知の自動化もしたいですね。