予測手法と予測システム ABD-Forecast (AutoML Bayes Deep Forecast)
while 1:
    target = np.array( [i%K    for i in range(2**depth) for j in range(n) ] )
    data = np.array([ [i +(n-j) for j in range(p)]  for i in range(n*2**depth)])
    reg = tree.DecisionTreeClassifier(max_depth = depth, splitter="random")
    reg.fit(data,target)
    if len(reg.tree_.threshold)==T:
        break
    
viz = dtreeviz(
    reg,
    data, 
    target,
    feature_names = ["X1","X2","X3","X4"],
    target_name = "y",
    class_names = [0, 1, 2]
) 
display(viz)
G cluster_legend node1 2020-12-11T12:23:52.125581 image/svg+xml Matplotlib v3.3.3, https://matplotlib.org/ node4 2020-12-11T12:23:52.244195 image/svg+xml Matplotlib v3.3.3, https://matplotlib.org/ leaf2 2020-12-11T12:23:52.429823 image/svg+xml Matplotlib v3.3.3, https://matplotlib.org/ node1->leaf2 leaf3 2020-12-11T12:23:52.467132 image/svg+xml Matplotlib v3.3.3, https://matplotlib.org/ node1->leaf3 leaf5 2020-12-11T12:23:52.506239 image/svg+xml Matplotlib v3.3.3, https://matplotlib.org/ node4->leaf5 leaf6 2020-12-11T12:23:52.538942 image/svg+xml Matplotlib v3.3.3, https://matplotlib.org/ node4->leaf6 node0 2020-12-11T12:23:52.358364 image/svg+xml Matplotlib v3.3.3, https://matplotlib.org/ node0->node1 < node0->node4 legend 2020-12-11T12:23:51.968199 image/svg+xml Matplotlib v3.3.3, https://matplotlib.org/

需要予測

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

日付時刻 顧客 製品 需要量 売上 付加情報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{eqnarray} &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{eqnarray} $$

乗法的季節変動:

$$ \begin{eqnarray} &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{eqnarray} $$

手法

  • 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[source]

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

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

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

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

需要がない日を自動的に追加し,その需要を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"
print(c,p)
agg_period ="1d"
df, future = make_forecast_df(demand_df, customer=c, product=p,promo_df= promo_df, agg_period="1w", forecast_periods = 10)
df.head()
仙台市 A
demand promo_0 Week promo_1 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 82 0 1 0 2019 1 6 6 6 False False False False False False 1546732800
1 91 0 2 0 2019 1 13 6 13 False False False False False False 1547337600
2 84 0 3 1 2019 1 20 6 20 False False False False False False 1547942400
3 88 0 4 0 2019 1 27 6 27 False False False False False False 1548547200
4 96 0 5 0 2019 2 3 6 34 False False False False False False 1549152000

全てのデータを用いて深層学習による予測をするためのデータフレームの生成関数

データ数が不足している場合には、顧客や製品もカテゴリー変数として扱うことによって十分なデータを確保することができる。

引数:

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

train_idx, valid_idx =prepare_index(df, horizon="2020/11/21")
len(train_idx), len(valid_idx)
(691, 39)

訓練

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

学習器を作成する関数 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);