トレタ データサイエンス研究所

Report
2018年3月29日

予約データから未来の客数を予測する – 実践編

この記事をシェアする

前回の記事では、飲食店における客数予測をどのようなステップで実施し、どの程度の予測精度が必要かについて議論しました。今回は、客数予測を自分の手でやってみたいという方に向けて、Pythonを使って予測モデルを構築する手順をご紹介しようと思います。

参考:予約データから未来の客数を予測する

予測モデルを構築する

今回は、Jupyter Notebook上で、疑似データ生成、モデル構築、可視化を行っています。

実行環境

  • Mac OS High Sierra 10.13.3
  • Python 3.6.4
  • jupyter 1.0.0
  • numpy 1.14.2
  • pandas 0.22.0
  • matplotlib 2.2.2
  • statsmodels 0.8.0

疑似データ生成

まずは、必要なライブラリをインポートします。

# 基本のライブラリ
import numpy as np
import pandas as pd

# グラフ描画
from matplotlib import pylab as plt
from matplotlib.pylab import rcParams
%matplotlib inline

# モデルのライブラリ
import statsmodels.api as sm

客数は全体の平均客数に対して、週(7日間)の周期的な変動と、季節(数十日単位)の変動が合わさって、一定のばらつきで正規分布に従って生成されると仮定し、疑似データを生成します。2017年の1年間を学習して、2018年1月の1ヶ月分を予測することを想定し、1年1ヶ月分のデータを生成します。
念のためですが、時系列データの各データは自己相関が発生することが想定されるため、正規分布で独立と仮定して疑似データを生成するのは、厳密には間違っています。ただ、今回は予測モデルを構築するための疑似データが欲しいだけですので、厳密な疑似データ生成の議論は置いておきます。

# 任意の期間を指定
start = pd.datetime(2017, 1, 1)  # 開始日
end = pd.datetime(2018, 1, 31)  # 終了日
date_index = pd.date_range(start=start, end=end)

# 疑似データを生成するパラメータ
mu = 40  # 客数の平均を40に設定
mu_week = [-2, -4, 0, 1, 7, 8, 3]  # 客数の平均40に対して、曜日ごとに平均が加減されるとする
s = 4  # 客数のばらつきを設定


# 季節項を任意の関数で表現
# 今回は365日の日付で線形の関数が選択されて、客数の平均に加減されるとする
def mu_season(x):
    if (x % 365) < 31:
        return -0.3*(x % 365) + 10
    elif 31 <= (x % 365) < 90:
        return 0.15*((x % 365)-31)
    elif 90 <= (x % 365) < 212:
        return -0.09*((x % 365)-90) + 10
    elif 212 <= (x % 365) < 243:
        return 0.3*((x % 365)-212)
    elif 243 <= (x % 365) < 273:
        return -0.3*((x % 365)-243)+10
    elif 273 <= (x % 365):
        return 0.11*((x % 365)-273)


# 疑似データを生成
Y = []
np.random.seed(1234)  # seedを設定して、実行ごとに同じ疑似データを生成するようにする
for i in range(len(date_index)):  # 指定した期間の日数だけ疑似データを生成する
    # 疑似データは正規分布に従うとし、平均は全体のmu、週変動のmu_week、季節変動のmu_seasonが加減され、分散はsで一定とする。
    Y.append(np.absolute(np.floor(np.random.normal(loc=mu + mu_week[i % 7] + mu_season(i), scale=s))))

# 疑似データに日付インデックスを付けて、pandasのDataFrameに格納
daily_customers = pd.DataFrame(data={'date': date_index, 'customers': Y}, index=date_index)

可視化

生成した疑似データが、設定した傾向になっているか確認します。まずは、1年間分のデータを可視化してみます。1週間の周期的な変動と、年間を通しての緩やかな季節変動が生成されていることが見て取れます。

# 生成したデータの1年間を抽出して可視化
ts_training = daily_customers['2017']
ax = ts_training['customers'].plot(figsize=(8, 4), grid=True, color='black')
ax.set_ylabel('customers')
fig = ax.get_figure()
fig.savefig('./ts_plot_1y.png')
客数推移(1年)

1年間の客数を日別に可視化

つぎに、1週間分のデータも可視化してみます。設定した通り、週の中で客数の少ない日が続く部分と、多い日が続く部分ができています。

# 生成したデータの1週間を抽出して可視化
ax = ts_training['2017-12-18':'2017-12-24']['customers'].plot(figsize=(8, 4), grid=True, color='black')
ax.set_ylabel('customers')
ax.set_ylim(30, 70)
fig = ax.get_figure()
fig.savefig('./ts_plot_1w.png')
客数推移(1週間)

1週間の客数を日別に可視化

予測モデル構築

今回の予測モデルは、「状態空間モデル」で構築します。状態空間モデルは時系列モデルの一種で、時間に沿って「状態」が変化し、その「状態」に誤差が乗って値が「観測」されるとするモデルです。状態空間モデルの解説は、それだけで1冊の本が書けてしまうほどですが、ライブラリを利用してモデル構築するだけであれば、わずかなコードで実行できます。
今回、「状態」は前日の状態と、1週間の周期的な状態の変化に影響されて決定するとしています。これ以外にも、「状態」は前日の状態だけに影響されるとか、「状態」は前日と1週間の周期と年間のトレンドに影響されるなど、異なるモデルを構築することができます。本来であれば、このような複数の状態空間モデルを作成し、モデル選択基準(例えば、赤池情報量規準(AIC))が一番良いモデルを採用するなどのモデル選定作業が必要ですが、今回は割愛して一つのモデルのみ構築します。

