予測手法と予測システム ABD-Forecast (AutoML Bayes Deep Forecast)

需要予測

サプライ・チェインにおける需要の予測のための入力は,一般的には以下の形式をもつデータである.

日付時刻 顧客 製品 需要量 売上 付加情報1 付加情報2 ...
2020年1月3日(何時何分) A商店 製品α 10 2600 特売日 ...
2020年1月3日 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関数を用いて生成したものである. 以下のように読み込んで使用する.

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$: 季節変動係数
  • $\hat{y}_{t+m}$: 期 $t$ における $m$期後の予測値
  • $L$ :季節変動の周期 (seasonal_periods)
  • $\alpha$: 平滑化定数 (smoothing_level)
  • $\beta$: 傾向変動に対する平滑化定数 (smoothing_slope)
  • $\gamma$: 季節変動に対する平滑化定数 (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: 単純指数平滑法( $\gamma$ のみ用いる)
  • Holt: Holt法 ($\beta$ と $\gamma$ を用いる)
  • ExponentialSmoothing: Holt-Wintersの指数平滑法(すべて用いる)

最後の2年間(24ヶ月)をテストデータとする。Holt-Winters法に対しては、パラメータは自動調整とし, Box-Cox変換をしてデータを正規化しておく(引数のuse_boxcoxをTrueにする).

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 + \phi_{1}y_{t-1} + \phi_{2}y_{t-2} + \cdots + \phi_{p}y_{t-p} + \varepsilon_{t} + \theta_{1}\varepsilon_{t-1} + \cdots + \theta_{q}\varepsilon_{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$

SARIMAX 外生変数 (eXogenous variables) の追加

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)
df = passengers.reset_index()
df.rename(columns={"Month": "ds", "#Passengers": "y"}, inplace=True)
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
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);
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: 集約を行う期。例えば1週間を1期とするなら"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: 顧客と製品のタプルをキーとし、モデルと予測オブジェクトのタプルを値とした辞書

forecast_demand[source]

forecast_demand(demand_df, cust_selected, prod_selected, agg_period='1d', forecast_periods=1, cumulative=False, lb=0.0, ub=None, **kwargs)

Prophetを用いた需要予測(顧客と製品に対して)

forecast_demand関数の使用例

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: 集約を行う期。例えば1週間を1期とするなら"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: 顧客と製品のタプルをキーとし、モデルと予測オブジェクトのタプルを値とした辞書

forecast_demand_with_features[source]

forecast_demand_with_features(demand_df, cust_selected, prod_selected, promo_df=None, features=None, agg_period='1d', forecast_periods=1, cumulative=False, lb=0.0, ub=None, **kwargs)

Prophetを用いた需要予測(顧客と製品に対して):特徴ベクトルの追加

forecast_demand_with_features関数の使用例

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']

集約したデータに対する予測

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
len(forecast)
730

交差検証

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.
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
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
plot_cross_validation_metric(df_cv, metric='mse');

平均や標準偏差も計算

過去の需要系列dfをもとに、Prophetによる予測を行う関数 forecast_prophet を準備する。 比較として、(需要の定常性を仮定して)平均と標準偏差も算出する。

引数:

  • demand_df: 需要の系列を表すデータフレーム。顧客Cust、製品Prod,日時Date、需要量demandの列をもつものと仮定する。
  • agg_period: 集約を行う期。例えば1週間を1期とするなら"1w"とする。既定値は1日("1d")
  • forecast_period: 予測を行う未来の期数。既定値は1

返値:

以下の情報を列情報にもつデータフレーム

  • mean: 需要の平均
  • std : 需要の標準偏差
  • Prophet[0..forecast_period-1]: Prophetによる予測から計算された未来のforecast_periods 期分の最大在庫量

forecast_prophet[source]

forecast_prophet(df, agg_period='1d', forecast_periods=1)

Prophetを用いた需要予測

forecast_prophet_cust_prod[source]

forecast_prophet_cust_prod(demand_df, cust_selected, prod_selected, agg_period='1d', forecast_periods=1, cumulative=False)

Prophetを用いた需要予測(顧客と製品に対して)

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())

forecast_prophet_cust[source]

forecast_prophet_cust(demand_df, cust_selected, agg_period='1d', forecast_periods=1)

各顧客に対する(すべての製品の)合計需要の予測

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: []

forecast_prophet_prod[source]

forecast_prophet_prod(demand_df, prod_selected, agg_period='1d', forecast_periods=1)

各製品に対する(すべての顧客の)合計需要の予測

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"とする.

make_future_dataframe[source]

make_future_dataframe(history_dates, periods, freq='1d', include_history=True)

未来の時系列データフレームを生成する関数

make_future_dataframe関数の使用例

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
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日が1期となる.
  • num_features: (休日を除く)特徴の数
  • holiday: 日本のカレンダーに基づく休日を表す特徴を追加するときTrue

