Prophetは,Facebookが公開している予測パッケージである.
from IPython.display import Image, YouTubeVideo
YouTubeVideo("23VJTLk1twM")

https://mikiokubo.github.io/analytics/15forecast/

インストール

Macの場合

conda install -c conda-forge fbprophet

Google Colabの場合

!pip install fbprophet

Windowsにインストールする際には,C言語のコンパイラが必要

https://github.com/facebook/prophet

参照

# インストールされていない場合には、以下を生かす。
#!pip install -U fbprophet
#!pip install -U vega_datasets

需要予測

筆者は,以前は最適化の仕事は引き受けても予測の仕事は引き受けないことにしていた.予測は現場で長年働いている人の経験を加味して行うべきものであり,門外漢がいくらテクニックを駆使してもそれを超えることは難しいかと考えていたからだ.

しかし最近になって,需要予測は当たらないといけないという考えを捨てて,誤差の管理を行うことだと割り切って考えることにして,予測に積極的に関与するようになった.多くの最適化モデルは,予測がある程度合っているという前提で構築される.予測をいいかげんにされると,最適化モデル自体が役に立たないものになる.ゴミを入れればゴミが出てくるからだ.

誤差を管理することによって,需要予測を「点」で行うのではなく,「範囲」で行うことが可能になる.また,需要の分布も特定できるようになる.範囲内での最適化はロバスト最適化,確率分布を仮定した最適化は確率最適化という枠組みで解決可能になる.

以下では,サプライ・チェイン最適化に関連した需要予測の基本と重要なモデルについて解説する.

予測の公理

最近,サプライ・チェインの現場において,予測に関する多くの誤りが浸透していることに気づいた.ここでは,このような誤用を減らして正しい予測手法を適用するための,予測の公理について述べる.

  1. 予測は予測のためにあらず

しばしば,予測を目的として仕事をしている人たちを見かける.予測は,ほかの重要な意思決定行うための基礎となる手段であり,予測そのものを目的としてはならない.実際には,予測よりもその誤差を評価することの方が重要である.誤差が増えているのか,減っているのか,その理由は何かを考えることが,需要予測の真の目的なのである.

たとえば,小売りの現場で需要予測を行うことは,在庫費用の削減や品切れ損出の回避などを目的としたものであり,予測の精度だけを問題にするのではなく,どの程度外れているのかという誤差の管理と,外れた場合の影響や,緊急発注などの回避手段とあわせて考える必要がある.

また,予測するだけでなく,なぜそのような値になったのかを究明することも重要である.需要が0という日が続いた場合には,それが,本当に需要がなかったのか,それとも在庫がないために売れなかったのか,陳列場所が悪かったために売れなかったのか,などをあらわせて原因を分析する必要がある.

多くのメーカーでは,需要予測をするためだけの部署を設けているが,これもナンセンスなのでできるだけ早くほかの部署と連携をとるように改めるべきである.需要は当てるものではなくコントロールするものであり,需要予測をあてることだけを目標としている部署は,廃止すべきである.

  1. 予測は外れるもの

しばしば予測が当たったとか外れたとかいう言葉を現場の人から聞くが,経営はギャンブルではないので的中というのはありえない.この人なら当たるとかいうのは迷信であり,たまたま当たったときに声を張り上げて宣伝しているか,誤差が大きいにもかかわらず当たったと宣伝しているかの何れかである.サプライ・チェインからはちょっと外れるが,地震の予測(予知ともいう)も似たようなものであり,日本中のどこかで地震が発生すると予測し,そのうち1つがあたると的中と宣伝していたりする.ましてや,株や競馬の予想的中などの宣伝は,たまたま当たったときの結果だけを掲載し,外れたときのものを消去して,予想的中の証拠として提出していたりする.いまだに,こんな宣伝にだまされる人がいるのかと感心するが,社内で需要予測が的中する人がいるという会社も似たような詐欺に遭っていると言える.

重要なことは,どの程度外れたのかを時系列的に管理することである.誤差が増えている際には,その原因を追及し,予測手法を改善するなり,在庫を増やして品切れを回避するなりの行動をとるのが正しい方法である.地震学者も予測だけでなく,学問として成立させるためには,誤差を時系列的に管理し,それを公表すべきであろう.

  1. 集約すれば精度があがる

