予測手法と予測システム BD-Forecast (Bayes Deep Forecast)
YouTubeVideo("D-L8I08CTxk")
YouTubeVideo("xYWGKiQ3Gzw")

指数平滑法

サプライ・チェインの古典手法に、指数平滑法と呼ばれる予測法がある。

航空機の乗客数の時系列データを用いて、様々な手法の比較を行う。

passengers = pd.read_csv("http://logopt.com/data/AirPassengers.csv")
#passengers["Month"] = pd.to_datetime(passengers.Month)
#passengers.set_index("Month", inplace=True) #月をインデックスとする
#passengers.index.freq ="MS" #月のはじめを頻度にする
passengers.head()
Month #Passengers
0 1949-01 112
1 1949-02 118
2 1949-03 132
3 1949-04 129
4 1949-05 121
df = passengers

freq = "MS" #月のはじめ
date_col = "Month"
target_col = "#Passengers"
feature_col = [ ]
horizon = 24 #最後の24ヶ月で検証する

try:
    df.reset_index(inplace=True)
except ValueError:
    pass
df = df.loc[:,[date_col, target_col]]
df.loc[:, date_col] = pd.to_datetime(df[date_col])
df.set_index(date_col, inplace=True) #月をインデックスとする
df.index.freq ="MS"

train = df[:-horizon]
#test = df[-horizon:]

fit1 = SimpleExpSmoothing(train).fit(smoothing_level=0.2,optimized=False)

start, end = 0, len(df)-1
predict1 = fit1.predict(start, end)

#predict1
passengers = pd.read_csv("http://logopt.com/data/AirPassengers.csv")
#passengers["Month"] = pd.to_datetime(passengers.Month)
#passengers.set_index("Month", inplace=True) #月をインデックスとする
#passengers.index.freq ="MS" #月のはじめを頻度にする
#passengers.head()

df = passengers

freq = "MS" #月のはじめ
date_col = "Month"
target_col = "#Passengers"
feature_col = [ ]
horizon = 24 #最後の24ヶ月で検証する

try:
    df.reset_index(inplace=True)
except ValueError:
    pass
df = df.loc[:,[date_col, target_col]]
df.loc[:, date_col] = pd.to_datetime(df[date_col])
df.set_index(date_col, inplace=True) #月をインデックスとする
df.index.freq ="MS"

train = df[:-horizon]

fit1 = SimpleExpSmoothing(train[target_col]).fit(smoothing_level=0.2,optimized=True)
fit2 = Holt(train[target_col]).fit(smoothing_level=0.8, smoothing_slope=0.2, optimized=True)
fit3 = ExponentialSmoothing(train[target_col], seasonal_periods=12, trend='add', seasonal='add').fit( optimized=True)
#start, end = len(df)-horizon, len(df)
start, end = 0, len(df)-1
predict1 = fit1.predict(start, end)
predict2 = fit2.predict(start, end)
predict3 = fit3.predict(start, end)

model_sarimax = SARIMAX(train[target_col],                                                                
            order=(1, 1, 1),
            seasonal_order=(1, 1, 1, 12),
            enforce_stationarity=False,
            enforce_invertibility=False)
fit4 = model_sarimax.fit()
predict4 = fit4.predict(start,end)
train.reset_index(inplace=True)
train.rename(columns={date_col: "ds", target_col: "y"}, inplace=True)
model_prophet = Prophet(daily_seasonality=False, weekly_seasonality=False, yearly_seasonality=True, seasonality_mode='multiplicative')
model_prophet.fit(train)
future = model_prophet.make_future_dataframe(periods=horizon, freq=freq)
forecast = model_prophet.predict(future)
df.loc[:,"Exp.Smooth."] = predict1.values
df.loc[:,"Holt"] = predict2.values
df.loc[:,"Holt-Winter"] = predict3.values
df.loc[:,"SARIMAX"] = predict4.values
df.loc[:,"Prophet"] = forecast.loc[:,"yhat"].values