返値:

  • 需要予測用のExcelのWorkbook

make_excel_forecast[source]

make_excel_forecast(start, finish, period=1, period_unit='日', num_features=0, holiday=True)

make_excel_forecast関数の使用例

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から検証開始から予測開始の直前までを抜き出したデータフレーム

forecast_excel[source]

forecast_excel(wb, valid_start, predict_start)

forecast_excel関数の使用例

#valid_start = "2019-09-01"
#predict_start = "2020-12-30"

#wb = load_workbook("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)
lightgbm Light Gradient Boosting Machine 1.3320 3.1986 1.7885 0.9354 0.2967 0.2199 0.3500
rf Random Forest Regressor 1.3561 3.4884 1.8677 0.9296 0.3334 0.1814 0.1800
et Extra Trees Regressor 1.3085 3.5258 1.8777 0.9288 0.2845 0.1732 0.1600
xgboost Extreme Gradient Boosting 1.2509 3.5396 1.8814 0.9285 0.2547 0.2001 0.5300
gbr Gradient Boosting Regressor 1.4922 3.9734 1.9933 0.9198 0.3078 0.3087 0.0700
dt Decision Tree Regressor 1.5694 7.2361 2.6900 0.8539 0.3852 0.2130 0.0100
ridge Ridge Regression 2.8140 11.1862 3.3446 0.7742 0.6221 0.7606 0.0100
br Bayesian Ridge 2.8264 11.3965 3.3759 0.7699 0.6325 0.7731 0.0100
tr TheilSen Regressor 2.8347 11.4245 3.3800 0.7694 0.6128 0.7581 2.5500
lar Least Angle Regression 2.8122 11.5554 3.3993 0.7667 0.6453 0.7861 0.0200
ada AdaBoost Regressor 5.1129 32.5647 5.7065 0.3426 0.8812 1.0945 0.0700
lasso Lasso Regression 5.5731 38.8488 6.2329 0.2158 1.0129 1.3597 0.0100
omp Orthogonal Matching Pursuit 5.8594 49.4031 7.0287 0.0027 0.6792 0.8511 0.0100
en Elastic Net 5.9584 50.0735 7.0763 -0.0108 1.1412 1.6917 0.0100
par Passive Aggressive Regressor 5.9702 53.3333 7.3030 -0.0766 1.1413 1.6373 0.0000
lr Linear Regression 6.4090 61.5741 7.8469 -0.2430 1.2076 1.9134 0.0100
svm Support Vector Regression 6.6225 64.9033 8.0563 -0.3102 1.2211 1.9716 0.0200
llar Lasso Least Angle Regression 6.8024 67.6916 8.2275 -0.3665 1.2349 2.0349 0.0000
ard Automatic Relevance Determination 6.8024 67.6916 8.2275 -0.3665 1.2349 2.0349 0.0100
huber Huber Regressor 6.8059 67.7414 8.2305 -0.3675 1.2349 2.0347 0.0100
ransac Random Sample Consensus 7.0263 72.1502 8.4941 -0.4565 1.0488 0.8727 0.0900
kr Kernel Ridge 7.3367 76.3726 8.7391 -0.5417 1.2718 2.2102 0.1400
knn K Neighbors Regressor 9.5194 117.0961 10.8211 -1.3638 1.3903 2.8196 0.0000
mlp MLP Regressor 420.1910 176642.1063 420.2881 -3564.8733 4.2362 91.7979 0.2500

予測に基づく動的発注量シートの生成関数 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

dynamic_inv_excel[source]

dynamic_inv_excel(wb, valid_start, predict_start, valid_df, demand_df, h, fc, vc=0.0, z=2.33)

dynamic_inv_excel関数の使用例

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: 訓練用と予測用の未来のデータを入れたデータフレーム

make_forecast_df[source]

make_forecast_df(demand_df, customer, product, promo_df=None, agg_period='1d', forecast_periods=1)

顧客と製品を与えると、訓練用データフレームと予測用の未来のデータフレームを生成する関数

make_forecast_df関数の使用例

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: 訓練用データフレーム

make_forecast_df_all[source]

make_forecast_df_all(demand_df, promo_df=None, agg_period='1d', forecast_periods=1)

全データに対する訓練用データフレームを生成する関数

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

find_horizon[source]

find_horizon(demand_df, ratio=0.9)

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: 検証データのインデックス

prepare_index[source]

prepare_index(df, horizon=None)

#len(train_idx), len(valid_idx)

訓練

需要を連続量とする場合と離散量とする場合で異なる予測を行う。連続量の場合には回帰問題となり、離散量の場合には分類問題になる。

学習器を作成する関数 make_learn