往々にして,個々のものの予測は難しいが,それをまとめたものの予測は容易になる.たとえば,特定の場所で特定のマグニチュードの地震が明日発生することを予測するのは難しいが,日本のどこかでマグニチュード4以上の地震が来年発生することは容易に予測できる.前者はほぼ0%であるが,後者はほぼ100%の精度で予測できる.これが集約の力である.

サプライ・チェインでも同様であり,商品を集約して商品群にすれば精度があがり,個々の店舗での売り上げでなく,地域内のすべての店舗に集約すれば精度があがる.時間軸でも同様であり,1時間以内の需要を予測精度は,日,週,月,年単位と集約しているにしたがってあがっていく.

予測精度だけを議論するのであれば,どんどん集約した方が得であるが,もちろんどんどん役に立たなくなる.前述したように,重要なことは,その予測を何に使うかであり,使用法にあわせて「適切に」集約を行うことである.

また,製品設計や在庫地点を考慮することによって,物理的に集約を行うこともできる.たとえば,製品のモジュール化や遅延差別化やリスク共同管理がこれに相当する.これらはモダンなサプライ・チェインの基本戦略であり,そのすべてがこの単純な公理(集約した方が予測精度があがる)に基づくものであることは興味深い.

  1. 目先の需要は当たりやすく,遠い未来の需要は予測しにくい

明日の天気はある程度予測できるが,1年後の天気は予測しにくい.不確実性は時間が経過するにしたがって増大するからである.需要予測も同様である.たとえば,カップラーメンなどは一部の定番品を除いて,来年店頭に並んでいるかどうかも怪しい.

近い未来の予測は,短期のオペレーショナルな意思決定に用いるので,比較的正確性が必要であるが,遠い未来の予測は長期のストラテジックな意思決定に用いるので,おおよそで構わない.さらに,長期の意思決定の際には,データは集約して行われるので,予測精度も上がる. 要は,目的のために適切な精度で管理できるように,時間軸を集約することが推奨される.たとえば,近い未来は予測精度がよいので,集約をせずに日単位で予測し,未来に行くに従って時間軸の集約を行い,週単位,月単位,年単位としていく訳である.このテクニックはテレスコーピングと呼ばれ,時間軸を含んだ実際問題を解くときによく用いられる.

予測手法と使い分け

古典的な予測手法は statsmodelsモジュールに入っている。種々の指数平滑法の拡張や、ARIMAなどは簡単にできる。 いずれも高速なので、補助的なデータがない時系列データに対して、簡易的な予測をしたい場合には便利だが、予測精度は期待すべきではない。

機械(深層)学習を用いた予測は、 scikit-kearnやfastaiのTabularモデルを用いることによって容易にできる。 補助的なデータが豊富な場合には、これらの手法が推奨される。単なる時系列データの場合には、深層学習のLSTMを使うことも考えられるが、精度は期待すべきでない。

補助的なデータが少ない時系列データの場合には、prophetを用いたBayes推論が推奨される。Bayes推論だと、予測を点でするのではなく、不確実性の範囲まで得られるので、 サプライ・チェインのモデルと相性が良い。例えば、在庫モデルにおいては、不確実性の情報が不可欠だからだ。

from IPython.display import Image
Image("../figure/how_to_select.PNG",width=1000,height=500)
Image("../figure/bayes1.PNG",width=1000,height=500)
Image("../figure/bayes2.PNG",width=1000,height=500)
Image("../figure/bayes3.PNG",width=1000,height=500)
Image("../figure/mcmc.PNG",width=1000,height=500)

諸パッケージのインポート

import pandas as pd
from fbprophet import Prophet
from vega_datasets import data
import plotly.express as px
import fbprophet.plot as fp
#import plotly.io as pio
#pio.renderers.default = "colab"

基本的な使い方

ProphetをPythonから呼び出して使う方法は,機械学習パッケージscikit-learnと同じである。

  1. Prophetクラスのインスタンスmodelを生成
  2. fitメソッドで学習(引数はデータフレーム)
  3. predictメソッドで予測(引数は予測したい期間を含んだデータフレーム)