#stacked_df = pd.DataFrame(df[-horizon:].stack(), columns=["num"]).reset_index()
stacked_df = pd.DataFrame(df[:].stack(), columns=["num"]).reset_index()
fig = px.line(stacked_df, x=date_col, y="num", color="level_1")
plotly.offline.plot(fig);
methods = ["Exp.Smooth.", "Holt", "Holt-Winter", "SARIMAX", "Prophet"]

metrics_list =[ [] for m in methods]
for i, m in enumerate(methods):
    y_pred = df[m][-horizon:] 
    y_true = df[target_col][-horizon:] 

    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 1217.831294 5719.201365 917.449416
mean_absolute 82.410342 270.694737 30.068210 69.895789 25.245934
mean_squared_log 0.055542 1.431378 0.005690 0.031673 0.003957
median_abusolute 51.102678 254.302112 30.110738 62.913468 24.692161
max 247.102678 517.826696 72.157078 129.229563 73.292586

Holt-Winterの季節変動と傾向変動を考慮した指数平滑法

$$ \begin{array}{ l l } & s_0 = x_0\\ & s_{t} = \alpha (x_{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 (x_{t}-s_{t-1}-b_{t-1})+(1-\gamma)c_{t-L}\\ & y_{t+m} = s_t + m b_t+c_{t-L+1+(m-1)\mod L} \end{array} $$

パラメータ

  • $L$ :周期 (seasonal_periods)
  • $\gamma$: 平滑化定数 (smoothing_level)
  • $\beta$: 傾向変動に対する平滑化定数 (smoothing_slope)
  • $\alpha$: 季節変動に対する平滑化定数 (smoothing_seasonal)

手法

  • SimpleExpSmoothing: 単純指数平滑法( $\gamma$ のみ用いる)
  • Holt: Holt法 ($\beta$ と $\gamma$ を用いる)
  • ExponentialSmoothing: Holt-Winterの指数平滑法(すべて用いる)

最後の2年間(24ヶ月)をテストデータとする。Holt-Winter法に対しては、パラメータは自動調整する。

passengers = pd.read_csv("http://logopt.com/data/AirPassengers.csv")
passengers["Month"] = pd.to_datetime(passengers.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
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.617864 91.776078 116.333985 115.617864 115.617864 -0.097102 -0.097102 -0.097102 -0.097102 -0.097102 -0.097102 0.0 0.0 0.0 104.391084
1 1949-02-01 117.308760 89.054286 112.738917 117.308760 117.308760 -0.141567 -0.141567 -0.141567 -0.141567 -0.141567 -0.141567 0.0 0.0 0.0 100.701744
2 1949-03-01 118.836021 107.984640 130.919967 118.836021 118.836021 0.004975 0.004975 0.004975 0.004975 0.004975 0.004975 0.0 0.0 0.0 119.427234
3 1949-04-01 120.526917 105.888404 128.954356 120.526917 120.526917 -0.031221 -0.031221 -0.031221 -0.031221 -0.031221 -0.031221 0.0 0.0 0.0 116.763905
4 1949-05-01 122.163268 107.848608 130.574211 122.163268 122.163268 -0.027509 -0.027509 -0.027509 -0.027509 -0.027509 -0.027509 0.0 0.0 0.0 118.802622
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
139 1960-08-01 450.295360 538.777453 565.676565 445.470906 455.116559 0.228100 0.228100 0.228100 0.228100 0.228100 0.228100 0.0 0.0 0.0 553.007938
140 1960-09-01 452.690181 463.239223 489.439250 447.449417 457.889348 0.051897 0.051897 0.051897 0.051897 0.051897 0.051897 0.0 0.0 0.0 476.183668
141 1960-10-01 455.007750 402.354480 428.188872 449.398529 460.641473 -0.088045 -0.088045 -0.088045 -0.088045 -0.088045 -0.088045 0.0 0.0 0.0 414.946765
142 1960-11-01 457.402572 352.902235 377.334503 451.362188 463.476228 -0.201518 -0.201518 -0.201518 -0.201518 -0.201518 -0.201518 0.0 0.0 0.0 365.227696
143 1960-12-01 459.720141 392.610280 419.363238 453.120206 466.237987 -0.116138 -0.116138 -0.116138 -0.116138 -0.116138 -0.116138 0.0 0.0 0.0 406.329050

144 rows × 16 columns

passengers["Exp.Smooth."] = predict1
passengers["Holt"] = predict2
passengers["Holt-Winter"] = predict3
passengers["SARIMAX"] = predict4
passengers["Prophet"] = forecast.loc[:,"yhat"].values

stacked_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);
#Image("../figure/sarimax.png")
passengers
#Passengers Exp.Smooth. Holt Holt-Winter SARIMAX Prophet
Month
1949-01-01 112 NaN NaN NaN NaN 104.391084
1949-02-01 118 NaN NaN NaN NaN 100.701744
1949-03-01 132 NaN NaN NaN NaN 119.427234
1949-04-01 129 NaN NaN NaN NaN 116.763905
1949-05-01 121 NaN NaN NaN NaN 118.802622
... ... ... ... ... ... ...
1960-08-01 606 374.897322 92.268388 554.606378 506.802984 553.007938
1960-09-01 508 374.897322 80.363471 455.607848 405.637677 476.183668
1960-10-01 461 374.897322 68.458554 409.052710 360.736984 414.946765
1960-11-01 390 374.897322 56.553638 361.330193 311.758931 365.227696
1960-12-01 432 374.897322 44.648721 400.129819 338.665680 406.329050

144 rows × 6 columns

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                   

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
  • kwargs: その他のProphetの引数

Prophetの引数と既定値:

  • growth='linear' :傾向変動の関数.規定値は線形.ロジスティック曲線にするに'logistic'に設定する.
  • changepoints=None : 傾向変更点のリスト
  • changepoint_range =0.8 : 傾向変化点の候補の幅(先頭から何割を候補とするか)
  • n_changepoints=25 : 傾向変更点の数
  • yearly_seasonality='auto' : 年次の季節変動を考慮するか否か
  • weekly_seasonality='auto' : 週次の季節変動を考慮するか否か
  • daily_seasonality='auto' : 日時の季節変動を考慮するか否か
  • holidays=None : 休日のリスト
  • seasonality_prior_scale=10.0 : 季節変動の事前分布のスケール値(パラメータの柔軟性をあら表す)
  • holidays_prior_scale=10.0 : 休日のの事前分布のスケール値(パラメータの柔軟性をあら表す)
  • changepoint_prior_scale=0.05 : 傾向変更点の事前分布のスケール値(パラメータの柔軟性をあら表す)
  • mcmc_samples=0 : MCMC法のサンプル数
  • interval_width=0.80 : 不確実性の幅
  • uncertainty_samples=1000 : 不確実性の幅を計算する際のサンプル数

返値:

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

  • dic: 顧客と製品のタプルをキーとし、モデルと予測オブジェクトのタプルを値とした辞書

forecast_demand[source]

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

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

forecast_demand関数の使用例

if FASTAI:
    demand_df = pd.read_csv(folder+"demand.csv")

    cust_selected = [ demand_df["cust"].iloc[350] ]
    prod_selected = [ demand_df["prod"].iloc[350] ]
    print(cust_selected, prod_selected)
    dic = forecast_demand(demand_df, cust_selected, prod_selected, agg_period ="1d", forecast_periods =10, cumulative=False,
                         daily_seasonality=False, weekly_seasonality=False, yearly_seasonality=False)
    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")
['徳島市'] ['A']
徳島市 A

特徴量を入れた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
  • kwargs: その他のProphetのキーワード引数

Prophetの引数と既定値:

  • growth='linear' :傾向変動の関数.規定値は線形.ロジスティック曲線にするに'logistic'に設定する.
  • changepoints=None : 傾向変更点のリスト
  • changepoint_range =0.8 : 傾向変化点の候補の幅(先頭から何割を候補とするか)
  • n_changepoints=25 : 傾向変更点の数
  • yearly_seasonality='auto' : 年次の季節変動を考慮するか否か
  • weekly_seasonality='auto' : 週次の季節変動を考慮するか否か
  • daily_seasonality='auto' : 日時の季節変動を考慮するか否か
  • holidays=None : 休日のリスト
  • seasonality_prior_scale=10.0 : 季節変動の事前分布のスケール値(パラメータの柔軟性をあら表す)
  • holidays_prior_scale=10.0 : 休日のの事前分布のスケール値(パラメータの柔軟性をあら表す)
  • changepoint_prior_scale=0.05 : 傾向変更点の事前分布のスケール値(パラメータの柔軟性をあら表す)
  • mcmc_samples=0 : MCMC法のサンプル数
  • interval_width=0.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, **kwargs)

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

forecast_demand_with_features関数の使用例

if FASTAI:
    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] ]
    prod_selected = [ demand_df["prod"].iloc[1] ]
    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,
                         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']

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

