Edited at

SIGNATE お弁当の売り上げをLightGBMで予測してみた



とある会社のカフェフロアで販売されているお弁当の販売数を予測するSIGNATEのコンペ「お弁当の需要予測」で作成したモデルを解説したいと思います。

こちら練習問題になっていますが、元ネタは2015年に実際に行われていたこちらのコンペのようです。

このモデルのRMSEは6.69380で、本投稿時6位相当でした。

データ分析の基本的な流れは網羅できていると思いますので、参考になれば幸いです。


目次

以下の順に沿って解説していきます。


  1. データの可視化

  2. トレンド抽出

  3. 特徴量エンジニアリング

  4. 特徴選択

  5. モデル作成

  6. テストデータの前処理と予測


1. データの可視化

まず目的変数である販売数の変動を可視化します。

train = pd.read_csv('../train.csv')

test = pd.read_csv('../test.csv')

train['y'].plot

ダウンロード.png


販売数をみて気になったこと:



  • チュートリアルにもある通り、初めの数十レコードは傾向が違う。ので、丸っと落としました

  • 明らかに右下がりのトレンドがあります。方針としては、トレンド抽出 → トレンド抜きのデータでモデル作成&予測 → トレンドを戻す

  • スパイクが多数存在します。スパイクを考慮せずにデータの平均を取ったりすると、おかしなことになりそう...


トレーニングデータとテストデータに大きな差異がないかも確認します

train.describe(include='all')

スクリーンショット 2019-03-27 23.22.00.png

test.describe(include='all')

スクリーンショット 2019-03-27 23.29.14.png


2. トレンド抽出

販売数のトレンドを抽出します。

statsmodelsのdecomposeを使う手もありますが、今回は最小二乗法によるカーブフィッティングでトレンドを抽出し、それを販売数から引き算することでトレンドを除いた目的変数を作成しました。

また、フィッティングに用いる関数は二次関数を使うより、ある値に漸近する指数関数を使う方が好ましいと考えました。

(販売数の減少量はどこかで飽和すると考えられるので)

from scipy.optimize import curve_fit

def func(x, a, b, c):
return a * np.exp(-b * x) + c

xs = train.index.values
ys = train['y'].values

popt, pcov = curve_fit(func, xs, ys)

a1, b1, c1 = popt
y_reg = a1 * np.exp(-b1 * xs) + c1

plt.figure()
plt.plot(train['y'])
plt.plot(y_reg)

ダウンロード (1).png

初めの数十レコードを落としたデータに当てはめました。良い感じです。

トレンドを除去します。

new_y = train['y'] - y_reg

new_y = pd.DataFrame({'new_y': new_y})

new_y.plot()

ダウンロード (2).png


3. 特徴量エンジニアリング

説明変数を一つ一つ見ていきます。


datetime

splitして月と日にちの部分だけ使います。

ただ、テストデータにしか存在しない月があるので、学習に月を使うのは微妙と思い今回は落としました。使うにしても加工が必要です。

train['day'] = train['datetime'].apply(lambda x : int(x.split("-")[2]))

「day」と「販売数」の関係性をplotlyで可視化します。

(インタラクティブなグラフなので本来なら触ると統計量を見れるのですが、今回は画像でご勘弁ください...)

import plotly.offline as offline

offline.init_notebook_mode(connected=False)

import plotly.graph_objs as go

layout = go.Layout(
autosize=False,
width=950,
height=500,
margin=go.layout.Margin(
l=80,
r=50,
b=50,
t=10,
pad=4),
xaxis=dict(
title='day',
titlefont=dict(
size=14,
),
showticklabels=True,
tickangle=0,
tickfont=dict(
size=12,
),
),
yaxis=dict(
title='トレンドを除いた販売数',
titlefont=dict(
size=14,
),
showticklabels=True,
tickfont=dict(
size=10,
),

))

data = [go.Box( x=train['day'], y=train['new_y'] )]

fig = go.Figure(data=data, layout=layout)
offline.iplot(fig, filename='example', show_link=False,
config={"displaylogo":False, "modeBarButtonsToRemove":["sendDataToCloud"]})

スクリーンショット 2019-03-27 23.52.08.png

plotlyで描画し、17日のBOXを触ってキャプチャしたのが上図です。

グラフを見ると何らかの周期性があり、paydayとの関係性も存在するはずなので残しておきます。


week

販売数の季節性に大きな影響を与えていそう。

「One hot / Label encoding」か「販売数の Mean/Median encoding」で数値に変換します。