例としてアメリカンフットボールプレーヤのPayton ManningのWikiアクセス数のデータを用いる。

df = pd.read_csv('http://logopt.com/data/peyton_manning.csv')
df.head()
ds y
0 2007-12-10 9.590761
1 2007-12-11 8.519590
2 2007-12-12 8.183677
3 2007-12-13 8.072467
4 2007-12-14 7.893572

Prophetモデルのインスタンスを生成し,fitメソッドで学習(パラメータの最適化)を行う.fitメソッドに渡すのは,上で作成したデータフレームである.このとき、ds(datestamp)列に日付(時刻)を、y列に予測したい数値を入れておく必要がある。 (この例題では,あらかじめそのように変更されている.)

#model = Prophet()
#model.fit(df)

make_future_dataframeメソッドで未来の時刻を表すデータフレームを生成する。既定値では、予測で用いた過去の時刻も含む。 引数は予測をしたい期間数periodsであり,ここでは、1年後(365日分)まで予測することにする。

#future = model.make_future_dataframe(periods=365)
#future.tail()

predict メソッドに予測したい時刻を含んだデータフレームfuture を渡すと、予測値を入れたデータフレームforecastを返す。このデータフレームは、予測値yhatの他に、予測の幅などの情報を含んだの列を含む。以下では,予測値yhatの他に,予測の上限と下限(yhat_lowerとyhat_upper)を表示している.

#forecast = model.predict(future)
#forecast[['ds', 'yhat', 'yhat_lower', 'yhat_upper']].tail()

matplotlibを用いた描画は,plotメソッドで行う.

#model.plot(forecast);

Prophetにおける予測は一般化加法モデルを用いて行われる.これは,傾向変動,季節変動,イベント情報などの様々な因子の和として予測を行う方法である.

一般化加法モデル

$$ y_t =g_t + s_t + h_t + \epsilon_t $$
  • $y_t$ : 予測値
  • $g_t$ : 傾向変動(trend);傾向変化点ありの線形もしくはロジスティック曲線
  • $s_t$ : 季節変動;年次,週次,日時の季節変動をsin,cosの組み合わせ(フーリエ級数)で表現
  • $h_t$ : 休日などのイベント項
  • $\epsilon_t$ : 誤差項

因子ごとに予測値の描画を行うには,plot_componentsメソッドを用いる.既定では,以下のように,上から順に傾向変動,週次の季節変動,年次の季節変動が描画される.また,傾向変動の図(一番上)には,予測の誤差範囲が示される.季節変動の誤差範囲を得る方法については,後述する.

#model.plot_components(forecast);

対話形式に,拡大縮小や範囲指定ができる動的な図も,Plotlyライブラリを用いて得ることができる.

#fig = fp.plot_plotly(model, forecast) 
#fig.show()

CO2排出量のデータ

データライブラリから読み込み,Plotly Expressで描画する.

co2 = data.co2_concentration()
co2.head()
Date CO2 adjusted CO2
0 1958-03-01 315.70 314.44
1 1958-04-01 317.46 315.16
2 1958-05-01 317.51 314.71
3 1958-07-01 315.86 315.19
4 1958-08-01 314.93 316.19
#px.line(co2,x="Date",y="CO2")

列名の変更には,データフレームのrenameメソッドを用いる.引数はcolumnsで,元の列名をキーとし,変更後の列名を値とした辞書を与える.また,元のデータフレームに上書きするために,inplace引数をTrueに設定しておく.

co2.rename(columns={"Date":"ds","CO2":"y"},inplace=True)
co2.head()
ds y adjusted CO2
0 1958-03-01 315.70 314.44
1 1958-04-01 317.46 315.16
2 1958-05-01 317.51 314.71
3 1958-07-01 315.86 315.19
4 1958-08-01 314.93 316.19

make_future_dataframeメソッドで未来の時刻を表すデータフレームを生成する。既定値では、(予測で用いた)過去の時刻も含む。 ここでは、200ヶ月先まで予測することにする。

そのために,引数periodsを200に,頻度を表す引数freqをMonthを表す'M'に設定しておく

predict メソッドに予測したい時刻を含んだデータフレームfuture を渡すと、予測値を入れたデータフレームforecastを返す。このデータフレームは、予測値yhatの他に、予測の幅などの情報を含んだの列を含む。

