はじめに
以前記事にしていますが、Titanicコンペでは基本的な機械学習のアルゴリズムを試すことができました。
初心者向けと言われるコンペでしたが、最初は苦労しましたね…
まだまだ経験が足りない!ということで、このたび新しいコンペに取り組むことになりました。
また、今回は共同開発?の経験とするため、2人でチームを組んで参加しております。
参加するコンペはすでに終了している、Predict Future Sales です。
タスクは、店舗と製品別の日毎の売り上げデータをもとに、次月(2015年12月)の店舗と製品別の売り上げを予測するというものです。
Couseraの最終課題として使われたことで有名なようです。
今回は先人のNoteBookを参考にしながら、ランダムフォレストとLightGBMでモデルを作成してみました。
データの確認からやっていきましょう!
データの確認
まずは、必要なものをimportしていきます。
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
%matplotlib inline
import sklearn
import sklearn.metrics
import datetime as dt
import lightgbm as lgb
import os
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_validate, KFold
from sklearn.preprocessing import StandardScaler
from sklearn.preprocessing import LabelEncoder
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error
from lightgbm import plot_importance
importしたものをどのように使うかは後述していきます。
データを読み込んでいきます。
from zipfile import ZipFile
zip_file = ZipFile('competitive-data-science-predict-future-sales.zip')
zip_file.namelist() # Zipの中身を確認
# 出力
['item_categories.csv',
'items.csv',
'sales_train.csv',
'sample_submission.csv',
'shops.csv',
'test.csv']
容量節約のためにzipファイルのまま処理していきたいと思います。
train = pd.read_csv(zip_file.open('sales_train.csv'))
test = pd.read_csv(zip_file.open('test.csv'))
submission = pd.read_csv(zip_file.open('sample_submission.csv'))
category = pd.read_csv(zip_file.open('item_categories.csv'))
items = pd.read_csv(zip_file.open('items.csv'))
shops = pd.read_csv(zip_file.open('shops.csv'))
zip_file.openで開封。
それぞれのファイルについて、.head()で確認していきます。(長くなるので出力結果は省略)
さらに、処理を軽くするためにdowncastを定義しておきます。
def downcast(df, verbose=True):
start_mem = df.memory_usage().sum() / 1024**2
for col in df.columns:
dtype_name = df[col].dtype.name
if dtype_name == 'object':
pass
elif dtype_name == 'bool':
df[col] = df[col].astype('int8')
elif dtype_name.startswith('int') or (df[col].round() == df[col]).all():
df[col] = pd.to_numeric(df[col], downcast='integer')
else:
df[col] = pd.to_numeric(df[col], downcast='float')
end_mem = df.memory_usage().sum() / 1024**2
if verbose:
print('{:.1f}% compressed'.format(100 * (start_mem - end_mem) / start_mem))
return df
all_df = [train, shops, items, category, test]
for df in all_df:
df = downcast(df)
前処理
外れ値の確認・削除
trainデータ内における、item_cnt_day(製品の個数)と、item_price(製品の価格)を箱ひげ図で確認してみます。
fig,ax = plt.subplots(2,1,figsize=(10,4))
# 製品の個数
# 尺度の調整
plt.xlim(-300, 3000)
# 箱ひげ図を描画
ax[0].boxplot((sales_train.item_cnt_day) , labels=['train.item_cnt_day'], vert=False)
# 製品の価格
# 尺度の調整
plt.xlim(-1000, 350000)
# 箱ひげ図を描画
ax[1].boxplot((sales_train.item_price) , labels=['train.item_price'], vert=False)
plt.show()
箱ひげ図より、外れ値を確認できたので処理していきます。
# 外れ値の除去
train = train[train['item_price'] > 0]
train = train[train['item_price'] < 50000]
train = train[train['item_cnt_day'] <1000]
# 返品分(item_cnt_day<0)は除去
train = train[train['item_cnt_day'] > 0]
shopsデータを確認したところ、同じ内容のものがありました。
print(shops['shop_name'][0], '||', shops['shop_name'][57])
print(shops['shop_name'][1], '||', shops['shop_name'][58])
print(shops['shop_name'][10], '||', shops['shop_name'][11])
print(shops['shop_name'][39], '||', shops['shop_name'][40])
# 出力
!Якутск Орджоникидзе, 56 фран || Якутск Орджоникидзе, 56
!Якутск ТЦ "Центральный" фран || Якутск ТЦ "Центральный"
Жуковский ул. Чкалова 39м? || Жуковский ул. Чкалова 39м²
РостовНаДону ТРК "Мегацентр Горизонт" || РостовНаДону ТРК "Мегацентр Горизонт" Островной
これは統一してしまってよさそうです。
train.loc[train['shop_id'] == 0 ,'shop_id'] =57
train.loc[train['shop_id'] == 1 ,'shop_id'] =58
train.loc[train['shop_id'] == 10 ,'shop_id'] =11
train.loc[train['shop_id'] == 39 ,'shop_id'] =40
test.loc[test['shop_id'] == 0 ,'shop_id'] =57
test.loc[test['shop_id'] == 1 ,'shop_id'] =58
test.loc[test['shop_id'] == 10 ,'shop_id'] =11
test.loc[test['shop_id'] == 39 ,'shop_id'] =40
エンコード
shops
shop_name から地域名を .split()を用いて分割し抽出。 city に格納します。
shops['city'] = shops['shop_name'].apply(lambda x : x.split()[0])
shops['city'].unique()
# 出力
array(['!Якутск', 'Адыгея', 'Балашиха', 'Волжский', 'Вологда', 'Воронеж',
'Выездная', 'Жуковский', 'Интернет-магазин', 'Казань', 'Калуга',
'Коломна', 'Красноярск', 'Курск', 'Москва', 'Мытищи', 'Н.Новгород',
'Новосибирск', 'Омск', 'РостовНаДону', 'СПб', 'Самара', 'Сергиев',
'Сургут', 'Томск', 'Тюмень', 'Уфа', 'Химки', 'Цифровой', 'Чехов',
'Якутск', 'Ярославль'], dtype=object)
shops.loc[shops['city'] == '!Якутск', 'city'] ='Якутск'
shops内はカテゴリ変数なので、ラベルエンコーディングしていきます。
le = LabelEncoder()
shops['city'] = le.fit_transform(shops['city'])
不要なカラムを削除し、中身を確認します。
shops.drop('shop_name', axis=1, inplace=True)
shops.head()
# 出力
shop_id city
0 0 29
1 1 29
2 2 0
3 3 1
4 4 2
items
まずは不要なカラムを削除します。
items = items.drop(['item_name'], axis=1)
trainデータのdate_block_numより、はじめて製品を購入した月のカラムを作成します。
items['first_sale_date'] = train.groupby('item_id').agg({'date_block_num':'min'})['date_block_num']
items.head()
欠損値を確認します。
items[items['first_sale_date'].isna()]
# 欠損値は368 rows
はじめて製品を購入した月が欠損しているということは、2015年10月より後に購入したものであることが考えられるので、欠損値を34として処理していきます。
items['first_sale_date'] = items['first_sale_date'].fillna(34)
category
item_category_nameは空白とハイフンで名前が区切られており、
アクセサリー – PS2 のように表記されています。
まずは空白で区切り、categoryというカラムを作成します。
category['category'] = category['item_category_name'].apply(lambda x : x.split()[0])
カテゴリー数が5未満のものが確認できたので、その他として一括りにしておきます。
def make_etc(x):
if len(category[category['category']==x]) >=5:
return x
else:
return 'etc'
category['category'] = category['category'].apply(lambda x : make_etc(x))
中身を確認しておきます。
category.head()
# 出力
item_category_name item_category_id category
0 PC - Гарнитуры/Наушники 0 etc
1 Аксессуары - PS2 1 Аксессуары
2 Аксессуары - PS3 2 Аксессуары
3 Аксессуары - PS4 3 Аксессуары
4 Аксессуары - PSP 4 Аксессуары
カテゴリー変数をラベルエンコーディングしていきます。
lr = LabelEncoder()
category['category'] = lr.fit_transform(category['category'])
del category['item_category_name']
category.head()
# 出力
item_category_id category
0 0 0
1 1 1
2 2 1
3 3 1
4 4 1
trainデータを組み合わせる
df_trainを新たに作成します。
同じ月のshop_idとitem_idをグループ化し、連結します。
# 簡略化のため、itertoolsを導入
from itertools import product
df_train = []
# 同じ月のshop_id, item_idをグループ化
for i in train['date_block_num'].unique():
all_shop = train.loc[train['date_block_num'] == i, 'shop_id'].unique()
all_item = train.loc[train['date_block_num'] == i, 'item_id'].unique()
df_train.append(np.array(list(product([i], all_shop, all_item))))
idx_features = ['date_block_num', 'shop_id', 'item_id']
df_train = pd.DataFrame(np.vstack(df_train), columns=idx_features)
.head() で確認します。
df_train.head()
trainデータのdate_block_numを月データに、item_priceは平均に集計します。
group = train.groupby(idx_features).agg({'item_cnt_day' : 'sum', 'item_price' : 'mean'}).reset_index()
group = group.rename(columns={'item_cnt_day' : 'item_cnt_month', 'item_price' : 'item_price_mean'})
df_train = df_train.merge(group, on=idx_features, how='left')
# 不要な変数を掃除してメモリを確保
import gc
del group
gc.collect();
item_cnt_dayをカウントします。
group = train.groupby(idx_features).agg({'item_cnt_day' : 'count'}).reset_index()
group = group.rename(columns={'item_cnt_day' : 'item_count'})
df_train = df_train.merge(group, on=idx_features, how='left')
del group, train
gc.collect()
df_trainとtestを結合します。
test['date_block_num'] = 34
all_data = pd.concat([df_train, test.drop('ID', axis=1)], ignore_index=True, keys=idx_features)
all_data.fillna(0, inplace=True)
all_data = all_data.merge(shops, on='shop_id', how='left')
all_data = all_data.merge(items, on='item_id', how='left')
all_data = all_data.merge(category, on='item_category_id', how='left')
all_data = downcast(all_data)
del shops, items, category
gc.collect();
可視化
月ごとの売上数を可視化してみる。
group_month_sum = all_data.groupby('date_block_num').agg({'item_cnt_month' : 'sum'}).reset_index()
plt.figure(figsize=(11,5))
sns.barplot(x='date_block_num', y='item_cnt_month', data=group_month_sum)
plt.show()
続いてアイテムカテゴリごとの売上数(1万個以上売れたものを抽出)
group_cat_sum = all_data.groupby('item_category_id').item_cnt_month.sum().reset_index()
group_cat_sum = group_cat_sum[group_cat_sum['item_cnt_month'] > 10000]
plt.figure(figsize=(11,5))
sns.barplot(x='item_category_id', y='item_cnt_month', data=group_cat_sum)
plt.show()
続いて、shop_idごとの売上数(1万個以上のものを抽出)
group_shop_sum = all_data.groupby('shop_id').item_cnt_month.sum().reset_index()
group_shop_sum = group_shop_sum[group_shop_sum['item_cnt_month'] > 10000]
plt.figure(figsize=(11,5))
sns.barplot(x='shop_id', y='item_cnt_month', data=group_shop_sum)
plt.show()
特徴量の生成
# df ->大元のデータ
# mean_features -> 新たな特徴量カラム(1つ後のセルを参照)
# idx_features -> 上記を生成するためのカラム
def add_mean_features(df, mean_features, idx_features):
# エラーが起きていないかチェック
assert (idx_features[0] == 'date_block_num') and len(idx_features) in [2,3]
if len(idx_features) == 2:
feature_name = idx_features[1] + '_mean_sales'
else:
feature_name = idx_features[1] + '_' + idx_features[2] + '_mean_sales'
# idx_features毎の月売上数の平均をgroupに格納し、name変更
group = df.groupby(idx_features).item_cnt_month.mean().reset_index()
group = group.rename(columns={'item_cnt_month' : feature_name})
# 上記で算出したものをdf と結合する
df = df.merge(group, on=idx_features, how='left')
# mean_feature('item_idmean_sales', 'item_id_city_mean_sales')にattend
mean_features.append(feature_name)
# 不必要なデータを削除
del group
gc.collect()
return df, mean_features
# 新たな特徴量カラムを入れるための箱を作成
item_mean_features = []
# 'date_block_num', 'item_id' 毎の売上個数の平均をall_dataへ格納
all_data, item_mean_features = add_mean_features(
df = all_data,
mean_features = item_mean_features,
idx_features = ['date_block_num', 'item_id'])
# 'date_block_num', 'item_id','city' 毎の売上個数の平均をall_dataへ格納
all_data, item_mean_features = add_mean_features(
df = all_data,
mean_features = item_mean_features,
idx_features = ['date_block_num', 'item_id', 'city'])
all_dataにカラム名が2つ追加されていることを確認します。
item_mean_features
# 出力
['item_id_mean_sales', 'item_id_city_mean_sales']
先ほどと同様にshop_mean_featuresとして空の箱を用意していきます。
# 先ほど同様に空の箱を用意
shop_mean_features = []
# 'date_block_num', 'shop_id', 'item_category_id'毎の売上個数の平均をall_dataへ格納
all_data, shop_mean_features = add_mean_features(
df = all_data,
mean_features = shop_mean_features,
idx_features = ['date_block_num', 'shop_id', 'item_category_id'])
カラムが作成できていることを確認します。
shop_mean_features
# 出力
['shop_id_item_category_id_mean_sales']
ラグ特徴量の生成
引き続き、特徴量エンジニアリングとしてラグ特徴量を作成していきます。
ラグ特徴量は、特徴量が時系列データである時に有効に使えます。
# 3month の ラグ特徴量を作成する関数
def add_lag_features(df, lag_features_to_clip, idx_features,
lag_feature, nlags=3, clip=False):
df_tmp = df[idx_features + [lag_feature]].copy()
# ラグ特徴量を生成(nlagsでラグの期間を設定する)
for i in range(1, nlags+1):
# ラグ特徴量の名前を作成
lag_feature_name = lag_feature + '_lag' + str(i)
# df_tmpのカラム名を作成
df_tmp.columns = idx_features + [lag_feature_name]
df_tmp['date_block_num'] += i
df = df.merge(df_tmp.drop_duplicates(), on=idx_features, how='left')
df[lag_feature_name] = df[lag_feature_name].fillna(0)
# clip = Trueにした場合のif 式
# (Falseにした場合はlag_features_to_clipにカラムが格納されない)
if clip:
lag_features_to_clip.append(lag_feature_name)
df = downcast(df, False)
del df_tmp
gc.collect()
return df, lag_features_to_clip
ラグ特徴量を生成する関数を定義したら、実際に入れ込んでいきます。
lag_features_to_clip = []
idx_features = ['date_block_num', 'shop_id', 'item_id']
# 3ヶ月分のラグ特徴量を生成
all_data, lag_features_to_clip = add_lag_features(
df = all_data,
lag_features_to_clip = lag_features_to_clip,
idx_features = idx_features,
lag_feature = 'item_cnt_month',
nlags = 3,
clip = True)
確認してみると、
lag_features_to_clip
# 出力
['item_cnt_month_lag1', 'item_cnt_month_lag2', 'item_cnt_month_lag3']
作成できていることがわかります。
先ほどと同様に、item_count、item_price_meanのラグ特徴量を生成します。
# 'item_count'のラグ特徴量を生成
all_data, lag_features_to_clip = add_lag_features(
df = all_data,
lag_features_to_clip = lag_features_to_clip,
idx_features = idx_features,
lag_feature = 'item_count',
nlags = 3)
# 'item_price_mean' のラグ特徴量を生成
all_data, lag_features_to_clip = add_lag_features(
df = all_data,
lag_features_to_clip = lag_features_to_clip,
idx_features = idx_features,
lag_feature = 'item_price_mean',
nlags = 3)
all_data = all_data.drop(shop_mean_features, axis=1)
date_blocl_numが3より小さいものをdropします。
all_data = all_data.drop(all_data[all_data['date_block_num'] < 3].index)
その他の特徴量
3ヶ月のラグ特徴量の平均値を生成し、’item_cnt_month_lag_mean’へ格納します。
all_data['item_cnt_month_lag_mean'] = all_data[['item_cnt_month_lag1',
'item_cnt_month_lag2',
'item_cnt_month_lag3']].mean(axis=1)
Evaluationに記載がありますが、データは(0, 20)でclipしておく必要があります。
all_data[lag_features_to_clip + ['item_cnt_month', 'item_cnt_month_lag_mean']]\
= all_data[lag_features_to_clip + ['item_cnt_month', 'item_cnt_month_lag_mean']].clip(0, 20)
下のコードでは、ゼロ除算の処理を行なっています。
pandasでゼロ除算があった場合は、infもしくは-infとなるのでそのままでは計算できません。
そのため、inf or -inf を nan に置き換えてから fillna で 0 として処理しています。
all_data['lag_grad1'] = all_data['item_cnt_month_lag1'] / all_data['item_cnt_month_lag2']
all_data['lag_grad1'] = all_data['lag_grad1'].replace([np.inf, -np.inf], np.nan).fillna(0)
all_data['lag_grad2'] = all_data['item_cnt_month_lag2'] / all_data['item_cnt_month_lag3']
all_data['lag_grad2'] = all_data['lag_grad2'].replace([np.inf, -np.inf], np.nan).fillna(0)
新規顧客のカラムを生成します。
all_data['brand_new'] = all_data['first_sale_date'] == all_data['date_block_num']
はじめて購入してからの期間を算出します。
all_data['duration_after_first_sale'] = all_data['date_block_num'] - all_data['first_sale_date']
all_data = all_data.drop('first_sale_date', axis=1)
%は余りの値を返すので、12で割って月を算出します。
all_data['month'] = all_data['date_block_num']%12
最後のデータ形成です。
all_data = all_data.drop(['item_price_mean', 'item_count'], axis=1)
all_data = downcast(all_data, False)
all_data.info()
# 出力
<class 'pandas.core.frame.DataFrame'>
Int64Index: 9904582 entries, 1122386 to 11026967
Data columns (total 24 columns):
# Column Dtype
--- ------ -----
0 date_block_num int8
1 shop_id int8
2 item_id int16
3 item_cnt_month int8
4 city int8
5 item_category_id int8
6 category int8
7 item_id_mean_sales float32
8 item_id_city_mean_sales float32
9 item_cnt_month_lag1 int8
10 item_cnt_month_lag2 int8
11 item_cnt_month_lag3 int8
12 item_count_lag1 int8
13 item_count_lag2 int8
14 item_count_lag3 int8
15 item_price_mean_lag1 float32
16 item_price_mean_lag2 float32
17 item_price_mean_lag3 float32
18 item_cnt_month_lag_mean float32
19 lag_grad1 float32
20 lag_grad2 float32
21 brand_new int8
22 duration_after_first_sale int8
23 month int8
dtypes: float32(8), int16(1), int8(15)
memory usage: 538.4 MB
モデル実装
今回はLightGBMを使って予測をしていきます。
ハイパーパラメータのチューニングはOptunaを使っています。
まずはtrain、validation、testデータに分けます。
X_train = all_data[all_data.date_block_num < 33].drop(['item_cnt_month'], axis=1)
y_train = all_data[all_data.date_block_num < 33]['item_cnt_month']
X_val = all_data[all_data.date_block_num == 33].drop(['item_cnt_month'], axis=1)
y_val = all_data[all_data.date_block_num == 33]['item_cnt_month']
X_test = all_data[all_data.date_block_num == 34].drop(['item_cnt_month'], axis=1)
評価されるのはRMSEなので、算出するための関数を定義します。
また、特徴量重要度を可視化させる準備もしておきます。
def rmse(y, y_pred):
return np.sqrt(mean_squared_error(y, y_pred))
def plot_features(booster, figsize):
fig, ax = plt.subplots(1,1,figsize=figsize)
return plot_importance(booster=booster, ax=ax)
LightGBMの準備をします。
lgb_train = lgb.Dataset(X_train, y_train, free_raw_data=False)
lgb_eval = lgb.Dataset(X_val, y_val, reference=lgb_train, free_raw_data=False)
import optuna
from sklearn.metrics import log_loss
def objective(trial):
params = {
'objective': 'regression',
'metric': 'rmse',
'verbosity': -1,
'boosting_type': 'gbdt',
'lambda_l1': trial.suggest_float("lambda_l1", 1e-8, 10.0, log=True),
'lambda_l2': trial.suggest_float("lambda_l2", 1e-8, 10.0, log=True),
'feature_fraction': trial.suggest_float("feature_fraction", 0.4, 1.0),
'bagging_fraction': trial.suggest_float("bagging_fraction", 0.4, 1.0),
'bagging_freq': trial.suggest_int("bagging_freq", 1, 7),
'min_child_samples': trial.suggest_int("min_child_samples", 5, 100),
'learning_rate': 0.05,
'num_leaves': trial.suggest_int('num_leaves', 2, 256),
}
cat_features = ['shop_id', 'item_category_id', 'category', 'date_block_num']
model = lgb.train(params,
train_set=lgb_train,
valid_sets=[lgb_train, lgb_eval],
early_stopping_rounds=15,
categorical_feature=cat_features,
verbose_eval=1)
y_pred = model.predict(X_val, num_iteration=model.best_iteration)
accuracy = rmse(y_val, y_pred)
return accuracy
学習させます。
study = optuna.create_study(direction='minimize')
study.optimize(objective, n_trials=5)
print('Number of finished trials:', len(study.trials))
print('Best trial:', study.best_trial.params)
# 出力(一部省略)
Number of finished trials: 5
Best trial: {'lambda_l1': 8.534179512307704, 'lambda_l2': 5.549945256079241e-06, 'feature_fraction': 0.6357223930020901, 'bagging_fraction': 0.524237510163727, 'bagging_freq': 4, 'min_child_samples': 29, 'num_leaves': 218}
パラメータをベストなもので修正して追加していきます。
best_params = study.best_trial.params
x = {"objective": "regression",
"metric" : "rmse",
"verbosity": -1,
"boosting_type": "gbdt"}
best_params.update(x)
best_params
# 出力
{'lambda_l1': 8.534179512307704,
'lambda_l2': 5.549945256079241e-06,
'feature_fraction': 0.6357223930020901,
'bagging_fraction': 0.524237510163727,
'bagging_freq': 4,
'min_child_samples': 29,
'num_leaves': 218,
'objective': 'regression',
'metric': 'rmse',
'verbosity': -1,
'boosting_type': 'gbdt'}
evals_result = {}
model = lgb.train(best_params,
lgb_train,
valid_sets=[lgb_train,lgb_eval],
evals_result=evals_result,
early_stopping_rounds=30, # 20
verbose_eval=1,
)
y_pred = model.predict(X_val)
rmse(y_val, y_pred)
# 出力
0.47423599685001055
特徴量重要度を確認します。
plot_features(model, (10,14))
shop_idやラグ特徴量の重要度が高いと思われます。
ax = lgb.plot_metric(evals_result, figsize=(10, 5))
# 予測モデルの確認し、提出
y_test = model.predict(X_test).clip(0, 20)
submission = pd.DataFrame({
"ID": test.index,
"item_cnt_month": y_test
})
submission.head(10)
submission.to_csv('submission.csv', index=False)
LeaderBoardでのscoreは1.23316でした。
RMSEとしては低くない値かなとは思ったのですが、予想以上にscoreが伸びませんでした。
特徴量の生成の仕方などを見直す必要があるかもしれません。