(次元を増やしたくない&スパイクの影響を受けずに周期性を反映させたいので Median encodingしました。)

week_encoded = train.groupby('week').new_y.median()

week
week_encoded


7.854857


-8.968654


-1.908993


-3.708110


-8.816003

train['week'] = train['week'].replace(week_encoded)


soldout

「売り切れた日」と「売れ残った日」では販売数の意味合いが違うので、モデルをそれぞれ分けることも考えましたが今回はいっしょくたに作成します。


name

まずは 「name」と「販売数」の関係性を可視化します。

その前に、なるべく純粋な人気度を見たいので、先ほど作成した「week」を引き算し、曜日がもたらす季節性を超簡易的に除去した販売数を用いて可視化します。

train['new_new_y'] = train['new_y'] - train['week']

スクリーンショット 2019-03-28 10.58.28.png


グラフから読み取れること:


  • 人気メニューと不人気メニューがある

  • カレーが人気メニューのようにみえる

それぞれについて見ていきます。


まずはカレーから

チュートリアルの通り、「カレー」というワードを含むnameを1、含まないnameを0にラベルづけします。

train['curry'] = train['name'].apply(lambda x : 1 if x.find("カレー") >=0 else 0)

スクリーンショット 2019-03-28 11.09.27.png


次に人気メニューと不人気メニューをラベル付けします

今回はベースライン(トレンドとweekを差し引いた販売数)から15以上大きいnameを人気メニュー、-15より小さいnameを不人気メニューとしました。

popular_menu = set(train[train['new_new_y']>15].name)

{'さわら焼味噌掛け',

'さんま辛味焼',

'ひやしたぬきうどん・炊き込みご飯',

'キーマカレー',

'サバ焼味噌掛け',

'チキンカレー',

'チンジャオロース',

'チーズハンバーグ',

'ハンバーグカレーソース',

'ポークカレー',

'ポークハヤシ',

'メンチカツ',

'ロコモコ丼',

'手作りロースカツ',

'海老クリーミ―クノーデル',

'牛丼',

'牛肉筍煮',

'豚味噌メンチカツ',

'酢豚',

'酢豚orカレー',

'鶏のピリ辛焼き',

'鶏の唐揚げおろしソース',

'鶏の親子煮',

'鶏チリソース'}

train['popular'] = train['name'].apply(lambda x : 1 if x in popular_menu else 0)

不人気メニューも同様に処理します。

unpopular_menu = set(train[train['new_new_y']<-15].name)

train['unpopular'] = train['name'].apply(lambda x : 1 if x in unpopular_menu else 0)


nameをアノテーションしてみる

次に、トレーニングデータの「name」とテストデータの「name」で積集合を取ってみると、実はテストデータにしか含まれていない「name」が多数あることがわかります。

そこで説明力を増すために、以下の変数で「name」を0/1でアノテーションします。

食材:   beef、pork、chicken、fish、other(その他魚介)、vegi(野菜)

テイスト: japanese(和風)、western(洋風)、chinese(中華)

調理法:  grilled、sauteed、stewed、fried、steamed

ラベル付けすると、以下のようになります。

スクリーンショット 2019-03-28 0.11.27.png


  • JanomeやMecabで分かち書き→Word2vecなどでベクトル化する方法もたくさんあると思いますが、マニュアルでアノテーションした方が今回はお弁当の特徴をうまく表現できると思いました...ただしデータ数が少ない場合にしかできませんが...


  • また、この弁当名のベクトル化は前述の「popular」「unpopular」よりも汎用的なエンコーディングになっているので、モデルの汎化性能を高めているはずです。



kcal

欠損値が存在するので、トレーニングデータ全体の平均値で補間します。

train['kcal'] = train['kcal'].fillna(train['kcal'].mean())


remarks と event

情報量は少ないですが、重要な情報がないか精査します。

まずは「remarks」と「販売数」の関係性を可視化します。

スクリーンショット 2019-03-28 1.26.43.png

お楽しみメニューが大人気なので、特徴量に加えます。

単純に「お楽しみメニュー」に合致するremarksを1、それ以外を0にラベルづけします。

train['fun'] = train['remarks'].apply(lambda x: 1 if x=="お楽しみメニュー" else 0)

eventも同様に可視化しましたが、特に情報がなかったので今回は落とします。


payday

給料日の前後で販売数の傾向が異なると予想されます。今回は0/1ラベルでそのまま使用しました。

train['payday'] = train['payday'].fillna(0)