# 周期変動ありのローカルレベルモデル
# 周期的な変動がseasonなので、1週間の周期的な変動をseasonで表す
week_local_level_model = sm.tsa.UnobservedComponents(ts_training['customers'].to_frame(), level='local level', seasonal=7)

# モデルのパラメータをデータから学習させる
result_week_local_level_model = week_local_level_model.fit()

# 推定されたパラメタ一覧
print(result_week_local_level_model.summary())

# 推定された状態・トレンド・週の影響の描画
rcParams['figure.figsize'] = 15, 20
fig = result_week_local_level_model.plot_components()
fig.savefig('./ts_train.png')

学習されたパラメータを表示します。モデル選択基準のひとつである赤池情報量規準(AIC)は右上の方に表示されます。

                           Unobserved Components Results                            
====================================================================================
Dep. Variable:                    customers   No. Observations:                  365
Model:                          local level   Log Likelihood               -1029.951
                   + stochastic seasonal(7)   AIC                           2065.901
Date:                      Mon, 26 Mar 2018   BIC                           2077.601
Time:                              13:46:24   HQIC                          2070.551
Sample:                          01-01-2017                                         
                               - 12-31-2017                                         
Covariance Type:                        opg                                         
====================================================================================
                       coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------------
sigma2.irregular    14.3980      1.092     13.189      0.000      12.258      16.538
sigma2.level         0.4415      0.182      2.424      0.015       0.084       0.799
sigma2.seasonal      0.0066      0.016      0.417      0.677      -0.024       0.038
===================================================================================
Ljung-Box (Q):                       35.23   Jarque-Bera (JB):                 4.98
Prob(Q):                              0.68   Prob(JB):                         0.08
Heteroskedasticity (H):               0.87   Skew:                            -0.21
Prob(H) (two-sided):                  0.47   Kurtosis:                         3.39
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

構築されたモデルの結果を下図に描画しています。
上段では、黒線が実際のデータで、青線が予測モデルの値です。薄い青色は予測モデルの95%信頼区間です。薄い青色の信頼区間の中に黒線が入っていることから、構築されたモデルのパラメータは大きく外れてはいなさそうです。
中段は、抽出されたトレンドです。トレンドは今回のモデルに明示的に組み込んでいないため、トレンドの状態の変化は「前日の状態」の変化で説明していることになります。設定した季節変動がここに描画されていることが見て取れます。
下段は週の周期的な変動成分です。周期的な変動成分の名前はデフォルトでseasonとなっていますが、今回は周期的な変動は週で発生しているとしてモデルに組み込んでいますので、週の周期変動を捉えていることになります。

客数予測モデルの構築

予測の実施

1年間のデータで学習した予測モデルで、学習期間より先の未来を予測させてみます。2017/12/15から2017/12/31は学習データの期間ですが、データがある期間と無い期間での予測の変化を見るために含めておきます。

# 予測
pred = result_week_local_level_model.predict(pd.to_datetime('2017-12-15'), pd.to_datetime('2018-01-31'))

# 実データと予測結果の図示
ax = daily_customers['customers'].plot(figsize=(8, 4), grid=True, color='black')
ax.set_ylabel('customers')
ax = pred.plot(figsize=(8, 4), grid=True, color='red')
fig = ax.get_figure()
fig.savefig('./ts_pred.png')

黒線が実際のデータで、赤線が予測した客数です。2017/12/15から2017/12/31については、データがある期間ですので、予測と実績は大きく乖離していませんが、2018/1/1以降の学習時にデータが無い期間は、だんだんと予測と実績が乖離していっている様子が見て取れます。今回の予測モデルでは、「状態」は前日の状態と、週の周期的な変動だけで説明できるとしているので、前日の状態が無くなった期間では、前日の状態は変化せずに週の周期変動のみがそのまま続くとして予測してしまっています。
予測精度を上げるためには、季節変動をモデルの説明に組み込む必要がありそうです。この他にも、その日の平均気温が客数に影響しているということであれば、気温データを持ってきてモデルに組み込むなどの精度向上の工夫ができそうです。
ただ、このモデルでも2018/1/1から1週間であれば予測は当たっていそうです。予測モデルの目的であるアクションの仮説検証が1週間の予測で成り立つのであれば、このモデルでも十分な予測精度の可能性はあります。

予測の実施

この記事をシェアする
島田 哲朗
この記事を書いた人

コンピュータを含む環境を上手く設計することで、個人やチーム能力を自然に発揮できないかと考え、東京大学大学院 暦本研究室にてモバイルデバイスでのARを研究。2012年に卒業後、アクセンチュア株式会社のアナリティクス部署にて、通信、金融、アパレル、小売における数千万人規模の顧客データを分析・モデリングして、マーケティング、営業支援、商品企画、出店開発などへのデータ主導意志決定に従事。2017年12月より、トレタにて飲食店データを活用したビジネス構築に携わる。

新着記事