最後にplotメソッドで表示する.

# model = Prophet()
# model.fit(co2)
# future = model.make_future_dataframe(periods=200, freq='M')
# forecast = model.predict(future)
# model.plot(forecast);

予測は一般化加法モデルを用いて行われる.

これは,傾向変動,季節変動,イベント情報などの様々な因子の和として予測を行う方法である.

上に表示されているように,週次と日時の季節変動は無視され,年次の季節変動のみ考慮して予測している.

因子ごとに予測値の描画を行うには,plot_componentsメソッドを用いる.既定では,以下のように,上から順に傾向変動,週次の季節変動,年次の季節変動が描画される.また,傾向変動の図(一番上)には,予測の誤差範囲が示される.季節変動の誤差範囲を得る方法については,後述する.

#model.plot_components(forecast);

Plotlyで描画

一部を拡大,期の選択などが可能になる.

#fig = fp.plot_plotly(model, forecast)
#fig.show()

航空機乗客数のデータ

季節変動を加法的でなく乗法的にする.

passengers = pd.read_csv("http://logopt.com/data/AirPassengers.csv")
passengers.head()
Month #Passengers
0 1949-01 112
1 1949-02 118
2 1949-03 132
3 1949-04 129
4 1949-05 121
#px.line(passengers,x="Month",y="#Passengers")
passengers.rename(inplace=True,columns={"Month":"ds","#Passengers":"y"})
passengers.head()
ds y
0 1949-01 112
1 1949-02 118
2 1949-03 132
3 1949-04 129
4 1949-05 121

Prophetの規定値では季節変動は加法的モデルであるが、問題によっては乗法的季節変動の方が良い場合もある。例として、航空機の乗客数を予測してみよう。最初に既定値の加法的季節変動モデルで予測し,次いで乗法的モデルで予測する.

そのためには,モデルのseasonality_mode'multiplicative'に設定する.

なお,以下のデータは月次のデータであるので,make_future_dataframeのfreq引数を'M'(Month)に設定する.

# model = Prophet().fit(passengers)
# future = model.make_future_dataframe(periods=20, freq='M')
# forecast = model.predict(future)
# model.plot(forecast);
# model = Prophet(seasonality_mode='multiplicative').fit(passengers)
# future = model.make_future_dataframe(periods=20, freq='M')
# forecast = model.predict(future)
# model.plot(forecast);
#model.plot_components(forecast);

問題

以下の,小売りの需要データを描画し,予測を行え. ただし,モデルは乗法的季節変動で,月次で予測せよ.

retail = pd.read_csv('http://logopt.com/data/retail_sales.csv')
retail.head()
ds y
0 1992-01-01 146376
1 1992-02-01 147079
2 1992-03-01 159336
3 1992-04-01 163669
4 1992-05-01 170068

1時間ごとの気温データ

date列に日付と時間(1時間ごと)が,temp列に気温データが入っている.

日付未満のデータ

日別でないデータも扱うことができる。データ形式は、日付を表すYYYY-MM-DDの後に時刻を表すHH:MM:SSが追加されている。 未来の時刻を表すデータフレームは、make_future_dataframeメソッドで生成するが、このとき引数freqで時間の刻みを指定する。 ここでは1時間を表す'H'を指定する。

climate = data.seattle_temps()
climate.head()
date temp
0 2010-01-01 00:00:00 39.4
1 2010-01-01 01:00:00 39.2
2 2010-01-01 02:00:00 39.0
3 2010-01-01 03:00:00 38.9
4 2010-01-01 04:00:00 38.8
climate["Date"] = pd.to_datetime(climate.date)
climate.rename(columns={"Date":"ds","temp":"y"},inplace=True)
#px.line(climate,x="ds",y="y")
# model = Prophet().fit(climate)
# future = model.make_future_dataframe(periods=200, freq='H')
# forecast = model.predict(future)
# model.plot(forecast);
INFO:fbprophet:Disabling yearly seasonality. Run prophet with yearly_seasonality=True to override this.
#model.plot_components(forecast);

問題

以下のサンフランシスコの気温データを描画し,時間単位で予測を行え.

