# 需要予測モジュールforecastと予測システム ABD-Forecast ![](06forecast_files/figure-commonmark/cell-3-output-1.jpeg) \#hide \## 解釈可能な予測（最適決定木） 解釈可能な人工知能の必要性から，決定木の作成に最適化を適用する． 特徴量が連続の場合の分類は，数理最適化モデルとして記述できる． \#hide \### 最適分類木 Optimal Trees for Prediction and Prescription by Jack William Dunn　（MIT PhD Thesis) 2018 変数 - d\[t\] : ブランチtで分枝をするとき１ - b\[t\] : ブランチtでの閾値 - a\[j,t\]: ブランチtが特徴jで分けるとき1 - l\[t\]: 葉が空でないとき1 - z\[i,t\]: データiが葉tに入るとき１ - N\[k,t\] : 葉tに入ったクラスkのデータ数 - c\[k,t\]: 葉tがクラスkに分類されるとき1 - NN\[t\]: 葉tに入ったデータの総数 - L\[t\]: 葉tにおける誤りの数 \#hide \### ２値分類のための最適決定木（すべてのデータはカテゴリーデータと仮定） Optimal Decision Trees for Categorical Data via Integer Programming Oktay Gu ̈nlu ̈k · Jayant Kalagnanam · Minhan Li · Matt Menickelly · Katya Scheinberg 変数 - v\[g,t\]: グループ g がブランチ t で分割対象となるとき 1 - z\[j,t\]: 特徴 j がブランチ t で左の子に行くとき 1 - L\[i,t\]: データ i がブランチ t で左に行くとき 1 - R\[i,t\]: データ i がブランチ t で右に行くとき 1 - c\[i,t\]: データ i が葉 t に行くとき １ SCOP C\[i\] domain Leaf 定式化 - 各ブランチでは1つのグループが選択される： ``` python for t in Branch: m.addConstr( quicksum(v[g,t] for g in G) == 1 ) ``` - 特徴 j　がブランチ t で選ばれるなら，その特徴を含んだグループが選ばれなければならない： ``` python for t in Branch: for j in range(p): m.addConstr( z[j,t] <= v[group[j], t] ) ``` - 追加制約 =\> グループが選ばれるなら，その中の特徴は1つ以上選択 ``` python for t in Branch: for g in G: m.addConstr( quicksum(z[j,t] for j in J[g]) >= v[g, t] ) ``` - 選択された特徴をもつデータ i は左の子に行き，それ以外は右の子に行く： ``` python for i in range(n): for t in Branch: m.addConstr( L[i,t] == quicksum( X[i][j]*z[j,t] for j in range(p)) ) m.addConstr( L[i,t] + R[i,t]== 1 ) ``` - データ i が葉 t に含まれるなら，その葉に行く先祖点において左（右）に行くデータでなければならない： ``` python for t in Leaf: for i in range(n): for ansestor in AL[t]: m.addConstr( c[i,t] <= L[i,ansestor] ) for ansestor in AR[t]: m.addConstr( c[i,t] <= R[i,ansestor] ) ``` - データ i は何れかの葉に含まれなければならない： ``` python for i in range(n): m.addConstr( quicksum(c[i,t] for t in Leaf) == 1 ) ``` - 目的関数： 正しいクラスの葉に割り当てられたデータの数の合計を最大化： ``` python m.setObjective( quicksum( c[i,t] for i in DataPos for t in LeafPos ) + quicksum( c[i,t] for i in DataNeg for t in LeafNeg ), GRB.MAXIMIZE ) ``` ## 需要予測 サプライ・チェインにおける需要の予測のための入力は，一般的には以下の形式をもつデータである．
日付時刻 顧客 製品 需要量 売上 付加情報1 付加情報2 …
2020年1月3日（何時何分） A商店 製品α 10 2600 特売日 …
2020年1月３日 B商店 製品β 6 6000 ポイント10倍 …
… … … … … … … …
最初の列は日付時刻型のデータが入る．通常は，日付だけだが，何時何分に購入されたかの情報も入れたい場合には，日付時刻型で入力する． この列をインデックスにすることによって，時間の単位にデータを集約することができる． 2番目と3番目の列には，それぞれ顧客と製品（商品）のデータを入れる． 次の列は，予測をしたいデータ（ターゲット，従属変数）を入れる．通常は，需要量（何個売れたか）を予測するが，それに販売金額を乗じた売上を予測する場合もある． 付加的な情報がある場合には，その次の列に入力する．たとえば，特売日やポイント付加などの情報を入力する． 販売金額が変動する場合には，価格を入れることもある． 予測を行う際には，顧客と製品別に集約し，時系列データに変換してから予測手法を適用する場合が多い． その場合には，（販売金額が一定なら），需要を予測しても，売上を予測しても同じである． 顧客・製品別に予測を行うのは，簡便性のためであるが，顧客ごとの特徴や製品ごとの特徴を考慮できないことがある． さらに，顧客や製品から抽出できる付加情報を，付加情報列に追加しても良い．たとえば，同じ系列の顧客の場合には系列名を追加しておくと，特徴が認識しやすくなる． 日付は，自動的に特徴を表す列に変換して扱うこともある．ここでは，以下の列を追加する．（時刻も考慮する場合には，さらに時間，分，秒の列を追加する．） - ‘Year’: 年 - ‘Month’: 月 - ‘Day’ : 日 - ‘Dayofweek’: 曜日 - ‘Dayofyear’: 年初からの経過日数 - ‘Is_month_end’: 月の最終日 - ‘Is_month_start’: 月の開始日 - ‘Is_quarter_end’: 四半期の終了日 - ‘Is_quarter_start’: 四半期の開始日 - ‘Is_year_end’: 年の最終日 - ‘Is_year_start’: 何の開始日 - ‘Elapsed’: 経過時間 ここで用いるサンプルデータは， データ生成モジュールのgenerate_demand_with_promo関数を用いて生成したものである． 以下のように読み込んで使用する． ``` python demand_with_promo_df = pd.read_csv(folder + "demand_with_promo.csv", index_col=0) demand_with_promo_df.head() ```
date cust prod promo_0 promo_1 demand
0 2019-01-01 札幌市 A 0 0 2
1 2019-01-01 札幌市 B 0 0 0
2 2019-01-01 札幌市 C 0 0 2
3 2019-01-01 札幌市 D 0 0 0
4 2019-01-01 札幌市 E 0 0 1
## 指数平滑法 サプライ・チェインの古典手法に、指数平滑法と呼ばれる予測法がある。 航空機の乗客数の時系列データを用いて、様々な手法の比較を行う。 ## Holt-Wintersの季節変動と傾向変動を考慮した指数平滑法 パラメータ - *y**t*: 期 *t* の実現値 - *s**t*: 基本水準係数 - *b**t*: 傾向変動係数 - *c**t*: 季節変動係数 - *ŷ**t* + *m*: 期 *t* における *m*期後の予測値 - *L* ：季節変動の周期 (seasonal_periods) - *α*: 平滑化定数 (smoothing_level) - *β*: 傾向変動に対する平滑化定数 (smoothing_slope) - *γ*: 季節変動に対する平滑化定数 (smoothing_seasonal) 加法的季節変動： $$ \begin{array}{lll} &s_0 & =& y_0\\ &s\_{t} & =& \alpha (y\_{t}-c\_{t-L}) + (1-\alpha)(s\_{t-1} + b\_{t-1})\\ &b\_{t} &=& \beta (s_t - s\_{t-1}) + (1-\beta)b\_{t-1}\\ &c\_{t} &=& \gamma (y\_{t}-s\_{t})+(1-\gamma)c\_{t-L}\\ &\hat{y}\_{t+m} &=& s_t + m b_t+c\_{t-L+1+(m-1)\mod L} \end{array} $$ 乗法的季節変動： $$ \begin{array}{lll} s_0 &=& y_0\\ s\_{t} &=& \alpha (y\_{t}/c\_{t-L}) + (1-\alpha)(s\_{t-1} + b\_{t-1})\\ b\_{t} &=& \beta (s_t - s\_{t-1}) + (1-\beta)b\_{t-1}\\ c\_{t} &=& \gamma (y\_{t}/s\_{t})+(1-\gamma)c\_{t-L}\\ \hat{y}\_{t+m} &=& (s_t + m b_t) c\_{t-L+1+(m-1)\mod L} \end{array} $$ 手法 - SimpleExpSmoothing: 単純指数平滑法（ *γ* のみ用いる） - Holt: Holt法 (*β* と *γ* を用いる） - ExponentialSmoothing: Holt-Wintersの指数平滑法（すべて用いる） 最後の2年間（24ヶ月）をテストデータとする。Holt-Winters法に対しては、パラメータは自動調整とし， Box-Cox変換をしてデータを正規化しておく（引数のuse_boxcoxをTrueにする）． ``` python passengers = pd.read_csv("http://logopt.com/data/AirPassengers.csv") passengers["Month"] = pd.to_datetime(passengers.Month) # Monthの列を日付時刻型に変換 passengers.set_index("Month", inplace=True) #月をインデックスとする passengers.index.freq ="MS" #月のはじめを頻度にする passengers.head() fit1 = SimpleExpSmoothing(passengers[:-24]).fit(smoothing_level=0.2,optimized=False) fit2 = Holt(passengers[:-24]).fit(smoothing_level=0.8, smoothing_slope=0.2, optimized=False) fit3 = ExponentialSmoothing(passengers[:-24],seasonal_periods=12, trend='add', seasonal='add').fit() start, end = len(passengers)-24, len(passengers) predict1 = fit1.predict(start, end) predict2 = fit2.predict(start, end) predict3 = fit3.predict(start, end) ``` ### SARIMAX で予測 ARIMA (autoregressive integrated moving average)モデル *y**t* − *y**t* − *d* = *c* + *ϕ*1*y**t* − 1 + *ϕ*2*y**t* − 2 + ⋯ + *ϕ**p**y**t* − *p* + *ε**t* + *θ*1*ε**t* − 1 + ⋯ + *θ**q**ε**t* − *q* パラメータ - 自己回帰 AR (autoregressive): *p*； 過去 *p* 期の値の多項式で予測 - 差分 I (integrated): *d*； *d* 期前との差分をとることによって傾向変動を予測 - 移動平均 MA (moving average): *q*； 過去 *q* 期の誤差の多項式で予測 季節変動の追加 SARIMA (seasonal ARIMA) パラメータ - 季節自己回帰 SAR: *P* - 季節差分 SI : *D* - 季節移動平均 SMA: *Q* - 周期: *L* SARIMA**X** 外生変数 (e**X**ogenous variables) の追加 ``` python mod = SARIMAX(passengers[:-24], order=(1, 1, 1), seasonal_order=(1, 1, 1, 12), enforce_stationarity=False, enforce_invertibility=False) fit4 = mod.fit() predict4 = fit4.predict(start,end) ``` ``` python df = passengers.reset_index() df.rename(columns={"Month": "ds", "#Passengers": "y"}, inplace=True) ``` ``` python model = Prophet(daily_seasonality=False, weekly_seasonality=False, yearly_seasonality=True,seasonality_mode='multiplicative') model.fit(df[:-24]) future = model.make_future_dataframe(periods=24, freq="MS") forecast = model.predict(future) forecast.head() ```
ds trend yhat_lower yhat_upper trend_lower trend_upper multiplicative_terms multiplicative_terms_lower multiplicative_terms_upper yearly yearly_lower yearly_upper additive_terms additive_terms_lower additive_terms_upper yhat
0 1949-01-01 115.732114 92.193672 116.702905 115.732114 115.732114 -0.097351 -0.097351 -0.097351 -0.097351 -0.097351 -0.097351 0.0 0.0 0.0 104.465421
1 1949-02-01 117.417721 89.359809 112.239922 117.417721 117.417721 -0.142170 -0.142170 -0.142170 -0.142170 -0.142170 -0.142170 0.0 0.0 0.0 100.724450
2 1949-03-01 118.940205 107.671109 130.977234 118.940205 118.940205 0.004902 0.004902 0.004902 0.004902 0.004902 0.004902 0.0 0.0 0.0 119.523276
3 1949-04-01 120.625813 104.280365 128.476356 120.625813 120.625813 -0.031312 -0.031312 -0.031312 -0.031312 -0.031312 -0.031312 0.0 0.0 0.0 116.848788
4 1949-05-01 122.257045 106.734658 130.131457 122.257045 122.257045 -0.027676 -0.027676 -0.027676 -0.027676 -0.027676 -0.027676 0.0 0.0 0.0 118.873491
``` python passengers["Exp.Smooth."] = predict1 passengers["Holt"] = predict2 passengers["Holt-Winter"] = predict3 passengers["SARIMAX"] = predict4 passengers["Prophet"] = forecast.loc[:,"yhat"].values tacked_df = pd.DataFrame(passengers[-24:].stack(),columns=["num"]).reset_index() fig = px.line(stacked_df,x="Month",y="num",color="level_1") plotly.offline.plot(fig); ``` ![](06forecast_files/figure-commonmark/cell-18-output-1.png) ``` python methods = ["Exp.Smooth.", "Holt", "Holt-Winter", "SARIMAX", "Prophet"] #passengers[-24:] target = "#Passengers" metrics_list =[ [] for m in methods] for i, m in enumerate(methods): y_pred = passengers[m][-24:] y_true = passengers[target][-24:] metrics_list[i] =[ mean_squared_error(y_true, y_pred), mean_absolute_error(y_true, y_pred), mean_squared_log_error(y_true, y_pred), median_absolute_error(y_true, y_pred), max_error(y_true, y_pred) ] metrics = pd.DataFrame({m:metrics_list[i] for i,m in enumerate(methods)}, index=["mean_squared", "mean_absolute", "mean_squared_log", "median_abusolute", "max"]) metrics ```
Exp.Smooth. Holt Holt-Winter SARIMAX Prophet
mean_squared 11560.124288 90391.493041 1278.707087 5719.201415 936.384334
mean_absolute 82.410342 270.694737 31.076413 69.895789 25.504026
mean_squared_log 0.055542 1.431378 0.005957 0.031673 0.004042
median_abusolute 51.102678 254.302112 29.663906 62.913468 24.663997
max 247.102678 517.826696 74.254745 129.229564 74.062530
## Prophet ### 顧客の需要予測 過去の需要系列dfをもとに、Prophetによる予測を行う関数 forecast_demand を準備する。 比較として、（需要の定常性を仮定して）平均と標準偏差も算出する。 引数： - demand_df: 需要の系列を表すデータフレーム。顧客Cust、製品Prod,日時Date、需要量demandの列をもつものと仮定する。 - cust_selected: 顧客のリスト - prod_selected: 製品のリスト - agg_period: 集約を行う期。例えば１週間を１期とするなら”1w”とする。既定値は1日（“1d”） - forecast_period: 予測を行う未来の期数。既定値は1 - cumulative: 累積需要量を用いて予測するとき True；既定値はFalse - lb: 下限値（傾向変動がロジスティック曲線の場合に有効になる．）既定値は0. - ub: 上限値（傾向変動がロジスティック曲線の場合に有効になる．）既定値は Noneで，その場合にはデータの最大値が設定される． - kwargs: その他のProphetの引数 Prophetの引数と既定値： - growth=‘linear’ :傾向変動の関数．規定値は線形．ロジスティック曲線にするには，引数を’logistic’に設定する． また傾向変動をなくすためには引数を ’flat’に設定する． - 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.80 : 不確実性の幅 - uncertainty_samples=1000 : 不確実性の幅を計算する際のサンプル数 返値： 以下の情報を列情報にもつデータフレーム - dic: 顧客と製品のタプルをキーとし、モデルと予測オブジェクトのタプルを値とした辞書 ------------------------------------------------------------------------ source ### forecast_demand ``` python def forecast_demand( demand_df, cust_selected, prod_selected, agg_period:str='1d', forecast_periods:int=1, cumulative:bool=False, lb:float=0.0, ub:NoneType=None, kwargs:VAR_KEYWORD ): ``` *Prophetを用いた需要予測（顧客と製品に対して）* ### forecast_demand関数の使用例 ``` python demand_df = pd.read_csv(folder+"demand.csv") #cust_selected = [ demand_df["cust"].iloc[350] ] #prod_selected = [ demand_df["prod"].iloc[350] ] cust_selected = ["前橋市"] prod_selected = ["F"] print(cust_selected, prod_selected) dic = forecast_demand(demand_df, cust_selected, prod_selected, agg_period ="1d", forecast_periods =10, cumulative=False, lb=0., ub=None, growth ="logistic", daily_seasonality=False, weekly_seasonality=True, yearly_seasonality=True) print(cust_selected[0], prod_selected[0]) if dic[ cust_selected[0], prod_selected[0] ] is not None: model, forecast = dic[ cust_selected[0], prod_selected[0] ] model.plot(forecast); else: print("too few data") ``` ['前橋市'] ['F'] 前橋市 F too few data ## 特徴量を入れたProphetによる需要予測 過去の需要系列demand_dfとプロモーション情報promo_dfを用いてProphetによる予測を行う関数 forecast_demand_with_featuresを準備する。 引数： - demand_df: 需要の系列を表すデータフレーム。顧客Cust、製品Prod,日時Date、需要量demandの列をもつものと仮定する。 - cust_selected: 顧客のリスト - prod_selected: 製品のリスト - promo_df: プロモーション情報を入れたデータフレーム；demand_dfと同じ日時と、予測する期間に対する日時の情報を含むものとする。 - feature: promo_df内で特徴量として用いる列名のリスト - agg_period: 集約を行う期。例えば１週間を１期とするなら”1w”とする。既定値は1日（“1d”） - forecast_period: 予測を行う未来の期数。既定値は1 - cumulative: 累積需要量を用いて予測するとき True；既定値はFalse - lb: 下限値（傾向変動がロジスティック曲線の場合に有効になる．）既定値は0. - ub: 上限値（傾向変動がロジスティック曲線の場合に有効になる．）既定値は Noneで，その場合にはデータの最大値が設定される． - kwargs: その他のProphetのキーワード引数 Prophetの引数と既定値： - growth=‘linear’ :傾向変動の関数．既定値は線形．ロジスティック曲線にするに’logistic’に設定する．また傾向変動をなくすためには引数を ’flat’に設定する． - 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.80 : 不確実性の幅 - uncertainty_samples=1000 : 不確実性の幅を計算する際のサンプル数 返値： 以下の情報を列情報にもつデータフレーム - dic: 顧客と製品のタプルをキーとし、モデルと予測オブジェクトのタプルを値とした辞書 ------------------------------------------------------------------------ source ### forecast_demand_with_features ``` python def forecast_demand_with_features( demand_df, cust_selected, prod_selected, promo_df:NoneType=None, features:NoneType=None, agg_period:str='1d', forecast_periods:int=1, cumulative:bool=False, lb:float=0.0, ub:NoneType=None, kwargs:VAR_KEYWORD ): ``` *Prophetを用いた需要予測（顧客と製品に対して）：特徴ベクトルの追加* ### forecast_demand_with_features関数の使用例 ``` python promo_df = pd.read_csv(folder+"promo.csv", index_col=0) demand_df = pd.read_csv(folder+"demand_with_promo.csv") cust_selected = [ demand_df["cust"].iloc[90],demand_df["cust"].iloc[101] ] prod_selected = [ demand_df["prod"].iloc[1] ] #cust_selected = ["前橋市"] #prod_selected = ["F"] print(cust_selected, prod_selected) dic, df = forecast_demand_with_features(demand_df, cust_selected, prod_selected, promo_df, features=['promo_0', 'promo_1'], agg_period ="1d", forecast_periods =0, cumulative=False, growth="linear", daily_seasonality=False, weekly_seasonality=True, yearly_seasonality=True) if dic[ cust_selected[0], prod_selected[0] ] is not None: model, forecast = dic[ cust_selected[0], prod_selected[0] ] model.plot(forecast); else: print("too few data") ``` ['札幌市', '仙台市'] ['B'] AttributeError: 'Prophet' object has no attribute 'stan_backend' --------------------------------------------------------------------------- AttributeError Traceback (most recent call last) /var/folders/69/5y96sdc94jxf6khgc8mlmxrr0000gn/T/ipykernel_94971/15301494.py in   7 #prod_selected = ["F"]  8 print(cust_selected, prod_selected) ----> 9 dic, df = forecast_demand_with_features(demand_df, cust_selected, prod_selected, promo_df, features=['promo_0', 'promo_1'],  10 agg_period ="1d", forecast_periods =0, cumulative=False, growth="linear",  11 daily_seasonality=False, weekly_seasonality=True, yearly_seasonality=True) /var/folders/69/5y96sdc94jxf6khgc8mlmxrr0000gn/T/ipykernel_94971/3772577099.py in forecast_demand_with_features(demand_df, cust_selected, prod_selected, promo_df, features, agg_period, forecast_periods, cumulative, lb, ub, **kwargs)  45 df["floor"] = lb  46  ---> 47 model = Prophet(**kwargs)  48   49 #プロモーションの特徴ベクトルを追加 ~/Library/Caches/pypoetry/virtualenvs/scmopt-KVX_UVCf-py3.8/lib/python3.8/site-packages/fbprophet/forecaster.py in __init__(self, growth, changepoints, n_changepoints, changepoint_range, yearly_seasonality, weekly_seasonality, daily_seasonality, holidays, seasonality_mode, seasonality_prior_scale, holidays_prior_scale, changepoint_prior_scale, mcmc_samples, interval_width, uncertainty_samples, stan_backend)  139 self.fit_kwargs = {}  140 self.validate_inputs() --> 141 self._load_stan_backend(stan_backend)  142   143 def _load_stan_backend(self, stan_backend): ~/Library/Caches/pypoetry/virtualenvs/scmopt-KVX_UVCf-py3.8/lib/python3.8/site-packages/fbprophet/forecaster.py in _load_stan_backend(self, stan_backend)  152 self.stan_backend = StanBackendEnum.get_backend_class(stan_backend)()  153  --> 154 logger.debug("Loaded stan backend: %s", self.stan_backend.get_type())  155   156 def validate_inputs(self): AttributeError: 'Prophet' object has no attribute 'stan_backend' ### 集約したデータに対する予測 ``` python demand_df = pd.read_csv(folder+"aggregate_demand.csv") cust_selected = [ demand_df["cust"].iloc[3] ] prod_selected = [ demand_df["prod"].iloc[3] ] print(cust_selected, prod_selected) dic, df = forecast_demand_with_features(demand_df, cust_selected, prod_selected, features=[], agg_period ="1d", forecast_periods =100, cumulative=False ,growth="logistic", daily_seasonality=False, weekly_seasonality=True, yearly_seasonality=True) print(cust_selected[0], prod_selected[0]) if dic[ cust_selected[0], prod_selected[0] ] is not None: model, forecast = dic[ cust_selected[0], prod_selected[0] ] model.plot(forecast); else: print("too few data") ``` ['Plnt_Odawara'] ['B'] Plnt_Odawara B ![](06forecast_files/figure-commonmark/cell-24-output-2.png) ``` python len(forecast) ``` 730 ## 交差検証 ``` python #df_cv = cross_validation(model, initial='220 days', period='100 days', horizon = '100 days') df_cv = cross_validation(model, initial='200 days', period='30 days', horizon = '60 days') df_cv.head() ``` INFO:fbprophet:Making 16 forecasts with cutoffs between 2019-08-08 00:00:00 and 2020-10-31 00:00:00 WARNING:fbprophet:Seasonality has period of 365.25 days which is larger than initial window. Consider increasing initial. 0%| | 0/16 [00:00
ds yhat yhat_lower yhat_upper y cutoff
0 2019-08-09 2.806365 1.621363 3.974812 2 2019-08-08
1 2019-08-10 4.407486 3.145626 5.640692 5 2019-08-08
2 2019-08-11 1.621587 0.445880 2.891954 0 2019-08-08
3 2019-08-12 1.618451 0.417103 2.789864 0 2019-08-08
4 2019-08-13 3.644485 2.507137 4.867713 3 2019-08-08
``` python df_p = performance_metrics(df_cv, rolling_window=0.1) df_p.head() ``` INFO:fbprophet:Skipping MAPE because y close to 0
horizon mse rmse mae mdape coverage
0 6 days 1.866793 1.366306 1.128297 0.556571 0.625000
1 7 days 1.823782 1.350475 1.101126 0.452552 0.645833
2 8 days 1.747038 1.321756 1.099327 0.506288 0.677083
3 9 days 1.669610 1.292134 1.050407 0.445765 0.697917
4 10 days 1.858351 1.363214 1.065813 0.410124 0.687500
``` python plot_cross_validation_metric(df_cv, metric='mse'); ``` ![](06forecast_files/figure-commonmark/cell-29-output-1.png) ## 平均や標準偏差も計算 過去の需要系列dfをもとに、Prophetによる予測を行う関数 forecast_prophet を準備する。 比較として、（需要の定常性を仮定して）平均と標準偏差も算出する。 引数： - demand_df: 需要の系列を表すデータフレーム。顧客Cust、製品Prod,日時Date、需要量demandの列をもつものと仮定する。 - agg_period: 集約を行う期。例えば１週間を１期とするなら”1w”とする。既定値は1日（“1d”） - forecast_period: 予測を行う未来の期数。既定値は1 返値： 以下の情報を列情報にもつデータフレーム - mean: 需要の平均 - std : 需要の標準偏差 - Prophet\[0..forecast_period-1\]: Prophetによる予測から計算された未来のforecast_periods 期分の最大在庫量 ------------------------------------------------------------------------ source ### forecast_prophet ``` python def forecast_prophet( df, agg_period:str='1d', forecast_periods:int=1 ): ``` *Prophetを用いた需要予測* ------------------------------------------------------------------------ source ### forecast_prophet_cust_prod ``` python def forecast_prophet_cust_prod( demand_df, cust_selected, prod_selected, agg_period:str='1d', forecast_periods:int=1, cumulative:bool=False ): ``` *Prophetを用いた需要予測（顧客と製品に対して）* ``` python # ProphetのWarningがたくさん出る場合には、以下を生かす。 import warnings warnings.simplefilter('ignore') demand_df = pd.read_csv(folder+"demand.csv") cust_selected = [ demand_df["cust"].iloc[0], demand_df["cust"].iloc[-1] ] prod_selected = [ demand_df["prod"].iloc[0], demand_df["prod"].iloc[-1] ] print(cust_selected, prod_selected) forecast_df = forecast_prophet_cust_prod(demand_df, cust_selected, prod_selected, agg_period ="1w", forecast_periods =10) print(forecast_df.head()) ``` ------------------------------------------------------------------------ source ### forecast_prophet_cust ``` python def forecast_prophet_cust( demand_df, cust_selected, agg_period:str='1d', forecast_periods:int=1 ): ``` *各顧客に対する（すべての製品の）合計需要の予測* ``` python demand_df = pd.read_csv(folder+"demand.csv") cust_selected = [ demand_df["cust"].iloc[0], demand_df["cust"].iloc[-1] ] forecast_df = forecast_prophet_cust(demand_df, cust_selected, agg_period ="1w", forecast_periods =10) print(forecast_df.head()) ``` Prophet Error while computing for 札幌市 Prophet Error while computing for 那覇市 Empty DataFrame Columns: [cust, mean, std, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9] Index: [] ------------------------------------------------------------------------ source ### forecast_prophet_prod ``` python def forecast_prophet_prod( demand_df, prod_selected, agg_period:str='1d', forecast_periods:int=1 ): ``` *各製品に対する（すべての顧客の）合計需要の予測* ``` python demand_df = pd.read_csv(folder+"demand.csv") prod_selected = [ demand_df["prod"].iloc[0], demand_df["prod"].iloc[-1] ] forecast_df = forecast_prophet_prod(demand_df, prod_selected, agg_period ="1w", forecast_periods =10) print(forecast_df.head()) ``` Prophet Error while computing for A Prophet Error while computing for J Empty DataFrame Columns: [prod, mean, std, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9] Index: [] ## 深層学習ライブラリ fastai を用いた需要予測 ### 未来の時系列データフレームを生成する関数 make_future_dataframe 過去の時系列データから未来のデータを生成する関数 - history_dates: 過去の時系列データの日付を表す - finish: 終了日 - periods: 生成する未来の期数 - freq: 生成する期の単位（頻度）を表す文字列; たとえば，1日単位の場合には”1d”とする． ------------------------------------------------------------------------ source ### make_future_dataframe ``` python def make_future_dataframe( history_dates, periods, freq:str='1d', include_history:bool=True ): ``` *未来の時系列データフレームを生成する関数* #### make_future_dataframe関数の使用例 ``` python #製品／顧客を切り出したデータフレームの保存 promo_df = pd.read_csv(folder+"promo.csv", index_col=0) demand_df = pd.read_csv(folder+"demand_with_promo.csv",index_col=0) customer = "仙台市" product = "A" agg_period ="1d" try: demand_df.reset_index(inplace=True) except ValueError: pass demand_df["date"] = pd.to_datetime(demand_df["date"]) demand_df.set_index("date", inplace=True) demand_grouped = demand_df.groupby( ["cust", "prod"]).resample(agg_period)["demand"].sum() if (customer,product) in demand_grouped: df = pd.DataFrame(demand_grouped[customer, product]) # df.reset_index(inplace=True) # if promo_df is not None: # promo_df["date"] = pd.to_datetime(promo_df["date"]) # df = pd.merge(df, promo_df,on="date", how="left") #df.to_csv("maebashi-f.csv") df.head() ```
demand
date
2019-01-01 13
2019-01-02 19
2019-01-03 17
2019-01-04 12
2019-01-05 19
``` python future = make_future_dataframe(history_dates=df.index, periods=10, freq="1w", include_history=False) future.head() ```
date
0 2021-01-03
1 2021-01-10
2 2021-01-17
3 2021-01-24
4 2021-01-31
### Excelのテンプレートを生成する関数 make_excel_forecast Excelインターフェイスのテンプレートを生成する関数． 引数： - start: 開始日 - finish: 終了日 - period: 期を構成する単位期間の数；既定値は 1 - period_unit: 期の単位 （時，日，週，月から選択）； 既定値は日； periodとあわせて期間の生成に用いる． たとえば，既定値だと１日が1期となる． - num_features: （休日を除く）特徴の数 - holiday: 日本のカレンダーに基づく休日を表す特徴を追加するときTrue 返値： - 需要予測用のExcelのWorkbook ------------------------------------------------------------------------ source ### make_excel_forecast ``` python def make_excel_forecast( start, finish, period:int=1, period_unit:str='日', num_features:int=0, holiday:bool=True ): ``` *Call self as a function.* #### make_excel_forecast関数の使用例 ``` python start = '2015-01-05' #開始日 finish = '2019-1-30' period=1 period_unit="週" wb = make_excel_forecast(start, finish, period, period_unit, num_features = 0, holiday=False) wb.save("forecast-templete.xlsx") ``` ### Excelファイルを読み込んで予測を行う関数 forecast_excel ExcelのWorkbookを与えると予測を記入して返す関数 引数： - wb: 日付，需要，特徴のデータを入れたExcel Workbook - valid_start: 検証開始日 - predict_start: 予測開始日 返値： - wb: 予測結果と誤差の列を追加したExcel Workbook - best: 最良の予測をしたモデルオブジェクト - result_df: 予測結果のデータフレーム - forecast: 最良の予測モデルを用いた予測結果のデータフレーム - demand_df: 予測結果から需要を抽出し，誤差列を追加したデータフレーム - valid_df: demand_dfから検証開始から予測開始の直前までを抜き出したデータフレーム ------------------------------------------------------------------------ source ### forecast_excel ``` python def forecast_excel( wb, valid_start, predict_start ): ``` *Call self as a function.* #### forecast_excel関数の使用例 ``` python # wb = load_workbook("forecast-ex2.xlsx", data_only=True) # valid_start = "2019-09-01" # predict_start = "2020-12-30" wb = load_workbook("forecast-ex1.xlsx", data_only=True) valid_start = "2021-02-24" predict_start = "2021-3-03" #wb = load_workbook(folder+"forecast-walmart.xlsx", data_only=True) #valid_start = "2011-8-19" #predict_start = "2011-12-30" #wb = load_workbook("forecast-ex5.xlsx", data_only=True) #valid_start = "2020-10-19" #predict_start = "2020-12-30" wb, best, result_df, forecast, demand_df, valid_df = forecast_excel(wb, valid_start, predict_start) ```
Model MAE MSE RMSE R2 RMSLE MAPE TT (Sec)
dummy Dummy Regressor 10.5139 228.3212 15.1103 -0.1966 0.1213 0.0811 0.0000
[DummyRegressor(constant=None, quantile=None, strategy='mean')] ### 予測に基づく動的発注量シートの生成関数 dynamic_inv_excel 引数： - wb: 在庫シミュレーションシートを入れるExcel Workbook - valid_start: 検証開始日 - predict_start: 予測開始日 - valid_df: demand_dfから検証開始から予測開始の直前までを抜き出したデータフレーム - demand_df: 予測結果から需要を抽出し，誤差列を追加したデータフレーム - h: 在庫費用（定数かリスト）; 関数内で需要と同じ長さのNumPy配列に変換される． - fc: 固定費用（定数かリスト）; 関数内で需要と同じ長さのNumPy配列に変換される． - vc: 変動費用（定数かリスト）; 関数内で需要と同じ長さのNumPy配列に変換される；既定値は0. - z: 安全在庫係数；初期安全在庫量を決めるために使われる定数（既定値は2.33） 返値： - wb: 在庫シミュレーションシートを追加したExcel Workbook ------------------------------------------------------------------------ source ### dynamic_inv_excel ``` python def dynamic_inv_excel( wb, valid_start, predict_start, valid_df, demand_df, h, fc, vc:float=0.0, z:float=2.33 ): ``` *Call self as a function.* #### dynamic_inv_excel関数の使用例 ``` python h = 1 fc = 100. vc = 0. z = 1.65 wb = dynamic_inv_excel(wb, valid_start, predict_start, valid_df, demand_df, h, fc, vc, z) wb.save("forecast-ex5-simu.xlsx") ``` 2995.525076709324 [44.36038177 0. 0. 65.51441261 0. 0. 0. 0. 61.73584662 0. 0. 49.47059588 0. 0. 0. 0. 57.16618104 0. 0. 0. 0. 0. 50.6413246 0. 0. 0. 29.14120463 0. 0. 0. 54.8727735 0. 0. 0. 0. 0. 60.59831692 0. 0. 0. 0. 0. 0. 51.70768912 0. 0. 0. 0. 0. 0. 53.16539093 0. 0. 0. 0. 0. 0. 48.67772177 0. 0. 0. 0. 0. 0. 52.11188276 0. 0. 0. 0. 0. 0. 0. 62.27572896 0. 0. 0. 0. 0. 53.92021115 0. 0. 42.08768594 0. 0. ] ### 顧客と製品を与えると、訓練用データフレームと予測用の未来のデータフレームを生成する関数 需要がない日を自動的に追加し，その需要を0と補完する． 引数： - demand_df: 需要データフレーム - customer: 予測を行いたい顧客 - product: 予測を行うたい製品 - promo_df: 追加特徴（プロモーション）情報を含んだデータフレーム（オプション）；予測期間に対しても追加情報を準備する必要がある． 予測期間内のデータが空 (NaN)の場合には，訓練データにも空 (NaN)のものがある必要がある． - agg_period: 期間を表す文字列（例えば1日のときは”1d”） - forecast_periods: 予測したい期間 返値： - df: 訓練用データフレーム - future: 訓練用と予測用の未来のデータを入れたデータフレーム ------------------------------------------------------------------------ source ### make_forecast_df ``` python def make_forecast_df( demand_df, customer, product, promo_df:NoneType=None, agg_period:str='1d', forecast_periods:int=1 ): ``` *顧客と製品を与えると、訓練用データフレームと予測用の未来のデータフレームを生成する関数* ### make_forecast_df関数の使用例 ``` python c = "仙台市" p = "A" agg_period ="1d" df, future = make_forecast_df(demand_df, customer=c, product=p,promo_df= promo_df, agg_period=agg_period, forecast_periods = 10) df.head() ```
demand promo_0 promo_1 Year Month Week Day Dayofweek Dayofyear Is_month_end Is_month_start Is_quarter_end Is_quarter_start Is_year_end Is_year_start Elapsed
0 13 0 0 2019 1 1 1 1 1 False True False True False True 1.546301e+09
1 19 1 0 2019 1 1 2 2 2 False False False False False False 1.546387e+09
2 17 0 0 2019 1 1 3 3 3 False False False False False False 1.546474e+09
3 12 0 0 2019 1 1 4 4 4 False False False False False False 1.546560e+09
4 19 0 0 2019 1 1 5 5 5 False False False False False False 1.546646e+09
### 実データでの分析 ### 全てのデータを用いて深層学習による予測をするためのデータフレームの生成関数 データ数が不足している場合には、顧客や製品もカテゴリー変数として扱うことによって十分なデータを確保することができる。 引数： - demand_df: 需要データフレーム - promo_df: 追加特徴（プロモーション）情報を含んだデータフレーム（オプション）；予測期間に対しても追加情報を準備する必要がある． 予測期間内のデータが空 (NaN)の場合には，訓練データにも空 (NaN)のものがある必要がある． - agg_period: 期間を表す文字列（例えば1日のときは”1d”） 返値： - demand_df: 訓練用データフレーム ------------------------------------------------------------------------ source ### make_forecast_df_all ``` python def make_forecast_df_all( demand_df, promo_df:NoneType=None, agg_period:str='1d', forecast_periods:int=1 ): ``` *全データに対する訓練用データフレームを生成する関数* ``` python demand_df = pd.read_csv(folder+"demand_with_promo.csv", index_col=0) all_df, future = make_forecast_df_all(demand_df, promo_df, agg_period="1d") all_df.head() ```
cust prod Week promo_0 promo_1 demand Year Month Day Dayofweek Dayofyear Is_month_end Is_month_start Is_quarter_end Is_quarter_start Is_year_end Is_year_start Elapsed
0 仙台市 A 1 0 0 13 2019 1 1 1 1 False True False True False True 1546300800
1 仙台市 A 1 1 0 19 2019 1 2 2 2 False False False False False False 1546387200
2 仙台市 A 1 0 0 17 2019 1 3 3 3 False False False False False False 1546473600
3 仙台市 A 1 0 0 12 2019 1 4 4 4 False False False False False False 1546560000
4 仙台市 A 1 0 0 19 2019 1 5 5 5 False False False False False False 1546646400
### 検証期間を設定するための関数 find_horizon ------------------------------------------------------------------------ source ### find_horizon ``` python def find_horizon( demand_df, ratio:float=0.9 ): ``` *Call self as a function.* ``` python demand_df = pd.read_csv(folder+"demand_with_promo.csv") find_horizon(demand_df, ratio=0.9) ``` Timestamp('2020-10-19 00:00:00') ### 訓練データと検証データのインデックスを準備する関数 prepare_index 引数： - df: 日付情報（Year列、Month列、Day列）を含んだデータフレーム - horizon: 検証データを作成するときの閾値の日付；この日付より後のデータは検証データになる。Noneのときには最後の100日を検証データとする。 返値： - train_idx: 訓練データのインデックス - valid_idx: 検証データのインデックス ------------------------------------------------------------------------ source ### prepare_index ``` python def prepare_index( df, horizon:NoneType=None ): ``` *Call self as a function.* ``` python #train_idx, valid_idx =prepare_index(df, horizon="2020/11/21") #len(train_idx), len(valid_idx) ``` ## 訓練 需要を連続量とする場合と離散量とする場合で異なる予測を行う。連続量の場合には回帰問題となり、離散量の場合には分類問題になる。 ### 学習器を作成する関数 make_learn 引数： - df: 日付 (date)と需要 (deman)の列を含んだデータフレーム - horizon: 検証データを作成するときの閾値の日付；この日付より後のデータは検証データになる。 - classification: 分類問題としてモデルを作成する場合にはTrue、回帰問題としてモデルを作成する場合には Falseを入れる。 返値： - learn: fastaiの学習器 ------------------------------------------------------------------------ source ### make_learn ``` python def make_learn( df, horizon:NoneType=None, classification:bool=False ): ``` *深層学習を用いた需要予測用のデータを生成* ### make_learn関数の使用例 ``` python learn = make_learn(df, horizon="2020-10-01", classification=False) #learn = make_learn(all_df, horizon="2016-10-01", classification=False) #learn = make_learn(all_df, horizon="2020-10-01", classification=False) #df.drop("level_0", axis=1, inplace=True) #learn = make_learn(df, horizon="2016-12-01", classification=False) ``` 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. ## 適正な学習率を見つけるための関数 lr_find ``` python lr_min= learn.lr_find() print(lr_min) ``` SuggestedLRs(valley=0.002511886414140463) ### lr_findの結果をPlotlyで描画する関数 plot_lr_find ------------------------------------------------------------------------ source ### plot_lr_find ``` python def plot_lr_find( learn ): ``` *Call self as a function.* ``` python fig = plot_lr_find(learn) plotly.offline.plot(fig); ``` ## 訓練を行う関数 train 引数： - learn: 学習器 - kwargs: 訓練を行う関数fit_one_cycleの任意のキーワード引数 > https://docs.fast.ai/callback.schedule#Learner.fit_one_cycle 参照 ------------------------------------------------------------------------ source ### train ``` python def train( learn, kwargs:VAR_KEYWORD ): ``` *fit_one_cycle法を用いて訓練を行う関数* ``` python train(learn, n_epoch = 30, lr_max =0.01, wd=2., cbs=ShowGraphCallback()) ```
epoch train_loss valid_loss _rmse mae msle time
0 3446.414062 7105.580566 84.294609 66.609749 4.805219 00:03
1 2210.760254 1613.519287 40.168633 26.273926 2.868655 00:03
2 1863.013672 656.649414 25.625174 13.304914 0.741323 00:03
3 1639.386475 976.197571 31.244160 14.256841 0.766732 00:03
4 1687.551270 648.089172 25.457594 12.055953 0.347193 00:03
5 1956.444458 1265.262207 35.570530 14.617914 0.479162 00:03
6 1698.531738 466.392212 21.596113 11.547545 1.536979 00:03
7 2314.256104 1495.722778 38.674576 16.015421 0.508270 00:03
8 1567.057861 456.358246 21.362545 17.110693 1.679504 00:03
9 1580.399536 583.629700 24.158430 16.468655 1.654195 00:03
10 1677.902344 658.983154 25.670668 14.227415 0.733416 00:03
11 1710.951782 3192.795898 56.504829 31.815025 1.787573 00:03
12 1314.445557 1115.158813 33.393997 17.580637 2.050577 00:03
13 1093.038452 407.138885 20.177683 11.029508 0.769333 00:03
14 954.789368 1112.169678 33.349209 16.030300 1.018057 00:03
15 1350.003052 617.209290 24.843697 12.732859 0.934637 00:03
16 1671.424316 484.320282 22.007277 11.749466 0.757320 00:03
17 800.749023 998.935669 31.605944 15.915786 0.790849 00:03
18 749.288147 683.855164 26.150621 13.572218 1.163126 00:03
19 858.682617 876.759644 29.610126 13.271565 0.624707 00:03
20 594.039185 418.807831 20.464794 8.975345 0.675335 00:03
21 476.235504 276.303314 16.622375 9.983256 0.377772 00:03
22 443.672668 179.968506 13.415236 6.672860 0.587208 00:03
23 423.546722 132.249161 11.499963 5.285238 0.299202 00:03
24 334.764893 209.028351 14.457813 6.568413 0.476533 00:03
25 208.162155 72.868652 8.536314 4.375165 0.419365 00:03
26 155.964676 41.782406 6.463932 3.183359 0.320197 00:04
27 119.822289 44.234440 6.650898 3.861167 0.297581 00:03
28 93.038322 25.414803 5.041310 2.986454 0.235738 00:03
29 86.351746 15.100967 3.885997 2.207297 0.232204 00:03
![](06forecast_files/figure-commonmark/cell-74-output-2.png) ### Plotlyで損出関数を描画する関数 plot_loss 引数： - learn: 学習器 - valid: テストデータも描画するときにTrue（既定値） 返り値： - fig: 損出関数のグラフ ------------------------------------------------------------------------ source ### plot_loss ``` python def plot_loss( learn, valid:bool=True ): ``` *Call self as a function.* ``` python fig = plot_loss(learn) plotly.offline.plot(fig); ``` ![](06forecast_files/figure-commonmark/cell-77-output-1.png) ## 深層学習を用いて予測を行う関数 predict_using_dl 引数： - learn:学習器 - df: 日付 (date)と需要 (deman)の列を含んだデータフレーム - future: 訓練用と予測用の未来のデータを入れたデータフレーム - classification: 分類問題のとき True、それ以外のとき False 返値： - fig: 予測結果を入れた図オブジェクト - predict: 予測結果を入れた配列 - error: 予測の2乗誤差 - r2: 決定係数 ------------------------------------------------------------------------ source ### predict_using_dl ``` python def predict_using_dl( learn, df, future, classification:bool=False ): ``` *深層学習を用いた予測* ### predict_using_dl関数の使用例 ``` python fig, yhat, error, r2 = predict_using_dl(learn, df, future, classification=False) print("error=", error, r2) plotly.offline.plot(fig); ``` error= 103425.1953125 0.023292243480682373 ![](06forecast_files/figure-commonmark/cell-80-output-1.png) ### 全てのデータで訓練した学習器による予測関数 predict_using_all_dl 引数： - learn: 学習器 - df: 該当する顧客と製品の需要データを入れたデータフレーム - furure: 予測したい全期間の需要以外のデータを入れたデータフレーム - customer: 顧客名 - product: 製品名 - classification: 分類問題のとき True、それ以外のとき False 返値： - fig: 予測結果を入れた図オブジェクト - predict: 予測結果を入れた配列 - error: 予測の2乗誤差 ------------------------------------------------------------------------ source ### predict_using_all_dl ``` python def predict_using_all_dl( learn, df, future, customer, product, classification:bool=True ): ``` *深層学習を用いた予測（全部で訓練）* ### predict_using_all_dl関数の使用例 ``` python # demand_df = pd.read_csv(folder+"demand_with_promo.csv", index_col=0) c = "前橋市" p = "B" df, future = make_forecast_df(demand_df, customer=c, product=p, promo_df= promo_df, agg_period="1d", forecast_periods =100) fig, predict, error, r2 = predict_using_all_dl(learn, df, future, customer=c, product = p, classification=False) print("error=", error, r2) plotly.offline.plot(fig); ``` error= 15.100969314575195 0.9979501962661743 ``` python #hide_iinput Image("../figure/predict_using_all_dl.png") ``` ![](06forecast_files/figure-commonmark/cell-84-output-1.png) ### 埋め込み層の描画関数 draw_embedding 引数： - learn: 学習器（訓練済み） - cust: 顧客を描画するときTrue, 製品を描画するとき False - pca: 主成分分析で次元削減するときTrue、TSNEでするとき False 返値： - fig: Plotlyの図オブジェクト ------------------------------------------------------------------------ source ### draw_embedding ``` python def draw_embedding( learn, cust:bool=True, pca:bool=True ): ``` *draw embedding* ### draw_embedding関数の使用例 ``` python fig = draw_embedding(learn, cust=False, pca=True) plotly.offline.plot(fig); ``` ``` python fig = draw_embedding(learn, cust=True, pca=True) plotly.offline.plot(fig); ``` ![](06forecast_files/figure-commonmark/cell-89-output-1.png) ## scikit learnのランダム森によって予測をする関数 random_forest 引数： - df: 日付 (date)と需要 (deman)の列を含んだデータフレーム - horizon: 検証データを作成するときの閾値の日付；この日付より後のデータは検証データになる。 - feature: 特徴ベクトルとして使用する列名のリスト； Noneのときには全ての列とする。 - classification: 分類問題のときTrue、回帰問題のときFalse 返値： - forest: ランダム森のモデルオブジェクト - importance: permutation_importanceによる特徴ベクトルの重要度の情報 ------------------------------------------------------------------------ source ### random_forest ``` python def random_forest( df, horizon:NoneType=None, feature:NoneType=None, classification:bool=False ): ``` *ランダム森による需要予測* ### random_forest関数の呼び出し例 ``` python # list(imp[imp.importance>0.01].features) #feature = ['promo_1', 'Week', 'Month', 'Dayofyear', 'promo_0', 'Dayofweek'] #feature = list(df.columns[1:]) feature = None forest, importance = random_forest(df, horizon="2020-10-01", feature=feature, classification=False) ``` ### 重要度のデータフレームと図を返す関数 show_importance 引数： - importance: 重要度（順列重要度を random_forest関数で計算したときの返値として得られる．） - feature: 特徴ベクトルの名前のリスト - df: データフレーム；特徴ベクトルの元の列名を得るためだけに使用する． 返値： imp_df: 重要度を入れたデータフレーム fig: 重要度の棒グラフ ------------------------------------------------------------------------ source ### show_importance ``` python def show_importance( importance, feature, df ): ``` *Call self as a function.* ``` python feature = None imp_df, fig = show_importance(importance, feature=feature, df=df) plotly.offline.plot(fig); ``` ![](06forecast_files/figure-commonmark/cell-94-output-1.png) ### ランダム森による予測 predict_using_rf 引数： - df: 日付 (date)と需要 (deman)の列を含んだデータフレーム - future: 訓練用と予測用の未来のデータを入れたデータフレーム - forest: ランダム森のモデルオブジェクト - horizon: 予測期間 返値： - fig: 予測結果を入れた図オブジェクト - predict: 予測結果を入れた配列 - error: 予測の2乗誤差 - r2: 決定係数 ------------------------------------------------------------------------ source ### predict_using_rf ``` python def predict_using_rf( df, future, forest, horizon:NoneType=None, feature:NoneType=None ): ``` *ランダム木を用いた予測* ``` python fig, yhat_rf, error_rf, r2 = predict_using_rf(df, future, forest, horizon="2020-10-01", feature=feature) plotly.offline.plot(fig); print("error", error_rf, r2) ``` error 1.0633877777777778 0.5972015993265993 error 1.1132505494505496 ![](06forecast_files/figure-commonmark/cell-97-output-2.png) ## AutoMLで予測をする関数 automl AutoMLパッケージPyCaretによる予測 引数： - df: 日付 (date)と需要 (demand)の列を含んだデータフレーム - horizon: 検証データを作成するときの閾値の日付；この日付より後のデータは検証データになる。Noneの場合には100日前からが検証データになる． - num_bests: 最良モデルの選択数 - sort: 良い順に並べるときの評価尺度名（既定値は”R2”） - model_name: 返値： - best: 最良の予測をしたモデルオブジェクト - result_df: 結果のデータフレーム ------------------------------------------------------------------------ source ### automl ``` python def automl( df, horizon:NoneType=None, num_bests:int=5, sort:str='R2', model_name:NoneType=None ): ``` *Prediction using pycaret regression models* ### automl関数の使用例 ``` python horizon="2020/10/01" best, result_df = automl(df, horizon, 5, model_name=None) ```
Model MAE MSE RMSE R2 RMSLE MAPE TT (Sec)
et Extra Trees Regressor 1.6620 4.9674 2.2288 0.9262 0.3713 0.2443 0.1500
rf Random Forest Regressor 2.3624 7.3216 2.7058 0.8913 0.4690 0.3809 0.1900
dt Decision Tree Regressor 2.4286 10.6044 3.2564 0.8425 0.5424 0.3973 0.0100
gbr Gradient Boosting Regressor 3.3336 14.5483 3.8142 0.7840 0.4328 0.5026 0.0800
ada AdaBoost Regressor 5.7564 40.3512 6.3523 0.4008 0.8888 1.1191 0.1200
br Bayesian Ridge 6.5873 61.8033 7.8615 0.0823 1.0977 1.7060 0.0100
ridge Ridge Regression 6.6831 63.5019 7.9688 0.0571 1.1030 1.7497 0.0100
huber Huber Regressor 6.4650 63.8005 7.9875 0.0527 1.0793 1.4744 0.0200
kr Kernel Ridge 6.5818 65.3509 8.0840 0.0296 1.0998 1.5512 0.1100
en Elastic Net 6.7439 67.1895 8.1969 0.0023 1.1457 1.7762 0.0100
lasso Lasso Regression 6.7995 68.3639 8.2682 -0.0151 1.1519 1.7897 0.0100
tr TheilSen Regressor 7.0673 69.5302 8.3385 -0.0324 1.1465 1.8813 1.6900
lr Linear Regression 7.0492 69.5646 8.3405 -0.0329 1.1389 1.8974 0.0000
lar Least Angle Regression 7.0492 69.5646 8.3405 -0.0329 1.1389 1.8974 0.0100
ard Automatic Relevance Determination 7.0675 69.8744 8.3591 -0.0375 1.1409 1.8141 0.0100
omp Orthogonal Matching Pursuit 7.1926 71.5753 8.4602 -0.0628 1.1591 1.8336 0.0000
svm Support Vector Regression 7.0575 72.9081 8.5386 -0.0826 1.1784 1.8887 0.0200
llar Lasso Least Angle Regression 7.2027 75.3612 8.6811 -0.1190 1.1951 1.9659 0.0000
dummy Dummy Regressor 7.2027 75.3612 8.6811 -0.1190 1.1951 1.9659 0.0000
ransac Random Sample Consensus 8.5524 105.8737 10.2895 -0.5721 1.2047 2.2048 0.0500
knn K Neighbors Regressor 9.4000 117.2127 10.8265 -0.7404 1.3451 2.7281 0.0000
par Passive Aggressive Regressor 8.5968 122.1332 11.0514 -0.8135 1.2091 0.8227 0.0000
mlp MLP Regressor 7927.3291 62842730.7068 7927.3407 -933111.6422 6.9995 1633.9990 0.2000
``` python #予測 result = predict_model(best[0], df) result.head() ```
demand promo_0 promo_1 Year Month Week Day Dayofweek Dayofyear Is_month_end Is_month_start Is_quarter_end Is_quarter_start Is_year_end Is_year_start Elapsed Label
0 13 0 0 2019 1 1 1 1 1 False True False True False True 1.546301e+09 13.0
1 19 1 0 2019 1 1 2 2 2 False False False False False False 1.546387e+09 19.0
2 17 0 0 2019 1 1 3 3 3 False False False False False False 1.546474e+09 17.0
3 12 0 0 2019 1 1 4 4 4 False False False False False False 1.546560e+09 12.0
4 19 0 0 2019 1 1 5 5 5 False False False False False False 1.546646e+09 19.0
``` python #plot_model(best[1], display_format="streamlit") ``` ### 複数の製品・顧客に対して予測を行う関数（のサブ関数） automl_all_sub サプライ・チェイン最適化のためには，すべての顧客と製品に対して一度に需要予測を行う必要がある．しかし，通常は顧客数，製品数ともに非常に多い場合には，計算時間が問題になる．顧客や製品に対するABC分析を事前に行い，A製品（A顧客）に対してだけ，需要予測を行うという経験的手法が従来使われていたが，ここではデータに基づいた手法を考える． 需要予測は，顧客・製品ごとに時系列データ（にプロモーション情報を付加したもの）を抽出し，AutoMLで複数の機械学習・統計分析モデルを適用して，最も良いものを使うことにする．評価尺度としては決定係数を使う．決定係数が負のものは平均を用いた単純な予測より当たらないことを意味する． 顧客・製品ごとに需要の合計量を被増加順に並べ，大きいものから順に予測を行っていく．一般に，zip分布にしたがう需要を想定すると，需要の合計が小さい場合には，どのような予測をしても当たらなくなる．複数の手法を試し，最も良い予測の決定係数が負のものが続く場合には，需要が小さすぎてこれ以上予測をしても無駄だと判断できる．したがって，決定係数が負のものが適当な数（たとえば10）より大きくなったら，予測をやめて，顧客・製品ごとの需要の平均値と標準偏差だけで管理するものとする． 引数： - demand_df: 顧客 (cust)・製品 (prod) ごとの日付 (date)と需要 (deman)の列を含んだデータフレーム - promo_df: 追加特徴（プロモーション）情報を含んだデータフレーム（オプション）；予測期間に対しても追加情報を準備する必要がある． 予測期間内のデータが空 (NaN)の場合には，訓練データにも空 (NaN)のものがある必要がある． - cust_list: 予測を行いたい顧客のリスト - prod_list 予測を行うたい製品のリスト - horizon: 検証データを作成するときの閾値の日付；この日付より後のデータは検証データになる． - agg_period: 期間を表す文字列（例えば1日のときは”1d”） - forecast_periods: 予測したい期間（既定値は1） - threshold: 需要の予測誤差の幅を決めるためのパラメータ；外れる確率がthresholdになるように下限と上限を決める． - max_neg_r2: 予測打ち切りのためのパラメータ；この回数だけ R2 が負になったら，それ以上予測する必要はないので， 終了する．　ただし，cust_list と prod_listは需要量の降順になっているものと仮定する． 返値： - summary: 予測結果をまとめたデータフレーム；以下の列をもつ． - cust: 顧客名 - prod: 製品名 - model: 最も良い予測モデル名 - rmse: 平均２乗誤差の平方根 - r2: 決定係数 - forecast: 未来の予測値 - lb: 予測誤差がlb以下になる確率がthresholdになるように設定 - ub: 予測誤差がub以上になる確率がthresholdになるように設定 ------------------------------------------------------------------------ source ### automl_all_sub ``` python def automl_all_sub( demand_df, promo_df, cust_list, prod_list, horizon, agg_period:str='1d', forecast_periods:int=1, threshold:float=0.05, max_neg_r2:int=10 ): ``` *Call self as a function.* ### 複数の製品・顧客に対して予測を行う関数 automl_all model 列には，最も良い（決定係数が1に近い）モデル名が保管されている．これを用いて，未来の予測をしたい場合には，そのモデルだけで予測をすれば良い． 予測値と実際の需要量が大きく外れた場合には，原因を探ったり，再予測を行ったりする必要がある．それには，予測誤差の下限と上限（lb,ub列）を用いれば良い．実際の運用ができるような，自動処理を行うシステムは，MLOpsと呼ばれる． ここで開発したものは，MLOpsの基礎として使うことができる． 引数： - demand_df: 顧客 (cust)・製品 (prod) ごとの日付 (date)と需要 (deman)の列を含んだデータフレーム - promo_df: 追加特徴（プロモーション）情報を含んだデータフレーム（オプション）；予測期間に対しても追加情報を準備する必要がある． 予測期間内のデータが空 (NaN)の場合には，訓練データにも空 (NaN)のものがある必要がある． - horizon: 検証データを作成するときの閾値の日付；この日付より後のデータは検証データになる． - agg_period: 期間を表す文字列（例えば1日のときは”1d”） - forecast_periods: 予測したい期間（既定値は1） - threshold: 需要の予測誤差の幅を決めるためのパラメータ；外れる確率がthresholdになるように下限と上限を決める． - max_neg_r2: 予測打ち切りのためのパラメータ；この回数だけ *R*2 が負になったら，それ以上予測する必要はないので， 終了する．　ただし，cust_list と prod_listは需要量の降順になっているものと仮定する． - abc_threshold: ABC分析のための閾値のリスト（合計が1.0になっている必要がある．） 返値： - summary: 予測結果をまとめたデータフレーム；以下の列をもつ．（ただし予測を打ち切った行以降は model以降の列はNaNになる．） - cust: 顧客名 - prod: 製品名 - mean: 平均需要量 - std: 需要の標準偏差 - rank: 需要の降順ランク - abc: 需要量のABC分析の結果 - model: 最も良い予測モデル名 - rmse: 平均２乗誤差の平方根 - r2: 決定係数 - forecast: 未来の予測値 - lb: 予測誤差がlb以下になる確率がthresholdになるように設定 - ub: 予測誤差がub以上になる確率がthresholdになるように設定 ------------------------------------------------------------------------ source ### automl_all ``` python def automl_all( demand_df, promo_df, horizon, agg_period:str='1d', forecast_periods:int=1, threshold:float=0.05, max_neg_r2:int=10, abc_threshold:list=[0.6, 0.2, 0.2] ): ``` *Call self as a function.* ### automl_all関数の使用例 ``` python #summary = automl_all(demand_df, promo_df, horizon="2020/10/19", agg_period="1d", forecast_periods =10, threshold= 0.05, max_neg_r2 = 1) ``` ``` python summary.r2.hist(); ``` ![](06forecast_files/figure-commonmark/cell-108-output-1.png) ``` python summary.to_csv("summary.csv") #summary = pd.read_csv("summary.csv", index_col=0) #summary.sort_values(by="r2",inplace=True) summary.r2.hist(); ``` ![](06forecast_files/figure-commonmark/cell-109-output-1.png) ### すべての顧客・需要の予測結果を用いて，予測を行い図を返す関数 predict_using_automl_all_fig 引数： - demand_df: 顧客 (cust)・製品 (prod) ごとの日付 (date)と需要 (deman)の列を含んだデータフレーム - promo_df: 追加特徴（プロモーション）情報を含んだデータフレーム（オプション）；予測期間に対しても追加情報を準備する必要がある． 予測期間内のデータが空 (NaN)の場合には，訓練データにも空 (NaN)のものがある必要がある． - summary: 予測結果をまとめたデータフレーム（automl_allの返値） - customer: 顧客名 - product: 製品名 - agg_period: 期間を表す文字列（例えば1日のときは”1d”） - forecast_periods: 予測したい期間（既定値は1） 返値： - forecast: 未来の予測値（リスト） ------------------------------------------------------------------------ source ### predict_using_automl_all_fig ``` python def predict_using_automl_all_fig( demand_df, promo_df, summary, customer, product, agg_period:str='1d', forecast_periods:int=10 ): ``` *Call self as a function.* ``` python agg_df,category = abc.abc_analysis_all(demand_df,[0.6, 0.2, 0.2]) agg_df.reset_index(inplace=True) count = 0 for row in agg_df.itertuples(): fig, error, result = predict_using_automl_all_fig(demand_df, None, summary, row.cust, row.prod, agg_period='1d', forecast_periods=10) plotly.offline.plot(fig) count +=1 if count>=10: break ```
MAE MSE RMSE R2 RMSLE MAPE
0 245.907394 123375.59375 351.248596 0.1176 2.8776 2.9054
Model MAE MSE RMSE R2 RMSLE MAPE
0 Ridge Regression 245.9074 123375.5938 351.2486 0.1176 2.8776 2.9054
### すべての顧客・需要の予測結果を用いて，予測を行う関数 predict_using_automl_all 引数： - demand_df: 顧客 (cust)・製品 (prod) ごとの日付 (date)と需要 (deman)の列を含んだデータフレーム - promo_df: 追加特徴（プロモーション）情報を含んだデータフレーム（オプション）；予測期間に対しても追加情報を準備する必要がある． 予測期間内のデータが空 (NaN)の場合には，訓練データにも空 (NaN)のものがある必要がある． - summary: 予測結果をまとめたデータフレーム（automl_allの返値） - customer: 顧客名 - product: 製品名 - agg_period: 期間を表す文字列（例えば1日のときは”1d”） - forecast_periods: 予測したい期間（既定値は1） 返値： - forecast: 未来の予測値（リスト） ------------------------------------------------------------------------ source ### predict_using_automl_all ``` python def predict_using_automl_all( demand_df, promo_df, summary, customer, product, agg_period:str='1d', forecast_periods:int=10 ): ``` *Call self as a function.* ### predict_using_automl_all関数の使用例 ``` python # customer = "長崎市" # product = "D" customer ="熊本市" product ="H" predict_using_automl_all(demand_df, promo_df, summary, customer, product, agg_period='1d', forecast_periods=10) ```
MAE MSE RMSE R2 RMSLE MAPE
0 1.863 6.137 2.4773 -9.6113 0.9757 1.4333
[2.0, 1.0, 4.0, 0.0, 0.0, 3.0, 4.0, 2.0, 2.0, 4.0] ### AutoMLによる予測関数 predict_using_automl 引数： - df: 日付 (date)と需要 (deman)の列を含んだデータフレーム - future: 訓練用と予測用の未来のデータを入れたデータフレーム - best: 予測モデルのオブジェクト - horizon: 予測期間 返値： - fig: 予測結果を入れた図オブジェクト - predict: 予測結果を入れた配列 - error: 予測の2乗誤差 ------------------------------------------------------------------------ source ### predict_using_automl ``` python def predict_using_automl( df, future, best, horizon:NoneType=None ): ``` *AutoMLを用いた予測（回帰）* ### predict_using_automl関数の使用例 ``` python #fig, error, result = predict_using_automl(df, future, best[0]) #plotly.offline.plot(fig); ``` ![](06forecast_files/figure-commonmark/cell-116-output-1.png) ### 結果の可視化 - ‘residuals’ - Residuals Plot - ‘error’ - Prediction Error Plot - ‘cooks’ - Cooks Distance Plot - ‘rfe’ - Recursive Feat. Selection - ‘learning’ - Learning Curve - ‘vc’ - Validation Curve - ‘manifold’ - Manifold Learning - ‘feature’ - Feature Importance - ‘feature_all’ - Feature Importance (All) - ‘parameter’ - Model Hyperparameter - ‘tree’ - Decision Tree ``` python plot_model(best[0], "feature_all", save=True) ``` INFO:logs:Saving 'Feature Importance (All).png' in current active directory INFO:logs:Visual Rendered Successfully INFO:logs:plot_model() succesfully completed...................................... 'Feature Importance (All).png' ### 分類モデル ``` python from pycaret.classification import * #分類関連の関数のインポート ``` ``` python #def automl(df, horizon=None, num_bests = 5, sort = "R2"): """ Prediction using pycaret classification models """ horizon = "2020/10/01" num_bests = 1 sort = 'Accuracy' train_idx, valid_idx = prepare_index(df, horizon) ratio = len(train_idx)/(len(train_idx)+len(valid_idx)) df["Week"] = df.Week.astype(int) df["demand"] = df.demand.astype('category') #df["demand"] = df.demand.astype(int) reg = setup(df, target = 'demand', session_id=123, data_split_shuffle=False, fold_strategy= 'timeseries', fold=2, train_size=ratio, silent=True, verbose=False) best = compare_models(n_select=num_bests, sort = sort, cross_validation=False) result_df = pull() #return best, result_df ```
Model Accuracy AUC Recall Prec. F1 Kappa MCC TT (Sec)
ridge Ridge Classifier 0.5165 0.0000 0.2857 0.2979 0.3735 0.2591 0.3110 0.0100
rf Random Forest Classifier 0.5055 0.0000 0.2465 0.4723 0.4601 0.2719 0.2848 0.2100
et Extra Trees Classifier 0.4725 0.0000 0.2212 0.4560 0.4504 0.2373 0.2428 0.1700
lr Logistic Regression 0.4505 0.0000 0.1429 0.2030 0.2799 0.0000 0.0000 0.0200
knn K Neighbors Classifier 0.4505 0.0000 0.1429 0.2030 0.2799 0.0000 0.0000 0.0000
nb Naive Bayes 0.4505 0.0000 0.1429 0.2030 0.2799 0.0000 0.0000 0.0000
qda Quadratic Discriminant Analysis 0.4505 0.0000 0.1429 0.2030 0.2799 0.0000 0.0000 0.0100
lightgbm Light Gradient Boosting Machine 0.4396 0.0000 0.2687 0.4085 0.4170 0.2088 0.2120 1.1100
lda Linear Discriminant Analysis 0.4066 0.0000 0.3163 0.4123 0.3898 0.1879 0.1927 0.0100
xgboost Extreme Gradient Boosting 0.3956 0.0000 0.1998 0.4212 0.3911 0.1701 0.1751 1.3900
ada Ada Boost Classifier 0.3407 0.0000 0.3163 0.3963 0.3462 0.1451 0.1545 0.0700
gbc Gradient Boosting Classifier 0.3407 0.0000 0.2381 0.3764 0.3365 0.1332 0.1452 0.9300
catboost CatBoost Classifier 0.3077 0.0000 0.1233 0.3130 0.3093 0.0510 0.0518 3.6100
dt Decision Tree Classifier 0.2857 0.0000 0.2113 0.3375 0.2947 0.0424 0.0438 0.0100
svm SVM - Linear Kernel 0.2527 0.0000 0.1429 0.0639 0.1020 0.0000 0.0000 0.0200
INFO:logs:create_model_container: 0 INFO:logs:master_model_container: 0 INFO:logs:display_container: 2 INFO:logs:RidgeClassifier(alpha=1.0, class_weight=None, copy_X=True, fit_intercept=True, max_iter=None, normalize=False, random_state=123, solver='auto', tol=0.001) INFO:logs:compare_models() succesfully completed...................................... ``` python #検証データに対する予測 result = predict_model(best) #誤差の計算 dem = result.demand.astype(int) #error_ = ((result.Label - result.demand)**2).mean() error_ = ((result.Label - dem)**2).mean() #実需要量のプロット trace2 = go.Scatter( x = df.index, y = df.demand, mode = 'markers', name ="実需要量 Demand" ) #訓練＋未来のデータに対する予測 result = predict_model(best, future) #予測プロット trace1 = go.Scatter( x = result.index, y = result.Label, mode = 'markers + lines', name = "予測 Forecast" ) data = [trace1, trace2] layout = go.Layout( title ="AutoMLによる予測 Forecasting by AutoML", ) fig = go.Figure(data, layout) #return fig, error_, result plotly.offline.plot(fig); ``` INFO:logs:Initializing predict_model() INFO:logs:predict_model(estimator=RidgeClassifier(alpha=1.0, class_weight=None, copy_X=True, fit_intercept=True, max_iter=None, normalize=False, random_state=123, solver='auto', tol=0.001), probability_threshold=None, encoded_labels=False, round=4, verbose=True, ml_usecase=MLUsecase.CLASSIFICATION, display=None) INFO:logs:Checking exceptions INFO:logs:Preloading libraries INFO:logs:Preparing display monitor
Model Accuracy AUC Recall Prec. F1 Kappa MCC
0 Ridge Classifier 0.5165 0 0.2857 0.2979 0.3735 0.2591 0.3110
INFO:logs:Initializing predict_model() INFO:logs:predict_model(estimator=RidgeClassifier(alpha=1.0, class_weight=None, copy_X=True, fit_intercept=True, max_iter=None, normalize=False, random_state=123, solver='auto', tol=0.001), probability_threshold=None, encoded_labels=False, round=4, verbose=True, ml_usecase=MLUsecase.CLASSIFICATION, display=None) INFO:logs:Checking exceptions INFO:logs:Preloading libraries INFO:logs:Preparing display monitor ### 結果の可視化 ``` python plot_model(best, plot = 'feature', verbose=True) #fig_param = plot_model(best, plot = 'parameter') #interpret_model(best, plot = 'correlation') #plot_model(best, plot = 'error') ``` ![](06forecast_files/figure-commonmark/cell-121-output-1.png) ``` python plot_model(best, plot = 'learning') ``` ![](06forecast_files/figure-commonmark/cell-122-output-1.png) ``` python plot_model(best, plot = 'cooks') ``` ![](06forecast_files/figure-commonmark/cell-123-output-1.png) ## 需要の集約 上では、顧客ごとに需要予測を行なった。ここでは、倉庫（配送センター、流通センター）ごとの予測や、工場ごとの予測について考える。 そのためには、ロジスティック・ネットワーク設計モデルで製品の流れを最適化した結果を用いる必要がある。 最適化した結果はflow.csvに保存されているので、そこから生成されたデータフレーム flow_df と需要データフレーム demand_df を入れると、 地点（工場、倉庫）ごとに集約された需要量を返す関数を準備する。 ## 需要を地点ごとに集約する関数 aggregate_demand 引数： - demand_df : 需要データ - flow_df : 製品のフローを保管したデータ 返値： - aggregated_demand_df : 地点（工場、倉庫）と製品ごとに集約された需要データ - G: 製品をキーとし、フロー量を枝上に保管したグラフオブジェクト - Prod: 製品のリスト ------------------------------------------------------------------------ source ### aggregate_demand ``` python def aggregate_demand( demand_df, flow_df ): ``` *最適化された製品フローデータと需要データから、地点（工場、倉庫）ごとの集約された需要量を計算する関数* ### aggregate_demand関数の使用例 ``` python demand_df = pd.read_csv(folder+"demand.csv", index_col=0) flow_df = pd.read_csv(folder+"flow.csv", index_col=0) agg_demand_df, G, Prod = aggregate_demand(demand_df, flow_df) agg_demand_df.head() ```
date cust prod demand
0 2019-01-01 DC_さいたま市 A 1.0
1 2019-01-01 Plnt_Odawara A 5.0
2 2019-01-01 DC_さいたま市 B 1.0
3 2019-01-01 Plnt_Odawara B 4.0
4 2019-01-01 DC_さいたま市 C 3.0
``` python #グラフ G の描画 p ="B" pos = G[p].layout() nx.draw(G[p], pos=pos, node_color="Blue", node_size=10) ``` ![](06forecast_files/figure-commonmark/cell-127-output-1.png) ## 集約した需要をもとに安全在庫量を計算 単一ソース供給でない地点に対しては、安全在庫量を計算するための方法が確立されていない。 ここでは、異なるリード時間をもつ複数の地点から供給される場合には、ロジスティック・ネットワーク設計問題で求めた製品の年間フロー量をもとに、以下の方法で計算するものとする。 地点 *j* が、複数の地点 *i*(*i* ∈ *S**j*) からフロー量 *f**i**j* で供給を受けているものとする。供給地点をフロー量の大きさに応じてランダムに選ぶものと仮定し、その確率を *λ**i**j* とする。 $$ \lambda\_{ij} = \frac{f\_{ij}}{\sum\_{i \in S_j} f\_{ij}} $$ 地点 *i* から地点 *j* へのリード時間を *L*′*i**j* としたとき、地点 *j* へのリード時間（確率変数） *L**j* は、確率 *λ**i**j* で値 *L*′*i**j* をとる離散確率分布になる。 *L**j* の期待値 *E*\[*L**j*\] と分散 *V**a**r*(*L**j*) は、scipyモジュールを使えば簡単に計算できる。 これらの値を使えば、平均 *μ*、標準偏差 *σ* の正規分布にしたがう需要に対するリード時間内の需要量 *D* の期待値 *E*\[*D*\] と分散 *V**a**r*(*D*) は、以下のように計算できる。 *E*\[*D*\] = *μ**E*\[*L**j*\] *V**a**r*(*D*) = *σ**E*\[*L**j*\] + *μ*2*V**a**r*(*L**j*) 導出については、「はじめての確率論」（近代科学社）pp.92-94を参照されたい。 安全在庫量は、安全在庫係数を *z* としたとき、以下のように計算される。 $$ E\[D\] + z \sqrt{ Var(D) } $$ ## 集約した需要をもとに安全在庫量を計算する関数 compute_safety_stock 引数： - demand_df : 需要データ - agg_demand_df : 地点（工場、倉庫）と製品ごとに集約された需要データ - trans_df : 地点間のリード時間を保管したデータ - plnt_prod_df : 工場ごとの製品の生産リード時間を保管したデータ - Prod: 製品のリスト - G: 製品をキーとし、フロー量を枝上に保管したグラフオブジェクト 返値： - safety_df : 地点・製品ごとの安全在庫量を保管したデータ ------------------------------------------------------------------------ source ### compute_safety_stock ``` python def compute_safety_stock( demand_df, agg_demand_df, trans_df, plnt_prod_df, Prod, G ): ``` *倉庫・工場・顧客に対する安全在庫量の計算* ### compute_safety_stock関数の使用例 ``` python #trans_df = pd.read_csv(folder + "trans.csv", index_col=0) #plnt_prod_df = pd.read_csv(folder + "Plnt-Prod.csv", index_col=0) safety_df = compute_safety_stock(demand_df, agg_demand_df, trans_df, plnt_prod_df, Prod, G) # safety_df.to_csv(folder+"safety.csv") ```