【Pythonで需要予測①】重回帰分析で時系列の出荷数予測をしてみた
需要予測手法としてはARIMAモデルや、それに季節性を考慮したSARIMAモデルなどの時系列分析が知られていますが、普通の重回帰分析でもできるはずです。
今回はある物流センターからの６ヶ月間に渡る日々の出荷データを使って、重回帰分析でどこまでの精度で需要予測できるかを試してみました。
出荷データをPythonに読み込んでグラフにする
560アイテムについて取得した6ヶ月間の日々の出荷データをデータフレームraw_dfに読み込みます。
raw_df = pd.read_csv('demand.csv')
raw_df['Date'] = pd.to_datetime(raw_df['Date']) # 日付データを日付型に変換
raw_df
今回も前回に続いて、この中のS00239のアイテムについて分析することにします。
data_df = raw_df[['Date', 'S00239']]
data_df = data_df.set_index('Date')
data_df
グラフにすると次のようになります。
plt.figure(figsize=(10, 5))
plt.plot(data_df)
plt.show()
曜日と週番号を付加する
このグラフを見ると、週変動と僅かなトレンドがあることがわかります。
従って、曜日と週番号を説明変数として出荷数を予測できるのではという仮説が立てられます。
そこで、日付データから曜日と週番号を取得して、新たな列として追加しておきましょう。
Pandasのdt.isocalendar()でどちらとも取得できます。
data_df['DOW'] = data_df.Date.dt.isocalendar().day
data_df['WK'] = data_df.Date.dt.isocalendar().week
data_df
カテゴリカルデータをダミー変数に変換する
これでDOWの列に１～７のデータが入りましたが、これらは曜日を表しているのであって、数字の大きさに意味はありません。
例えば１は月曜日を、５は金曜日を表していますが、金曜日の方が月曜日より５倍凄いわけではありません。
これは１～５２の数字を持つ週番号でも同じことです。
重回帰分析においてはそのような誤解を避けるために、ダミー変数に変換します。
曜日であれば１から７のダミー変数の列を作って、例えば月曜日であれば１の列に１を、金曜日であれば５の列に１を入れます。
週番号であれば１から５２のダミー変数を作ります。
このように質的データ（カテゴリカルデータ）をダミー変数を使って変換することをOne-Hotエンコーディングといいます。
このような変換はPandasのget_dummies()で一度に行うことができますが、はじめに変換前のデータタイプを’category’にしておくことに注意です。
data_df['DOW'] = data_df['DOW'].astype('category')
data_df['WK'] = data_df['WK'].astype('category')
data_df2 = pd.get_dummies(data_df)
data_df2
説明変数と目的変数に分ける
今回の重回帰分析は曜日と週番号から出荷数を予測しようとしているので、曜日と週間号を説明変数ｘに、出荷数を目的変数ｙに入れます。
x = data_df2.drop(['Date', 'S00239'], axis = 1).values
x
...
array([[ 1. , 0. , 0. , ..., 0. , 436. , -5.25],
[ 0. , 1. , 0. , ..., 0. , 732. , -29.25],
[ 0. , 0. , 1. , ..., 0. , 804. , -202. ],
...,
[ 0. , 0. , 0. , ..., 0. , 844. , -511.5 ],
[ 0. , 0. , 0. , ..., 0. , 516. , -331.75],
[ 0. , 0. , 0. , ..., 0. , 228. , -129. ]])y = data_df2.S00239.values
y
...array([ 430.75 , 702.75 , 602. , 817.5 ,
701. , 335.5 , 244.25 , 488.75 ,
748.5 , 764.25 , 1549. , 1153. ,
811. , 348.5 , 404. , 512.5 ,
… 817.25 , 877.25 , 840.25 , 1070.5 ,
955.5 , 560.25 , 519.5 , 1478.25 ,
1044.5 , 846.25 , 572. , 332.5 ,
184.25 , 99. ])
線形回帰を実行する
次にsklearn.linear_modelのLinearRegression()にｘとｙを入れて実行させます。
これにより重回帰式の係数が決まります。
model = LinearRegression()
model.fit(x, y)
予測値を表示する
できた重回帰式を使った予測値はpredictionで取得できます。
data_df2['Prediction'] = prediction
data_df2
この予測値を実績値と一緒に表示させておきます。
plt.plot(data_df2.Date, data_df2.S00239, label = 'Actual')
plt.plot(data_df2.Date, data_df2.Prediction, label = 'Prediction')
plt.legend()
予測誤差を評価する
予測値と実績値を見比べてみると週変動もトレンドも捉えていて、重回帰式のモデルが良くフィットしているように見えます。
これを定量的に調べてみます。
二乗平均平方根誤差
予測値と実測値の差である残差のデータは１８２個あります。
この残差をそのまま足し合わせてもプラスとマイナスで打ち消し合ってしまうので、残差の二乗を足し合わせます。
しかしそうすると合計値はかなり大きな値になってしまいますので、データ個数である１８２で割って一旦平均を出し、更にこれの平方根を取ります。
こうすることで、予測値や実績値と同じくらいの指標に変換することができるというのが二乗平均平方根誤差（Root Mean Squared Error：RMSE）の考え方です。
これを計算してみると次のようになりました。
np.sqrt(mean_squared_error(data_df2.S00239, data_df2.Prediction))
…
239.10932937924403２３９と計算されました。
実績値の平均を計算してみると９１７ですので、残差はその約１／４であることがわかります。
正規性の検定
次に残差がどのような分布になっているかを調べます。
【SCM分析3】Pythonで出荷データから周期変動と長期トレンドを分離する
の記事でもふれましたが、出荷数のような時系列データは周期／トレンド／残差に分解できることが多いので、残差は周期やトレンドとして拾い切れなかったゴミ、つまり小さな不規則変動データになっていることが望まれます。
確率分布でいえば、正規分布になっていることが望まれます。
これを確かめるためにKS検定を行ってみます。
【PythonのkstestでKS検定】出荷数が少なくても正規分布になることを検証する
まず、予測値－実績値で残差を計算します。
data_df2['residuals'] = data_df2.Prediction - data_df2.S00239
data_df2['residuals'].plot()
０を中心とした不規則変動になっているように見えます。
次にこれをヒストグラムにしてみます。
data_df2['residuals'].hist()
正規分布っぽく見えます。
よりわかりやすく表示するために、正規分布近似曲線と重ねて表示してみましょう。
Pythonで「ヒストグラム＋正規分布の近似曲線」を作成する４通りの方法
xmin, xmax = data_df2.residuals.min(), data_df2.residuals.max()
mean, std = data_df2.residuals.mean(), data_df2.residuals.std()
x_axis = np.linspace(xmin, xmax, 100)
pdf = norm.pdf(x_axis, mean, std)
pdf_graph = pd.DataFrame({"X":x_axis, "norm":pdf}).set_index("X")
pdf_graph.plot()
data_df2['residuals'].hist(density = True)
かなり正規分布に近い分布になっていることがわかります。
最後にKS検定してみましょう。
from scipy import stats
from scipy.stats import norm
para = norm.fit(data_df2.residuals)
stats.kstest(data_df2.residuals, 'norm', args = para)
…
KstestResult(statistic=0.098083181116975, pvalue=0.05627906261792881)KS値が0.098、ｐ値が0.056になりました。
ｐ値が５％より少し大きいので、正規分布であるという帰無仮説を棄却できません。
従って、この重回帰式は周期やトレンドをうまく捉えて予測しており、捉えきれなかった分が不規則変動として残差として残っているとみることができるでしょう。
SARIMAモデルによる需要予測
このように重回帰で需要予測することもできますが、SARIMAモデルを使う方が精度良い予測ができます。
またPythonのライブラリーを使うことができるので簡単です。