if FASTAI:
    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=True,
                         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
#fig = fp.plot_plotly(model, forecast)
#plotly.offline.plot(fig);
#fig.show()

交差検証

if FASTAI:
    df_cv = cross_validation(model, initial='220 days', period='100 days', horizon = '100 days')
    df_cv.head()
WARNING:fbprophet:Seasonality has period of 365.25 days which is larger than initial window. Consider increasing initial.
INFO:fbprophet:Making 1 forecasts with cutoffs between 2019-09-22 00:00:00 and 2019-09-22 00:00:00
ds yhat yhat_lower yhat_upper y cutoff
0 2019-09-23 1356.827158 1353.908975 1359.876367 1364.0 2019-09-22
1 2019-09-24 1356.157077 1353.016292 1359.255343 1368.0 2019-09-22
2 2019-09-25 1355.360231 1352.457403 1358.378218 1371.0 2019-09-22
3 2019-09-26 1353.303071 1350.493928 1356.395551 1377.0 2019-09-22
4 2019-09-27 1348.004275 1345.100149 1351.135613 1381.0 2019-09-22
if FASTAI:
    df_p = performance_metrics(df_cv, rolling_window=0.1)
    df_p.head()
horizon mse rmse mae mape coverage
0 10 days 2958.844066 54.395258 44.799240 0.032194 0.0
1 11 days 4279.223385 65.415773 55.595097 0.039868 0.0
2 12 days 6061.847382 77.857867 67.814721 0.048522 0.0
3 13 days 8284.903307 91.021444 81.242462 0.058006 0.0
4 14 days 11097.661465 105.345439 95.810634 0.068273 0.0
if FASTAI:
    plot_cross_validation_metric(df_cv, metric='rmse');

