Prophetは,Facebookが公開しているBayes推論を用いた予測パッケージである.

インストール

Macの場合

conda install -c conda-forge fbprophet

Google Colabの場合

!pip install fbprophet

Windowsにインストールする際には,C言語のコンパイラが必要;以下を参照

https://github.com/facebook/prophet

#!pip install pystan==2.19.1.1
#!pip install -U fbprophet
#!pip install -U vega_datasets

需要予測

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

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

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

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

予測の公理

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

  • 予測は予測のためならず

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

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

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

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

  • 予測は外れるもの

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

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

  • 集約すれば精度があがる

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

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

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

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

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

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

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

予測手法と使い分け

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

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

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

ベイズの公式

  • $P(A | B)$ : $B$ という事象が発生したきに $A$ が起こる条件付き確率
  • $P( A \cap B)$: 事象 $A$ と $B$ が同時に発生する確率
$$ P(A|B) = \frac{P( A \cap B)}{P(B)} $$

同様に, $$ P(B|A) = \frac{P( A \cap B)}{P(A)} $$

事前分布と事後分布

  • $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)

事後分布は,尤度と事前分布の積に比例する.

ベイズ推論

生成モデル(シミュレーション)ができれば,ベイズの公式から,事後分布の推定が可能

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

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

# 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)
INFO:numexpr.utils:Note: NumExpr detected 16 cores but "NUMEXPR_MAX_THREADS" not set, so enforcing safe limit of 8.
INFO:numexpr.utils:NumExpr defaulting to 8 threads.
INFO:fbprophet:Disabling daily seasonality. Run prophet with daily_seasonality=True to override this.
Initial log joint probability = -19.4685
<fbprophet.forecaster.Prophet at 0x7fe330a0eeb0>
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
      99        7975.9   0.000255252       131.906           1           1      131   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     199       7994.17    0.00448229        246.32           1           1      248   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     299       7996.79   0.000338229       123.492      0.8331      0.8331      370   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     313        7997.2   6.05911e-05       182.322   2.151e-07       0.001      428  LS failed, Hessian reset 
     364       7998.16   7.11253e-05       192.403   4.219e-07       0.001      530  LS failed, Hessian reset 
     399       7998.39   0.000591238       110.388           1           1      572   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     418       7999.13   5.70978e-05       151.423   1.809e-07       0.001      635  LS failed, Hessian reset 
     499       8001.17   0.000213455       76.8988           1           1      733   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     574        8001.8   3.01012e-05         98.37   3.835e-07       0.001      862  LS failed, Hessian reset 
     599        8001.8   4.20467e-06       72.0693           1           1      895   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     613       8002.31   0.000207208       407.536   8.267e-07       0.001      951  LS failed, Hessian reset 
     699       8003.12    0.00217504       91.3187           1           1     1056   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     799       8003.68   2.97403e-05       62.4093      0.2543      0.9659     1178   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     899       8003.89   3.28902e-06       66.6888           1           1     1308   
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
     903       8003.89   7.21402e-07       54.3237      0.0652           1     1315   
Optimization terminated normally: 
  Convergence detected: relative gradient magnitude is below tolerance

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

future = model.make_future_dataframe(periods=365)
future.tail()
ds
3265 2017-01-15
3266 2017-01-16
3267 2017-01-17
3268 2017-01-18
3269 2017-01-19

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

forecast = model.predict(future)
forecast[["ds", "yhat", "yhat_lower", "yhat_upper"]].tail()
ds yhat yhat_lower yhat_upper
3265 2017-01-15 8.214006 7.500984 8.940732
3266 2017-01-16 8.539078 7.809182 9.259313
3267 2017-01-17 8.326518 7.610270 9.085418
3268 2017-01-18 8.159179 7.449574 8.847648
3269 2017-01-19 8.171119 7.453488 8.911960

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)