予測手法と予測システム ABD-Forecast (AutoML Bayes Deep Forecast)
需要予測
サプライ・チェインにおける需要の予測のための入力は，一般的には以下の形式をもつデータである．
| 日付時刻 | 顧客 | 製品 | 需要量 | 売上 | 付加情報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関数を用いて生成したものである． 以下のように読み込んで使用する．
demand_with_promo_df = pd.read_csv(folder + "demand_with_promo.csv", index_col=0)
demand_with_promo_df.head()
パラメータ
- $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()
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
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: 顧客と製品のタプルをキーとし、モデルと予測オブジェクトのタプルを値とした辞書
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")
特徴量を入れた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: 顧客と製品のタプルをキーとし、モデルと予測オブジェクトのタプルを値とした辞書
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")
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")
len(forecast)
df_cv = cross_validation(model, initial='200 days', period='30 days', horizon = '60 days')
df_cv.head()
df_p = performance_metrics(df_cv, rolling_window=0.1)
df_p.head()
plot_cross_validation_metric(df_cv, metric='mse');
平均や標準偏差も計算
過去の需要系列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 期分の最大在庫量
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())
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())
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())
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()
future = make_future_dataframe(history_dates=df.index, periods=10, freq="1w", include_history=False)
future.head()
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から検証開始から予測開始の直前までを抜き出したデータフレーム
# 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)
予測に基づく動的発注量シートの生成関数 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
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")
顧客と製品を与えると、訓練用データフレームと予測用の未来のデータフレームを生成する関数
需要がない日を自動的に追加し，その需要を0と補完する．
引数：
- demand_df: 需要データフレーム
- customer: 予測を行いたい顧客
- product: 予測を行うたい製品
- promo_df: 追加特徴（プロモーション）情報を含んだデータフレーム（オプション）；予測期間に対しても追加情報を準備する必要がある． 予測期間内のデータが空 (NaN)の場合には，訓練データにも空 (NaN)のものがある必要がある．
- agg_period: 期間を表す文字列（例えば1日のときは"1d"）
- forecast_periods: 予測したい期間
返値：
- df: 訓練用データフレーム
- future: 訓練用と予測用の未来のデータを入れたデータフレーム
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_df: 需要データフレーム
- promo_df: 追加特徴（プロモーション）情報を含んだデータフレーム（オプション）；予測期間に対しても追加情報を準備する必要がある． 予測期間内のデータが空 (NaN)の場合には，訓練データにも空 (NaN)のものがある必要がある．
- agg_period: 期間を表す文字列（例えば1日のときは"1d"）
返値：
- demand_df: 訓練用データフレーム
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()
demand_df = pd.read_csv(folder+"demand_with_promo.csv")
find_horizon(demand_df, ratio=0.9)
#len(train_idx), len(valid_idx)
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)
lr_min= learn.lr_find()
print(lr_min)
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参照
train(learn, n_epoch = 30, lr_max =0.01, wd=2., cbs=ShowGraphCallback())
fig = plot_loss(learn)
plotly.offline.plot(fig);
fig, yhat, error, r2 = predict_using_dl(learn, df, future, classification=False)
print("error=", error, r2)
plotly.offline.plot(fig);
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);
Image("../figure/predict_using_all_dl.png")
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による特徴ベクトルの重要度の情報
#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)
feature = None
imp_df, fig = show_importance(importance, feature=feature, df=df)
plotly.offline.plot(fig);
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)
horizon="2020/10/01"
best, result_df = automl(df, horizon, 5, model_name=None)
result = predict_model(best[0], df)
result.head()
複数の製品・顧客に対して予測を行う関数（のサブ関数） 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になるように設定
複数の製品・顧客に対して予測を行う関数 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になるように設定
summary.r2.hist();
summary.to_csv("summary.csv")
#summary = pd.read_csv("summary.csv", index_col=0)
#summary.sort_values(by="r2",inplace=True)
summary.r2.hist();
すべての顧客・需要の予測結果を用いて，予測を行い図を返す関数 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: 未来の予測値（リスト）
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
すべての顧客・需要の予測結果を用いて，予測を行う関数 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: 未来の予測値（リスト）
# product = "D"
customer ="熊本市"
product ="H"
predict_using_automl_all(demand_df, promo_df, summary, customer, product, agg_period='1d', forecast_periods=10)
#plotly.offline.plot(fig);
結果の可視化
- '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
plot_model(best[0], "feature_all", save=True)
from pycaret.classification import * #分類関連の関数のインポート
"""
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
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);
plot_model(best, plot = 'feature', verbose=True)
#fig_param = plot_model(best, plot = 'parameter')
#interpret_model(best, plot = 'correlation')
#plot_model(best, plot = 'error')
plot_model(best, plot = 'learning')
plot_model(best, plot = 'cooks')
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()
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) } $$#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")