引数:

  • df: 日付 (date)と需要 (deman)の列を含んだデータフレーム
  • horizon: 検証データを作成するときの閾値の日付;この日付より後のデータは検証データになる。
  • classification: 分類問題としてモデルを作成する場合にはTrue、回帰問題としてモデルを作成する場合には Falseを入れる。

返値:

  • learn: fastaiの学習器

make_learn[source]

make_learn(df, horizon=None, classification=False)

深層学習を用いた需要予測用のデータを生成

make_learn関数の使用例

#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)

適正な学習率を見つけるための関数 lr_find

lr_min, lr_steep = learn.lr_find() 
print(lr_min, lr_steep)
0.33113112449646 0.005248074419796467

lr_findの結果をPlotlyで描画する関数 plot_lr_find

plot_lr_find[source]

plot_lr_find(learn)

fig = plot_lr_find(learn)
plotly.offline.plot(fig);

訓練を行う関数 train

引数:

train[source]

train(learn, **kwargs)

fit_one_cycle法を用いて訓練を行う関数

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

Plotlyで損出関数を描画する関数 plot_loss

引数:

  • learn: 学習器
  • valid: テストデータも描画するときにTrue(既定値)

返り値:

  • fig: 損出関数のグラフ

plot_loss[source]

plot_loss(learn, valid=True)

fig = plot_loss(learn)
plotly.offline.plot(fig);

深層学習を用いて予測を行う関数 predict_using_dl

引数:

  • learn:学習器
  • df: 日付 (date)と需要 (deman)の列を含んだデータフレーム
  • future: 訓練用と予測用の未来のデータを入れたデータフレーム
  • classification: 分類問題のとき True、それ以外のとき False

返値:

  • fig: 予測結果を入れた図オブジェクト
  • predict: 予測結果を入れた配列
  • error: 予測の2乗誤差
  • r2: 決定係数

predict_using_dl[source]

predict_using_dl(learn, df, future, classification=False)

深層学習を用いた予測

predict_using_dl関数の使用例

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

全てのデータで訓練した学習器による予測関数 predict_using_all_dl

引数:

  • learn: 学習器
  • df: 該当する顧客と製品の需要データを入れたデータフレーム
  • furure: 予測したい全期間の需要以外のデータを入れたデータフレーム
  • customer: 顧客名
  • product: 製品名
  • classification: 分類問題のとき True、それ以外のとき False

返値:

  • fig: 予測結果を入れた図オブジェクト
  • predict: 予測結果を入れた配列
  • error: 予測の2乗誤差

predict_using_all_dl[source]

predict_using_all_dl(learn, df, future, customer, product, classification=True)

深層学習を用いた予測(全部で訓練)

predict_using_all_dl関数の使用例

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
Image("../figure/predict_using_all_dl.png")

埋め込み層の描画関数 draw_embedding

引数:

  • learn: 学習器(訓練済み)
  • cust: 顧客を描画するときTrue, 製品を描画するとき False
  • pca: 主成分分析で次元削減するときTrue、TSNEでするとき False

返値:

  • fig: Plotlyの図オブジェクト

draw_embedding[source]

draw_embedding(learn, cust=True, pca=True)

draw embedding

draw_embedding関数の使用例

fig = draw_embedding(learn, cust=False, pca=True)
plotly.offline.plot(fig);
fig = draw_embedding(learn, cust=True, pca=True)
plotly.offline.plot(fig);

scikit learnのランダム森によって予測をする関数 random_forest

引数:

  • df: 日付 (date)と需要 (deman)の列を含んだデータフレーム
  • horizon: 検証データを作成するときの閾値の日付;この日付より後のデータは検証データになる。
  • feature: 特徴ベクトルとして使用する列名のリスト; Noneのときには全ての列とする。
  • classification: 分類問題のときTrue、回帰問題のときFalse

返値:

  • forest: ランダム森のモデルオブジェクト
  • importance: permutation_importanceによる特徴ベクトルの重要度の情報

random_forest[source]

random_forest(df, horizon=None, feature=None, classification=False)

ランダム森による需要予測

random_forest関数の呼び出し例

#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: 重要度の棒グラフ

show_importance[source]

show_importance(importance, feature, df)

feature = None
imp_df, fig = show_importance(importance, feature=feature, df=df)
plotly.offline.plot(fig);

ランダム森による予測 predict_using_rf

引数:

  • df: 日付 (date)と需要 (deman)の列を含んだデータフレーム
  • future: 訓練用と予測用の未来のデータを入れたデータフレーム
  • forest: ランダム森のモデルオブジェクト
  • horizon: 予測期間

返値:

  • fig: 予測結果を入れた図オブジェクト
  • predict: 予測結果を入れた配列
  • error: 予測の2乗誤差
  • r2: 決定係数

predict_using_rf[source]

predict_using_rf(df, future, forest, horizon=None, feature=None)

ランダム木を用いた予測

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