weather

「One hot / Label encoding」か「販売数の Mean/Median encoding」で数値に変換します。

次元が増えると嫌なので、「快晴」「晴れ」「曇」「雨およびその他」の4種で Median encoding しました。

weather_word = ['快晴','晴','曇','雨','雪','雷']

row_index = []
number = 1

for x in weather_word:

row_index = train['weather'].str.contains(x, na=False)
train.loc[row_index, ['weather']] = number

if number < 4:
number += 1
else:
number = 4

row_index = []

weather_encoded = train.groupby('weather').new_y.median()
train['weather'] = train['weather'].replace(weather_encoded)


precipitation

言わずもがなweatherと相関があるので、今回は落としました。

train = train.drop(columns = ['precipitation'])


temperature

そのままでも使えそうですが、前日との差分や月平均との差分を取ったほうが説明力が高いと考えました。

そこで今回は月平均との差分を「temp」として説明変数に加えました。

(「暖かい日が続いたが、ある日突然冷え込むとその日は外に出たくなくなる」みたいなことがあるので、販売数との相関がありそうと考えた。)

train['month'] = train['datetime'].apply(lambda x : int(x.split("-")[1]))

temp_mean = train.groupby('month').temperature.mean()
train['month'] = train['month'].replace(temp_mean)
train['temp'] = train['temperature'] - train['month']


最後に、pairplotを作成します

外れ値の有無や、目的変数との相関を可視化して確認します。

今回は極端な外れ値がなかったため、そのまま使用します。

ダウンロード (3).png


4. 特徴選択

どの特徴量をモデルに採用するか取捨選択します。

やり方は色々あります。


  • 説明変数と目的変数の相関をみて選択する(特徴量同士の相関が強い場合、どちらか落とすことも)

  • RandomForest系の学習器を用いて重要度を算出し、選択する

  • SBSアルゴリズムを用いて選択する

  • 直感的に重要だと感じる特徴量を入れてみて試行錯誤する

今回はRandomForestの重要度を用いて特徴量を選定しました。

こちらの投稿が大変参考になります。

(XGBoostやLightGBMを用いて特徴量の重要度を算出することもできますが、それらの重要度に従って特徴量を採用しても今回はあまり精度が上がりませんでした。この辺りは試行錯誤かなと思います...)

base = ['soldout','kcal','day','payday','temperature','temp']

weather = ['weather']
week = ['week']
name = ['fun','curry','popular','unpopular']
annotation = ['beef','pork','chicken','fish','other','vegi','japanese','chinese',
'western','grilled','sauteed','stewed','fried','steamed']

feature_x = base + weather + week + name + annotation
feature_y = ['new_y']

data_x = train[feature_x]
data_y = train[feature_y]

from sklearn.ensemble import RandomForestRegressor

rf_reg = RandomForestRegressor(n_estimators=10000)
rf_reg = rf_reg.fit(data_x, data_y)

fti = rf_reg.feature_importances_

dic_arr = {'importance':fti, 'feature':feature_x}
pd.DataFrame(dic_arr).sort_values('importance', ascending=False).reset_index(drop=True).loc[:15, :]

feature
importance

0
popular
0.247882

1
fun
0.178017

2
week
0.114100

3
unpopular
0.113638

4
temperature
0.086833

5
day
0.057744

6
temp
0.047977

7
curry
0.032029

8
kcal
0.030388

9
beef
0.015847

10
western
0.011213

11
weather
0.009339

12
japanese
0.007453

13
soldout
0.006305

14
payday
0.006071

15
pork
0.005581

今回はシンプルに、RandomForestの重要度が上位16の特徴量をそのまま入力しました。

(今回のRandomForestによる重要度算出では n_estimatorsを10000以上に設定すると、だいたいいつも同じ特徴量が顔を出すようになりました。)


5. モデル作成

学習器ですが、今回はオーソドックスにXGBoostとLightGBMでモデルを作成しました。

結果的にLightGBMの方が精度が良かったため、そちらを採用しました。

(XGBoostは今回のトレーニングデータに容易にoverfitしてしまい、それを抑えるのが難しかったです...)

コーディングはこちらのサイトを参考にしています。


※ 本来ならば、時系列データの交差検証で ShuffleSplit してしまうのは宜しくないと思います。しかし、今回はトレンドを除去し、季節性を説明する「week」や「day」も特徴に加えたので、シャッフルしてもなんとかなると考えました...



パラメータ調整

