Facebookが開発した時系列データ分析ツール「Prophet」を使って東京都新規感染者数の予測をしてみよう

データ分析の醍醐味でもあり目的ともいうべき「予測」にトライする

データ分析のアプローチ方法はいくつもあるのだが、こと「時系列データ」の分析に関しては、どちらかというと「Machine Learningアプローチ」よりもトラディショナルな「統計学アプローチ」の方が精度が高いのではないかという説が、最近気になっている。 おそらくMachine Learningの本業はClassification(つまりクラス分け)であって、ファイナンスデータを予測するには、トラディショナルな統計学アプローチの方が最短ルートなのではないだろうか。

予測するのは東京都新規感染者数

と言いつつ予測するのはファイナンスデータではなく、しかもMLアプローチで攻めていく。
データが正規分布にもとづいているのか、定常性条件を満たしているのかは検証せずに、生きた時系列データを用いて将来を予測する。

hiroshimaeasygogo.hatenadiary.com

東京都の新規感染者数のグラフを描くだけの記事は以前掲載した。
今回はこのコードを生かしていく。 そして、Facebookが開発したお手軽時系列データ分析・予測ツール「Prophet」を使用する。 当然、R派もPython派もプロもアマも使える便利仕様だ。
コードを描く前に是非ともコマンドで pip install fbprophet を叩いて準備しておいてもらいたい。

コードスニペット

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(style="darkgrid")
import pandas as pd
import numpy as np
import statistics
import datetime
import matplotlib.pyplot as plt
import seaborn as sns
plt.style.use('ggplot')


row_data = pd.read_csv('https://stopcovid19.metro.tokyo.lg.jp/data/130001_tokyo_covid19_patients.csv')
by_day = row_data.groupby('公表_年月日')
[print(str(bd[0]) + ': ' + str(len(bd[1])) + '人') for bd in by_day]


dates = []
counts = []
for bd in by_day:
    dates.append(bd[0])
    counts.append(len(bd[1]))
    
# column名を指定する必要があるのに注意
daily_cases = pd.DataFrame({
    'ds':dates,
    'y':counts
})

# 日付をdatetimeに変換する
datetimes = []
for d in daily_cases.ds.values:
    datetimes.append(datetime.datetime.strptime(d, '%Y-%m-%d'))
daily_cases.ds = datetimes

# Prophetを利用した予測
from fbprophet import Prophet
from fbprophet.plot import add_changepoints_to_plot

# MLでお馴染みパラメータチューニング次第で分析の精度が変わる
# トレンドが変化するタイミングを調整することができる
model = Prophet(
    changepoint_range=1,
    changepoint_prior_scale=5,
    daily_seasonality=False,
    weekly_seasonality=True,
    yearly_seasonality=False
)
model.fit(daily_cases)

# 7日間の予測をしてみる
future_df = model.make_future_dataframe(7, freq='D')
forecast_df = model.predict(future_df)

# プロットする
model.plot(forecast_df)
plt.show()

7日分の予測がプロットされている

# パラメータチューニングで調整した「トレンドの転換点」を表す
fig = model.plot(forecast_df)
a = add_changepoints_to_plot(fig.gca(), model, forecast_df)

結構細かめにトレンドの変化を察知してくれたようだ

# 全体のトレンドと曜日によるトレンドに分解したものを図示
model.plot_components(forecast_df)

周知の事実、月曜の報告数が少なくて週末にかけて多くなる傾向

結局明日は何人報告される予測なのか

直近10日の予測データをのぞいてみる。

forecast_df[-10:]

予測値である最終column、"yhat"の値は60.122961人とのこと。

6/30のyhatは60.122961

一方 columns[['yhat_lower', 'yhat_upper']] から予測値の95%信頼区間も示されており、39.442601人〜78.499673人だそう(バッファがありすぎでは...)。 95%信頼区間については、あくまでもその95%の区間内に入る平均値から母集団を推計する為のものなので、少し数値とイメージにギャップがあるが。

いずれにせよ、この自分が作ったモデルが合っているか、答え合わせは随時行う必要がある。 また、今回はファイナンスデータではないデータを用いたが、これがファイナンスデータ(夢はもちろん株価ですね)に応用できるか、日々チェックが必要になりそうである。
先に予想しておくと、先人たちが幾度と試みて挫折してきた歴史をみると、多分株価をはじめとする「ランダムウォーク」するデータの予測は不可能と思われる。
株価の予測は何かしらの多変量分析に定性情報となるイベント情報を組み合わせて、伝統的なテクニカル分析を組み合わせるしかないと思われ、しかも都度都度人の手でパラメータチューニングを行う必要がありそうではあるが...
果たしてどうだろうか。