Prophetとは
Prophet https://facebook.github.io/prophet/ は需要予測のためのパッケージである． ここでは，Prophetを用いた需要予測の方法について述べる． また，需要予測の基本原理と， Prophetの基礎になるベイズ推論をPyMCパッケージを用いて解説する．
需要予測
筆者は，以前は最適化の仕事は引き受けても予測の仕事は引き受けないことにしていた．予測は現場で長年働いている人の経験を加味して行うべきものであり，門外漢がいくらテクニックを駆使してもそれを超えることは難しいかと考えていたからだ．
しかし最近になって，需要予測は当たらないといけないという考えを捨てて，誤差の管理を行うことだと割り切って考えることにして，予測に積極的に関与するようになった．多くの最適化モデルは，予測がある程度合っているという前提で構築される．予測をいいかげんにされると，最適化モデル自体が役に立たないものになる．ゴミを入れればゴミが出てくるからだ．
誤差を管理することによって，需要予測を「点」で行うのではなく，「範囲」で行うことが可能になる．また，需要の分布も特定できるようになる．範囲内での最適化はロバスト最適化，確率分布を仮定した最適化は確率最適化という枠組みで解決可能になる．
以下では，サプライ・チェイン最適化に関連した需要予測の基本と重要なモデルについて解説する．
予測の公理
最近，サプライ・チェインの現場において，予測に関する多くの誤りが浸透していることに気づいた．ここでは，このような誤用を減らして正しい予測手法を適用するための，予測の公理について述べる．
- 予測は予測のためならず
しばしば，予測を目的として仕事をしている人たちを見かける．予測は，ほかの重要な意思決定を行うための基礎となる手段であり，予測そのものを目的としてはならない．実際には，予測よりもその誤差を評価することの方が重要である．誤差が増えているのか，減っているのか，その理由は何かを考えることが，需要予測の真の目的なのである．
たとえば，小売りの現場で需要予測を行うことは，在庫費用の削減や品切れ損出の回避などを目的としたものであり，予測の精度だけを問題にするのではなく，どの程度外れているのかという誤差の管理と，外れた場合の影響や，緊急発注などの回避手段とあわせて考える必要がある．
また，予測するだけでなく，なぜそのような値になったのかを究明することも重要である．需要が０という日が続いた場合には，それが，本当に需要がなかったのか，それとも在庫がないために売れなかったのか，陳列場所が悪かったために売れなかったのか，などをあわせて原因を分析する必要がある．
多くのメーカーでは，需要予測をするためだけの部署を設けているが，これもナンセンスなのでできるだけ早くほかの部署と連携をとるように改めるべきである．需要は当てるものではなくコントロールするものであり，需要予測を当てることだけを目標としている部署は，廃止すべきである．
- 予測は外れるもの
しばしば予測が当たったとか外れたとかいう言葉を現場の人から聞くが，経営はギャンブルではないので的中というのはありえない．この人なら当たるとかいうのは迷信であり，たまたま当たったときに声を張り上げて宣伝しているか，誤差が大きいにもかかわらず当たったと宣伝しているかの何れかである．サプライ・チェインからはちょっと外れるが，地震の予測（予知ともいう）も似たようなものであり，日本中のどこかで地震が発生すると予測し，そのうち1つが当たると的中と宣伝していたりする．ましてや，株や競馬の予想的中などの宣伝は，たまたま当たったときの結果だけを掲載し，外れたときのものを消去して，予想的中の証拠として提出していたりする．いまだに，こんな宣伝にだまされる人がいるのかと感心するが，社内で需要予測が的中する人がいるという会社も似たような詐欺に遭っていると言える．
重要なことは，どの程度外れたのかを時系列的に管理することである．誤差が増えている際には，その原因を追求し，予測手法を改善するなり，在庫を増やして品切れを回避するなりの行動をとるのが正しい方法である．
- 集約すれば精度が上がる
往々にして，個々のものの予測は難しいが，それをまとめたものの予測は容易になる．たとえば，特定の場所で特定のマグニチュードの地震が明日発生することを予測するのは難しいが，日本のどこかでマグニチュード4以上の地震が来年発生することは容易に予測できる．前者はほぼ$0$％であるが，後者はほぼ$１００$%の精度で予測できる．これが集約の力である．
サプライ・チェインでも同様であり，商品を集約して商品群にすれば精度があがり，個々の店舗での売上でなく，地域内のすべての店舗の売上を集約すれば精度が上がる． 時間軸でも同様であり，１時間以内の需要を予測精度は，日，週，月，年単位と集約していくにしたがってあがっていく．
予測精度だけを議論するのであれば，どんどん集約した方が得であるが，もちろんどんどん役に立たなくなる．前述したように，重要なことは，その予測を何に使うかであり，使用法にあわせて「適切に」集約を行うことである．
また，製品設計や在庫地点を考慮することによって，物理的に集約を行うこともできる．たとえば，製品のモジュール化や遅延差別化やリスク共同管理がこれに相当する．これらはモダンなサプライ・チェインの基本戦略であり，そのすべてがこの単純な公理（集約した方が予測精度が上がる）に基づくものであることは興味深い．
- 目先の需要は当たりやすく，遠い未来の需要は予測しにくい
明日の天気はある程度予測できるが，1年後の天気は予測しにくい．不確実性は時間が経過するにしたがって増大するからである．需要予測も同様である．たとえば，カップラーメンなどは一部の定番品を除いて，来年店頭に並んでいるかどうかも怪しい．
近い未来の予測は，短期のオペレーショナルな意思決定に用いるので，比較的正確性が必要であるが，遠い未来の予測は長期のストラテジックな意思決定に用いるので，おおよそで構わない．さらに，長期の意思決定の際には，データは集約して行われるので，予測精度も上がる． 要は，目的のために適切な精度で管理できるように，時間軸を集約することが推奨される．たとえば，近い未来は予測精度がよいので，集約をせずに日単位で予測し，未来に行くに従って時間軸の集約を行い，週単位，月単位，年単位としていく訳である．このテクニックはテレスコーピングと呼ばれ，時間軸を含んだ実際問題を解くときによく用いられる．
予測手法と使い分け
古典的な予測手法は statsmodelsに入っている。種々の指数平滑法の拡張や、ARIMAなどは簡単にできる。 いずれも高速なので、補助的なデータがない時系列データに対して、簡易的な予測をしたい場合には便利だが、予測精度は期待すべきではない。
機械（深層）学習を用いた予測は、 scikit-learn（最近ではpycaret）やfastaiのTabularモデルやLSTMを用いることによって容易にできる。 補助的なデータが豊富な場合には、これらの手法が推奨される。
補助的なデータが少ない時系列データの場合には、prophetを用いたベイズ推論が推奨される。ベイズ推論だと、予測を点でするのではなく、不確実性の範囲まで得られるので、 サプライ・チェインのモデルと相性が良い。例えば、在庫モデルにおいては、不確実性の情報が不可欠だからだ。
ベイズ推論
ここでは， ベイズ推論 (Bayesian inference) をPyMCパッケージ https://www.pymc.io/ を用いて解説する． PyMCはGoogle Colabにプレインストールされている．
事前分布と事後分布と尤度関数
- $P(\bar{B})$: 事象 $B$ が発生しない確率
- $P(B) + P(\bar{B}) =1$
ベイズの公式から， $$ P(B|A) = \frac{P(A|B) P(B) }{P(A)} = \frac{P( A|B) p(B) }{P(A|B) P(B) + P(A|\bar{B}) P(\bar{B}) } $$
$B$ をモデルのパラメータ $\theta$， $A$ を（観測された）データとすると， $$ P(\theta| data) = \frac{P( data|\theta) P(\theta) }{P(data)} $$
- $P(\theta)$: 事前分布 (prior)
- $P( data|\theta)$ : パラメータ $\theta$ に対する尤度関数 $L(\theta)$ (likelihood)
- $P(\theta| data)$ : 事後分布 (posterior)
事後分布は，尤度と事前分布の積に比例する．
$$ P(\theta| data) \propto P( data|\theta) P(\theta) = Likelyhood \times Prior $$import arviz as az
import numpy as np
import pymc as pm
alpha, sigma = 1, 1
beta = [1, 2.5]
size = 100
X1 = np.random.randn(size)
X2 = np.random.randn(size) * 0.2
Y = alpha + beta[0] * X1 + beta[1] * X2 + np.random.normal(size=size) * sigma
model = pm.Model()
with model:
# 事前分布
alpha = pm.Normal("alpha", mu=0, sigma=10)
beta = pm.Normal("beta", mu=0, sigma=10, shape=2)
sigma = pm.HalfNormal("sigma", sigma=1)
# 生成モデル
mu = alpha + beta[0] * X1 + beta[1] * X2
# 尤度関数
Y_obs = pm.Normal("Y_obs", mu=mu, sigma=sigma, observed=Y)
MCMC (Markov chain Monte Carlo)法
実際の分布の推定には， MCMC (Markov chain Monte Carlo)法を利用する． <!-- 詳細は，以下を参照．
https://en.wikipedia.org/wiki/Markov_chain_Monte_Carlo -->
MCMC法は，以下に示すように改悪も許した確率的な探索を行うことによって，分布のサンプリングを行う．
ベイズ線形回帰の事後分布をサンプリングによって生成する．
with model:
# 事後分布を 1000 個サンプル（時間がかかるのでしばらく待つ）
idata = pm.sample(draws=1000)
az.plot_trace(idata, combined=True, figsize=(10,10));
az.summary(idata, round_to=2)
元のデータは，以下のコードによって生成していたので，良い近似になっていることが確認できる．
alpha, sigma = 1, 1
beta = [1, 2.5]
#import plotly.io as pio
#pio.renderers.default = "colab"
#!pip install prophet
#!pip install -U vega_datasets
import pandas as pd
from prophet import Prophet
from vega_datasets import data
import plotly.express as px
import prophet.plot as fp
import plotly
df = pd.read_csv("http://logopt.com/data/peyton_manning.csv")
df.head()
Prophetモデルのインスタンスを生成し，fitメソッドで学習（パラメータの最適化）を行う．fitメソッドに渡すのは，上で作成したデータフレームである．このとき、ds(datestamp)列に日付（時刻）を、y列に予測したい数値を入れておく必要がある （この例題では，あらかじめそのように変更されている）．
model = Prophet()
model.fit(df)
make_future_dataframeメソッドで未来の時刻を表すデータフレームを生成する。既定値では、予測で用いた過去の時刻も含む。 引数は予測をしたい期間数periodsであり，ここでは、１年後（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)
plotly.offline.plot(fig);
co2 = data.co2_concentration()
co2.head()
fig = px.line(co2,x="Date",y="CO2")
plotly.offline.plot(fig);
列名の変更には，データフレームのrenameメソッドを用いる．引数はcolumnsで，元の列名をキーとし，変更後の列名を値とした辞書を与える．また，元のデータフレームに上書きするために，inplace引数をTrueに設定しておく．
co2.rename(columns={"Date":"ds","CO2":"y"},inplace=True)
co2.head()
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)
plotly.offline.plot(fig);
passengers = pd.read_csv("http://logopt.com/data/AirPassengers.csv")
passengers.head()
fig = px.line(passengers,x="Month",y="#Passengers")
plotly.offline.plot(fig);
passengers.rename(inplace=True,columns={"Month":"ds","#Passengers":"y"})
passengers.head()
季節変動を乗法的に変更するには， モデルの 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);
結果から，乗法的季節変動の方が，良い予測になっていることが確認できる．
retail = pd.read_csv(`http://logopt.com/data/retail_sales.csv`)
retail.head()
climate = data.seattle_temps()
climate.head()
このデータは， date 列に日付と1時間ごとの時刻が， temp 列に気温データが入っている．
Prophetは， 日別でないデータも扱うことができる。 date列のデータ形式は、日付を表すYYYY-MM-DDの後に時刻を表すHH:MM:SSが追加されている。
未来の時刻を表すデータフレームは、make_future_dataframeメソッドで生成するが、このとき引数freqで時間の刻みを指定する。
ここでは1時間を表す H を指定する。
climate["Date"] = pd.to_datetime(climate.date)
climate.rename(columns={"Date":"ds","temp":"y"},inplace=True)
model = Prophet().fit(climate)
future = model.make_future_dataframe(periods=200, freq="H")
forecast = model.predict(future)
model.plot(forecast);
因子ごとに予測値を描画すると，傾向変動と週次の季節変動の他に，日次の季節変動（1日の気温の変化）も出力される．
model.plot_components(forecast);
sf = data.sf_temps()
sf.head()
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);
変化点の数を制御するための引数は 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);
傾向変化点のリストを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);
sp500 = data.sp500()
sp500.tail()
sp500.rename(inplace=True,columns={"date":"ds","price":"y"})
model = Prophet(changepoint_prior_scale=0.5, changepoint_range=0.95,yearly_seasonality=5).fit(sp500)
future = model.make_future_dataframe(periods=200, freq=`D`)
forecast = model.predict(future)
model.plot(forecast);
model.plot_components(forecast);
from prophet.plot import add_changepoints_to_plot
fig = model.plot(forecast)
a = add_changepoints_to_plot(fig.gca(), model, forecast)
stocks = data.stocks()
stocks.tail()
fig = px.line(stocks,x="date",y="price",color="symbol")
plotly.offline.plot(fig);
以下では，マイクロソフトの株価を予測してみる．
msft = stocks[ stocks.symbol == "MSFT"]
msft.head()
msft = msft.rename(columns={"date":"ds","price":"y"})
msft.head()
model = Prophet(changepoint_prior_scale=0.5, changepoint_range=0.95,yearly_seasonality=5).fit(msft)
future = model.make_future_dataframe(periods=200, freq=`D`)
forecast = model.predict(future)
model.plot(forecast);
model.plot_components(forecast);
stocks = data.stocks()
stocks.head()
logistic = pd.read_csv("http://logopt.com/data/logistic.csv")
model = Prophet(growth="logistic")
logistic["cap"] = 8.5
model.fit(logistic)
future = model.make_future_dataframe(periods=1826)
future[`cap`] = 8.5
forecast = model.predict(future)
model.plot(forecast);
outliers1 = pd.read_csv("http://logopt.com/data/outliers1.csv")
model = Prophet()
model.fit(outliers1)
future = model.make_future_dataframe(periods=1096)
forecast = model.predict(future)
model.plot(forecast);
2010年のデータを除外することによって、予測が改善される。
outliers1.loc[(outliers1[`ds`] > `2010-01-01`) & (outliers1[`ds`] < `2011-01-01`), `y`] = None
model =Prophet()
model.fit(outliers1)
forecast = model.predict(future)
model.plot(forecast);
上では、外れ値を除外することによって予測が改善されたが、これがいつでも成立するとは限らない。以下の例では2015年6月に外れ値が観察される。
outliers2 = pd.read_csv("http://logopt.com/data/outliers2.csv")
model = Prophet()
model.fit(outliers2)
future = model.make_future_dataframe(periods=1096)
forecast = model.predict(future)
model.plot(forecast);
上の予測では2015年6月のデータを予測に用いず，外れ値として処理しているので，外れ値を除外すると予測の幅が広がる。
outliers2.loc[(outliers2[`ds`] > `2015-06-01`) & (outliers2[`ds`] < `2015-06-30`), `y`] = None
model = Prophet()
model.fit(outliers2)
future = model.make_future_dataframe(periods=1096)
forecast = model.predict(future)
model.plot(forecast);
df = pd.read_csv(`http://logopt.com/data/peyton_manning.csv`)
df.head()
playoffs = pd.DataFrame({
`holiday`: `playoff`,
`ds`: pd.to_datetime([`2008-01-13`, `2009-01-03`, `2010-01-16`,
`2010-01-24`, `2010-02-07`, `2011-01-08`,
`2013-01-12`, `2014-01-12`, `2014-01-19`,
`2014-02-02`, `2015-01-11`, `2016-01-17`,
`2016-01-24`, `2016-02-07`]),
`lower_window`: 0,
`upper_window`: 1,
})
superbowls = pd.DataFrame({
`holiday`: `superbowl`,
`ds`: pd.to_datetime([`2010-02-07`, `2014-02-02`, `2016-02-07`]),
`lower_window`: 0,
`upper_window`: 1,
})
holidays = pd.concat((playoffs, superbowls))
holidays.head()
引数holidaysで休日を表すデータフレームを与えることによって、特別なイベントを考慮した予測を行うことができる。
model= Prophet(holidays=holidays)
model.fit(df)
future = model.make_future_dataframe(periods=365)
forecast = model.predict(future)
プレーオフやスーパーボールなどのイベント効果がある日だけ抜き出してデータフレームを表示する．
forecast[(forecast[`playoff`] + forecast[`superbowl`]).abs() > 0][
[`ds`, `playoff`, `superbowl`]][-10:]
因子別に描画を行うと，イベントによって変化した量が描画される（上から2番目）．
model.plot_components(forecast);
model = Prophet(holidays=holidays)
model.add_country_holidays(country_name=`US`)
model.fit(df)
追加された休日名をtrain_holiday_names属性で確認しておく．
model.train_holiday_names
米国の休日を考慮して予測を行い，因子別に描画してみる．上から2番目が，休日に対する影響を表している．
forecast = model.predict(future)
model.plot_components(forecast);
def nfl_sunday(ds):
date = pd.to_datetime(ds)
if date.weekday() == 6 and (date.month > 8 or date.month < 2):
return 1
else:
return 0
df = pd.read_csv(`http://logopt.com/data/peyton_manning.csv`)
df[`nfl_sunday`] = df[`ds`].apply(nfl_sunday)
model = Prophet()
model.add_regressor(`nfl_sunday`)
model.fit(df)
future = model.make_future_dataframe(periods=365)
future[`nfl_sunday`] = future[`ds`].apply(nfl_sunday)
forecast = model.predict(future)
model.plot_components(forecast);
model = Prophet(weekly_seasonality=False)
model.add_seasonality(name=`monthly`, period=30.5, fourier_order=5)
model.fit(df)
future = model.make_future_dataframe(periods=365)
forecast = model.predict(future)
model.plot_components(forecast);
def is_nfl_season(ds):
date = pd.to_datetime(ds)
return (date.month > 8 or date.month < 2)
df[`on_season`] = df[`ds`].apply(is_nfl_season)
df[`off_season`] = ~df[`ds`].apply(is_nfl_season)
model = Prophet(weekly_seasonality=False)
model.add_seasonality(name=`weekly_on_season`, period=7, fourier_order=3, condition_name=`on_season`)
model.add_seasonality(name=`weekly_off_season`, period=7, fourier_order=3, condition_name=`off_season`)
model.fit(df)
future = model.make_future_dataframe(periods=365)
future[`on_season`] = future[`ds`].apply(is_nfl_season)
future[`off_season`] = ~future[`ds`].apply(is_nfl_season)
forecast =model.predict(future)
model.plot_components(forecast);
df = pd.read_csv(`http://logopt.com/data/peyton_manning.csv`)
model= Prophet(holidays=holidays, holidays_prior_scale=0.05)
model.fit(df)
future = model.make_future_dataframe(periods=365)
forecast = model.predict(future)
forecast[(forecast[`playoff`] + forecast[`superbowl`]).abs() > 0][
[`ds`, `playoff`, `superbowl`]][-10:]
スーパーボール（superbowl）の効果が抑制されていることが見てとれる．
同様に，季節変動の影響はseasonality_prior_scaleを小さくすることによって抑制できる。
co2 = data.co2_concentration()
co2.rename(columns={"Date":"ds","CO2":"y"},inplace=True)
co2.head()
model = Prophet(interval_width=0.95)
model.fit(co2)
future = model.make_future_dataframe(periods=200, freq=`M`)
forecast = model.predict(future)
model.plot(forecast);
forecast = Prophet(interval_width=0.5).fit(co2).predict(future)
model.plot(forecast);
季節変動に対する不確実性を予測するためには、マルコフ連鎖モンテカルロ法を行う必要がある、そのためには、引数 mcmc_samples をシミュレーションの反復回数に設定する。
このパラメータの既定値は $0$ である。
これによって、既定値の最大事後確率（MAP）推定の代わりにマルコフ連鎖モンテカルロ法によるサンプリングが行われる。これは、非常に時間がかかることもある。 要因別に図を描画してみると、季節変動に対しても不確実性の幅が示されていることが確認できる。
model = Prophet(mcmc_samples=100)
forecast = model.fit(co2).predict(future)
model.plot_components(forecast);
検証と誤差の評価
Prophetでは、予測の精度を検証するための仕組みが組み込まれている。例として、Peyton Manningのデータセットを用いる。 このデータセットは、全部で $2905$ 日分のデータで構成されている。
交差検証のためにはcross_validationを用いる。
引数は以下の通り．
- model: 予測を行うモデル；事前にfitメソッドで学習しておく必要がある．
- horizon : 計画期間（予測を行う期間）
- period : 予測の間隔；省略するとhorizonの半分が代入される．
- initial : 交差検証を開始する最初の期；省略するとhorizonの3倍が代入される．
以下の例では、initialが $730$ 日なので、$729$ 日までの情報を用いて、その後の $365$(horizon)日の予測を行い、本当の値との誤差を評価し、 次いで $730+180$(period)日までの情報を用いて、その後の $365$ 日の予測を行い評価し、という手順を最後の日まで繰り返す。 $(2905-730-365)/180 = 10.05$ であるので、 $11$ 回の予測を行い評価することになる。cross_validationは、交差検証用のデータフレームを返す。
最初の検証は $730$ 日後である 2010-2-15(cutoff)までのデータを用いて，2010-2-16から $365$ (horizon)日分の予測で行われ、次の検証はその $180$(period)日後である2010-08-14日から行われる。最後の検証は2015-01-20日までのデータを用いて2016-01-20日まで行われる。
df = pd.read_csv(`http://logopt.com/data/peyton_manning.csv`)
model = Prophet()
model.fit(df);
from prophet.diagnostics import cross_validation
df_cv = cross_validation(model, initial=`730 days`, period=`180 days`, horizon = `365 days`)
df_cv.head()
performance_metrics を用いてメトリクス（評価尺度）を計算する。評価尺度は、 平均平方誤差(mean squared error: MSE), 平均平方誤差の平方根 (root mean squared error: RMSE), 平均絶対誤差 (mean absolute error: MAE), 平均絶対パーセント誤差 (mean absolute percent error : MAPE), yhat_lower とyhat_upper の間に入っている割合（被覆率: coverage) である。
既定値では予測期間の最初の$10$%は除外して示される。これは、引数rolling_window によって変更できる。
from prophet.diagnostics import performance_metrics
df_p = performance_metrics(df_cv, rolling_window=0.1)
df_p.head()
評価尺度は plot_cross_validation_metricで可視化できる。以下では平均絶対パーセント誤差(MAPE)を描画している．
from prophet.plot import plot_cross_validation_metric
plot_cross_validation_metric(df_cv, metric=`mape`);
主なパラメータと既定値
以下にProphetの主要なパラメータ（引数）とその既定値を示す．
- growth=
linear:傾向変動の関数．既定値は線形．ロジスティック曲線にするにはlogisticに設定 - changepoints=None : 傾向変更点のリスト
- changepoint_range = $0.8$ : 傾向変化点の候補の幅（先頭から何割を候補とするか）
- n_changepoints=25 : 傾向変更点の数
- yearly_seasonality=
auto: 年次の季節変動を考慮するか否か - weekly_seasonality=
auto: 週次の季節変動を考慮するか否か - daily_seasonality=
auto: 日次の季節変動を考慮するか否か - holidays=None : 休日のリスト
- seasonality_prior_scale= $10.0$ : 季節変動の事前分布のスケール値（パラメータの柔軟性を表す）
- holidays_prior_scale= $10.0$ : 休日のの事前分布のスケール値（パラメータの柔軟性を表す）
- changepoint_prior_scale= $0.05$ : 傾向変更点の事前分布のスケール値（パラメータの柔軟性を表す）
- mcmc_samples= $0$ : MCMC法のサンプル数
- interval_width= $0.8$ : 不確実性の幅
- uncertainty_samples= $1000$ : 不確実性の幅を計算する際のサンプル数