勾配ブースティング木系学習器チューニングのベストプラクティスは有名なものがあるので、そちらを参考に調整しました。

手順は以下の通りです。

1. 大きめの学習率(0.1くらい)で最適な木の数(n_estimators)を決定する

2. 木に関するパラメータ(colsample_bytree, max_depth, subsample など)をチューニングする

3. 正則化パラメータ(reg_alpha, reg_lambda)をチューニングする

4. より小さい学習率で最適な木の数を決定する

詳しくはこちらのリンクをご覧ください!

LGBMR = lgb.LGBMRegressor()

cv_split = model_selection.ShuffleSplit(n_splits = 10, test_size = .3, train_size = .7,
random_state = 0 )

base_results = model_selection.cross_validate(LGBMR, data_x[feature_x], data_y[feature_y],
scoring = 'mean_squared_error', cv = cv_split)
LGBMR.fit(data_x[feature_x], data_y[feature_y])

print("Training score mean: {:.2f}". format
(np.sqrt(abs(tune_model.cv_results_['mean_train_score'][tune_model.best_index_]))))
print("Test score mean: {:.2f}". format
(np.sqrt(abs(tune_model.cv_results_['mean_test_score'][tune_model.best_index_]))))

パラメータ調整前のRMSE:

Training score mean: 7.86

Test score mean: 12.74


param_grid = {'boosting_type': ['gbdt'],
'class_weight': [None],
'colsample_bytree': [0.9],
'importance_type': ['split'],
'learning_rate': [0.05],
'max_depth': [2],
'min_child_samples': [2],
'min_child_weight': [0],
'min_split_gain': [0],
'n_estimators': [200],
'n_jobs': [-1],
'num_leaves': [3],
'objective': [None],
'random_state': [None],
'reg_alpha': [0.01],
'reg_lambda': [0.01],
'silent': [True],
'subsample': [0.1],
'subsample_freq': [0]}

tune_model = model_selection.GridSearchCV(lgb.LGBMRegressor(), param_grid=param_grid,
scoring = 'mean_squared_error', cv = cv_split)
tune_model.fit(data_x[feature_x], data_y[feature_y])

print("Training score mean: {:.5f}". format
(np.sqrt(abs(tune_model.cv_results_['mean_train_score'][tune_model.best_index_]))))
print("Test score mean: {:.5f}". format
(np.sqrt(abs(tune_model.cv_results_['mean_test_score'][tune_model.best_index_]))))

パラメータ調整後のRMSE:

Training score mean: 6.75292

Test score mean: 10.72868


試しに、作成したモデルをトレーニングデータに当てはめてみます

train['new_y'] = tune_model.predict(data_x[feature_x])

original = train['y']
pred_y = train['new_y']
final_train = y_reg + pred_y

plt.figure()
plt.plot(original, label='original', linestyle='-')
plt.plot(final_train, label='trend + pred_y', linestyle='-')
plt.legend(loc='best')

ダウンロード (4).png


6. テストデータの前処理と予測

テストデータをトレーニングデータと同様に前処理します。


「week」のようにカテゴリカルデータを数値にエンコーディングした場合は、以下のようにトレーニングデータと同じエンコード対応表を用いて処理します。


test['week'] = test['week'].replace(week_encoded)

test['weather'] = test['weather'].replace(weather_encoded)
test['popular'] = test['name'].apply(lambda x : 1 if x in popular_menu else 0)
test['unpopular'] = test['name'].apply(lambda x : 1 if x in unpopular_menu else 0)

最後に、作成したモデルをテストデータに適用します。

data_x = test

test['y'] = tune_model.predict(data_x[feature_x])

トレンドの関数も外挿する感じでテストデータに当てはめ、テストデータのトレンド成分を予測、最後にそれを足し算します。

test['index_new'] = test.index + train.index.max() + 1

xs = test['index_new'].values
y_reg = a1 * np.exp(-b1 * xs) + c1

test['y'] = test['y'] + y_reg

あとは提出用データを作成するだけです!


まとめ

怪しげなアノテーション以外はベーシックな手法だけで予測できていると思います。

参考になれば幸いです。

※ 怪しい解説に対するご指摘もお待ちしております。


参考文献

お弁当の需要予測チュートリアル

ランダムフォレスト系ツールで特徴量の重要度を測る

A Data Science Framework: To Achieve 99% Accuracy

Complete Guide to Parameter Tuning in XGBoost (with codes in Python)

[Python] Plotlyでぐりぐり動かせるグラフを作る