平均や標準偏差も計算

過去の需要系列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を用いた需要予測(顧客と製品に対して)

# ProphetのWarningがたくさん出る場合には、以下を生かす。
import warnings
warnings.simplefilter('ignore')
if FASTAI:
    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)

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

if FASTAI:
    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())

forecast_prophet_prod[source]

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

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

if FASTAI:
    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())

深層学習ライブラリ fastai を用いた需要予測

make_future_dataframe[source]

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

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

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

引数:

  • demand_df: 需要データフレーム
  • customer: 予測を行いたい顧客
  • product: 予測を行うたい製品
  • promo_df: 追加特徴(プロモーション)情報を含んだデータフレーム(オプション)
  • 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関数の使用例

if FASTAI:
    promo_df = pd.read_csv(folder+"promo.csv", index_col=0)
    demand_df = pd.read_csv(folder+"demand_with_promo.csv")
    #demand_df = pd.read_csv(folder+"liondemand.csv")
    c = demand_df["cust"].iloc[90] 
    p = demand_df["prod"].iloc[1] 
    print(c,p)
    agg_period ="1d"
    forecast_periods = 10
    df, future = make_forecast_df(demand_df, customer=c, product=p,promo_df= promo_df, agg_period="1d", forecast_periods =1)
    df.tail()
