【Pythonで需要予測①】重回帰分析で時系列の出荷数予測をしてみた

2024年5月18日

需要予測手法としてはARIMAモデルや、それに季節性を考慮したSARIMAモデルなどの時系列分析が知られていますが、普通の重回帰分析でもできるはずです。

今回はある物流センターからの6ヶ月間に渡る日々の出荷データを使って、重回帰分析でどこまでの精度で需要予測できるかを試してみました。

 

当ページのリンクには広告やプロモーションが含まれています。

出荷データを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の列に1~7のデータが入りましたが、これらは曜日を表しているのであって、数字の大きさに意味はありません。

例えば1は月曜日を、5は金曜日を表していますが、金曜日の方が月曜日より5倍凄いわけではありません。

これは1~52の数字を持つ週番号でも同じことです。

 

重回帰分析においてはそのような誤解を避けるために、ダミー変数に変換します。

曜日であれば1から7のダミー変数の列を作って、例えば月曜日であれば1の列に1を、金曜日であれば5の列に1を入れます。

週番号であれば1から52のダミー変数を作ります。

このように質的データ(カテゴリカルデータ)をダミー変数を使って変換することを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に、出荷数を目的変数yに入れます。

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()にxとyを入れて実行させます。

これにより重回帰式の係数が決まります。

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

 

予測誤差を評価する

予測値と実績値を見比べてみると週変動もトレンドも捉えていて、重回帰式のモデルが良くフィットしているように見えます。

これを定量的に調べてみます。

 

二乗平均平方根誤差

予測値と実測値の差である残差のデータは182個あります。

この残差をそのまま足し合わせてもプラスとマイナスで打ち消し合ってしまうので、残差の二乗を足し合わせます。

しかしそうすると合計値はかなり大きな値になってしまいますので、データ個数である182で割って一旦平均を出し、更にこれの平方根を取ります。

こうすることで、予測値や実績値と同じくらいの指標に変換することができるというのが二乗平均平方根誤差(Root Mean Squared Error:RMSEの考え方です。

これを計算してみると次のようになりました。

np.sqrt(mean_squared_error(data_df2.S00239, data_df2.Prediction))
…
239.10932937924403

239と計算されました。

実績値の平均を計算してみると917ですので、残差はその約1/4であることがわかります。

 

正規性の検定

次に残差がどのような分布になっているかを調べます。

【SCM分析3】Pythonで出荷データから周期変動と長期トレンドを分離する

の記事でもふれましたが、出荷数のような時系列データは周期トレンド残差に分解できることが多いので、残差は周期やトレンドとして拾い切れなかったゴミ、つまり小さな不規則変動データになっていることが望まれます。

確率分布でいえば、正規分布になっていることが望まれます。

これを確かめるためにKS検定を行ってみます。

【PythonのkstestでKS検定】出荷数が少なくても正規分布になることを検証する

 

まず、予測値-実績値で残差を計算します。

data_df2['residuals'] = data_df2.Prediction - data_df2.S00239
data_df2['residuals'].plot()

0を中心とした不規則変動になっているように見えます。

次にこれをヒストグラムにしてみます。

data_df2['residuals'].hist()

正規分布っぽく見えます。

よりわかりやすく表示するために、正規分布近似曲線と重ねて表示してみましょう。

Pythonで「ヒストグラム+正規分布の近似曲線」を作成する4通りの方法

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、p値が0.056になりました。

p値が5%より少し大きいので、正規分布であるという帰無仮説を棄却できません。

従って、この重回帰式は周期やトレンドをうまく捉えて予測しており、捉えきれなかった分が不規則変動として残差として残っているとみることができるでしょう。

 

SARIMAモデルによる需要予測

このように重回帰で需要予測することもできますが、SARIMAモデルを使う方が精度良い予測ができます。

またPythonのライブラリーを使うことができるので簡単です。

【Pythonで需要予測②】SARIMAモデルの自動構築で出荷数量を予測してみた