sf = data.sf_temps()
sf.head()
temp date
0 47.8 2010-01-01 00:00:00
1 47.4 2010-01-01 01:00:00
2 46.9 2010-01-01 02:00:00
3 46.5 2010-01-01 03:00:00
4 46.0 2010-01-01 04:00:00

Githubのアクセス数データ

誤差が正規分布でない場合の取り扱い

github = data.github()
github.tail()
time count
950 2015-05-29 17:00:00 1
951 2015-05-29 19:00:00 1
952 2015-05-30 00:00:00 10
953 2015-05-30 09:00:00 1
954 2015-05-30 11:00:00 2
#px.line(github,x="time",y="count")

累積量を計算

cumsumで累積和(cumulative sum)が計算できる.

これを用いて予測する.

誤差が正規分布に近いことが確認できる.

他の方法(予測したいものの対数をとるなどの方法)だとうまくいかないことを確認せよ.

github["cumsum"] = github["count"].cumsum()
github.head()
time count cumsum
0 2015-01-01 01:00:00 2 2
1 2015-01-01 04:00:00 3 5
2 2015-01-01 05:00:00 1 6
3 2015-01-01 08:00:00 1 7
4 2015-01-01 09:00:00 3 10
github.rename(inplace=True,columns={"time":"ds","cumsum":"y"})
github.head()
ds count y
0 2015-01-01 01:00:00 2 2
1 2015-01-01 04:00:00 3 5
2 2015-01-01 05:00:00 1 6
3 2015-01-01 08:00:00 1 7
4 2015-01-01 09:00:00 3 10
# model = Prophet().fit(github)
# future = model.make_future_dataframe(periods=200, freq='H')
# forecast = model.predict(future)
# model.plot(forecast);
INFO:fbprophet:Disabling yearly seasonality. Run prophet with yearly_seasonality=True to override this.
#model.plot_components(forecast);
# forecast["error"] = forecast.yhat - github.y
# forecast["error"].plot()
<matplotlib.axes._subplots.AxesSubplot at 0x7fd1e15af6d0>
# forecast["error"].hist()
<matplotlib.axes._subplots.AxesSubplot at 0x7fe3d7c8b2b0>

傾向変化点

「上昇トレンドの株価が,下降トレンドに移った」というニュースをよく耳にするだろう.このように,傾向変動は,時々変化すると仮定した方が自然なのだ.Prophetでは,これを傾向の変化点として処理する.再び,Peyton Manningのデータを使う.

add_changepoints_to_plotを使うと、変化した点(日次)と傾向変動を図に追加して描画できる。引数は軸(axis),モデル(model),予測データフレーム(forecast)であり, 軸は図オブジェクトのgca(get current axis)メソッドで得る.

# df = pd.read_csv('http://logopt.com/data/peyton_manning.csv')
# model = Prophet().fit(df)
# future = model.make_future_dataframe(periods=366)
# forecast = model.predict(future)
# fig = model.plot(forecast)
# a = fp.add_changepoints_to_plot(fig.gca(),model,forecast);
INFO:fbprophet:Disabling daily seasonality. Run prophet with daily_seasonality=True to override this.

変化点の数を制御するためのパラメータはchangepoint_prior_scaleであり、既定値は0.05である。これを増やすと変化点が増え、予測の自由度が増すため予測幅が大きくなる。

# model = Prophet(changepoint_prior_scale=0.5).fit(df)
# future = model.make_future_dataframe(periods=366)
# forecast = model.predict(future)
# fig = model.plot(forecast)
# fp.add_changepoints_to_plot(fig.gca(), model, forecast);
INFO:fbprophet:Disabling daily seasonality. Run prophet with daily_seasonality=True to override this.

傾向変化点のリストをchangepoints引数で与えることもできる。以下の例では、1つの日だけで変化するように設定している。

# model = Prophet(changepoints=['2014-01-01']).fit(df)
# future = model.make_future_dataframe(periods=366)
# forecast = model.predict(future)
# fig = model.plot(forecast)
# fp.add_changepoints_to_plot(fig.gca(), model, forecast);
INFO:fbprophet:Disabling daily seasonality. Run prophet with daily_seasonality=True to override this.