前橋市 B
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
725 1 0 0 2020 12 52 26 5 361 False False False False False False 1608940800
726 0 0 0 2020 12 52 27 6 362 False False False False False False 1609027200
727 0 0 0 2020 12 53 28 0 363 False False False False False False 1609113600
728 1 0 0 2020 12 53 29 1 364 False False False False False False 1609200000
729 0 1 1 2020 12 53 30 2 365 False False False False False False 1609286400
if FASTAI:
    demand_df = pd.read_csv(folder+"aggregate_demand.csv")

    c = demand_df["cust"].iloc[1] 
    p = demand_df["prod"].iloc[1] 

    print(c,p)
    agg_period ="1d"
    forecast_periods = 10
    df, future = make_forecast_df(demand_df, customer=c, product=p,promo_df=None, agg_period="1d", forecast_periods =1)
    df.tail()
Plnt_Odawara A
demand 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
360 1.0 2019 12 52 27 4 361 False False False False False False 1577404800
361 4.0 2019 12 52 28 5 362 False False False False False False 1577491200
362 4.0 2019 12 52 29 6 363 False False False False False False 1577577600
363 5.0 2019 12 1 30 0 364 False False False False False False 1577664000
364 5.0 2019 12 1 31 1 365 True False True False True False 1577750400

訓練

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

学習器を作成する関数 make_learn

引数:

  • df: 日付 (date)と需要 (deman)の列を含んだデータフレーム
  • horizon: 予測期間
  • classification: 分類問題としてモデルを作成する場合にはTrue、回帰問題としてモデルを作成する場合には Falseを入れる。

返値:

  • learn: fastaiの学習器

make_learn[source]

make_learn(df, horizon=0, classification=True)

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

make_learn関数の使用例

if FASTAI:
    learn = make_learn(df, horizon=100, classification=False)

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

lr_find[source]

lr_find(learn)

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

if FASTAI:
    lr_find(learn)

訓練を行う関数 train

引数:

train[source]

train(learn, **kwargs)

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

if FASTAI:
    train(learn, cyc_len=10, max_lr=0.01)
epoch train_loss valid_loss r2_score time
0 9.242335 8.273754 -0.730178 00:00
1 8.155198 8.342188 -0.748181 00:00
2 6.721432 7.617674 -0.587057 00:00
3 5.626612 5.892003 -0.185719 00:00
4 4.653332 5.509470 -0.098495 00:00
5 3.901257 5.292382 -0.025324 00:00
6 3.280552 5.328099 -0.031117 00:00
7 2.809872 5.278369 -0.027509 00:00
8 2.450302 5.287752 -0.029664 00:00
9 2.165404 5.292372 -0.032019 00:00

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

引数:

  • df: 日付 (date)と需要 (deman)の列を含んだデータフレーム
  • future: 訓練用と予測用の未来のデータを入れたデータフレーム
  • classification: 分類のとき True、回帰のとき False

返値:

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

predict_using_dl[source]

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

深層学習を用いた予測

predict_using_dl関数の使用例

if FASTAI:
    fig, yhat, error = predict_using_dl(df, future, learn, classification=True)
    print("error=", error)
    #plotly.offline.plot(fig);
    #fig.show()
Image("../figure/predict_using_dl.png")
error= 21.17

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

引数:

  • df: 日付 (date)と需要 (deman)の列を含んだデータフレーム
  • horizon: 予測期間

返値:

  • forest: ランダム森のモデルオブジェクト

random_forest[source]

random_forest(df, horizon)

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

random_forest関数の呼び出し例

if FASTAI:
    horizon = 100
    forest = random_forest(df, horizon)

ランダム森による予測

引数:

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

返値:

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

predict_using_rf[source]

predict_using_rf(df, future, forest, horizon)

ランダム木を用いた予測

if FASTAI:
    fig, yhat_rf, error_rf =  predict_using_rf(df, future, forest, horizon)
    #fig.show()
Image("../figure/predict_using_rf.png")

需要の集約

上では、顧客ごとに需要予測を行なった。ここでは、倉庫(配送センター、流通センター)ごとの予測や、工場ごとの予測について考える。

そのためには、ロジスティック・ネットワーク設計モデルで製品の流れを最適化した結果を用いる必要がある。

最適化した結果はflow.csvに保存されているので、そこから生成されたデータフレーム flow_df と需要データフレーム demand_df を入れると、 地点(工場、倉庫)ごとに集約された需要量を返す関数を準備する。

需要を地点ごとに集約する関数 aggregate_demand

引数:

  • demand_df : 需要データ
  • flow_df : 製品のフローを保管したデータ

返値:

  • agg_demand_df : 地点(工場、倉庫)と製品ごとに集約された需要データ
  • G: 製品をキーとし、フロー量を枝上に保管したグラフオブジェクト
  • Prod: 製品のリスト

aggregate_demand[source]

aggregate_demand(demand_df, flow_df)

最適化された製品フローデータと需要データから、地点(工場、倉庫)ごとの集約された需要量を計算する関数

aggregate_demand関数の使用例

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 0.0
1 2019-01-01 Plnt_Odawara A 1.0
2 2019-01-01 DC_さいたま市 B 13.0
3 2019-01-01 Plnt_Odawara B 48.0
4 2019-01-01 DC_さいたま市 C 7.0
agg_demand_df.to_csv(folder+"aggregate_demand.csv")
#グラフ G の描画
p ="B"
pos = G[p].layout()
nx.draw(G[p], pos=pos,  node_color="Blue", node_size=10)

集約した需要をもとに安全在庫量を計算

単一ソース供給でない地点に対しては、安全在庫量を計算するための方法が確立されていない。 ここでは、異なるリード時間をもつ複数の地点から供給される場合には、ロジスティック・ネットワーク設計問題で求めた製品の年間フロー量をもとに、以下の方法で計算するものとする。

地点 $j$ が、複数の地点 $i (i \in S_j)$ からフロー量 $f_{ij}$ で供給を受けているものとする。供給地点をフロー量の大きさに応じてランダムに選ぶものと仮定し、その確率を $\lambda_{ij}$ とする。

$$ \lambda_{ij} = \frac{f_{ij}}{\sum_{i \in S_j} f_{ij}} $$

地点 $i$ から地点 $j$ へのリード時間を $L'_{ij}$ としたとき、地点 $j$ へのリード時間(確率変数) $L_j$ は、確率 $\lambda_{ij}$ で値 $L'_{ij}$ をとる離散確率分布になる。

$L_j$ の期待値 $E[L_j]$ と分散 $Var(L_j)$ は、scipyモジュールを使えば簡単に計算できる。 これらの値を使えば、平均 $\mu$、標準偏差 $\sigma$ の正規分布にしたがう需要に対する、リード時間内の需要量 $D$ の期待値 $E[D]$ と分散 $Var(D)$ は、以下のように計算できる。

$$ E[D] = \mu E[L_j] $$$$ Var(D) = \sigma E[L_j] + \mu^2 Var(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 : 地点・製品ごとの安全在庫量を保管したデータ

compute_safety_stock[source]

compute_safety_stock(demand_df, agg_demand_df, trans_df, plnt_prod_df, Prod, G)

倉庫・工場・顧客に対する安全在庫量の計算

compute_safety_stock関数の使用例

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