はじめに
Kaggleは、英語のページしかありません。そこで、日本語で読みたい方向けに記事を作成しました。
Predict Future Salesの概要の説明が前のブログ記事で済みました。
実際にGoogle Colabを動かしてみましょう！
Google Colaboratoryについては知りたい方は、以下のブログ記事を参考にしてください。
>>Google Colaboratoryとは？ いつできた？ mount、ファイルの読み込み等の使い方
>>Google Colaboratoryよく使う便利なショートカットキー
こちらの方のコードを参考に解説していきます。
https://www.kaggle.com/code/ekrembayar/store-sales-ts-forecasting-a-comprehensive-guide
データのインポート
import datetime
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from xgboost import XGBRegressor
from xgboost import plot_importance
import warnings
warnings.filterwarnings('ignore')
%matplotlib inline
from sklearn.metrics import mean_squared_error
from sklearn.linear_model import LinearRegression
from sklearn.neighbors import KNeighborsRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn.preprocessing import StandardScaler, MinMaxScaler
sns.set()
pd.set_option('display.float_format', lambda x: '%.2f' % x)
pd.options.mode.use_inf_as_na = True
sales = pd.read_csv('/content/drive/MyDrive/kaggle/Predict Future Sales/sales_train.csv', parse_dates=['date'],
dtype={'date': 'str', 'date_block_num': 'int32', 'shop_id': 'int32',
'item_id': 'int32', 'item_price': 'float32', 'item_cnt_day': 'int32'})
test = pd.read_csv('/content/drive/MyDrive/kaggle/Predict Future Sales/test.csv', dtype={'ID': 'int32', 'shop_id': 'int32',
'item_id': 'int32'})
items = pd.read_csv('/content/drive/MyDrive/kaggle/Predict Future Sales/items.csv', dtype={'item_name': 'str', 'item_id': 'int32',
'item_category_id': 'int32'})
item_categories = pd.read_csv('/content/drive/MyDrive/kaggle/Predict Future Sales/item_categories.csv',
dtype={'item_category_name': 'str', 'item_category_id': 'int32'})
shops = pd.read_csv('/content/drive/MyDrive/kaggle/Predict Future Sales/shops.csv', dtype={'shop_name': 'str', 'shop_id': 'int32'})
train = sales.join(items, on='item_id', rsuffix='_').join(shops, on='shop_id', rsuffix='_').join(item_categories, on='item_category_id', rsuffix='_').drop(['item_id_', 'shop_id_', 'item_category_id_'], axis=1)
trainデータの確認
trainの統計情報を見てましょう。
次にtrainデータの期間を見てみます。
print('trainの開始日:%s' % train['date'].min().date())
print('trainの最終日:%s' % train['date'].max().date())
データリークについては、testに登場する「shop_id」「item_id」だけを使うことにします。
test_shop_ids = test['shop_id'].unique()
test_item_ids = test['item_id'].unique()
lk_train = train[train['shop_id'].isin(test_shop_ids)]
lk_train = lk_train[lk_train['item_id'].isin(test_item_ids)]データリーク前と後のデータを見てみます。
print('before date leakage:', train.shape[0])
print('after date leakage:', lk_train.shape[0])
trainの中で価格が入っているアイテムだけを抽出します。
train = train.query('item_price > 0')データの前処理
不要なカラムを削除します。trainの中の-nameカラムは、-idカラムに代替できるで、使用しません(item_nameはitem_idに代替できる)。
このコンペは翌月の全商品、全店舗の総売上高を予測するよう求められていますが、データは日ごとに与えられているので、不要な列を削除して月ごとにデータを集計してみましょう。
train_monthly = lk_train[['date', 'date_block_num', 'shop_id', 'item_category_id', 'item_id', 'item_price', 'item_cnt_day']]月ごとにグループ化して、特徴量を集計します。
train_monthly = train_monthly.sort_values('date').groupby(['date_block_num', 'shop_id', 'item_category_id', 'item_id'], as_index=False)
train_monthly = train_monthly.agg({'item_price':['sum', 'mean'], 'item_cnt_day':['sum', 'mean','count']})
train_monthly.columns = ['date_block_num', 'shop_id', 'item_category_id', 'item_id', 'item_price', 'mean_item_price', 'item_cnt', 'mean_item_cnt', 'transactions']
次に、欠損値の処理をするために、‘date_block_num’,’shop_id’,’item_id’の欠損値のデータフレームを作成します。
for文では、‘date_block_num’のユニーク値である34回繰り返し処理を行います。
shop_ids = train_monthly['shop_id'].unique()
item_ids = train_monthly['item_id'].unique()
empty_df = []
for i in range(34):
for shop in shop_ids:
for item in item_ids:
empty_df.append([i, shop, item])
empty_df = pd.DataFrame(empty_df, columns=['date_block_num','shop_id','item_id'])
train_monthlyとempty_dfを結合し、欠損値を０で埋めます。
train_monthly = pd.merge(empty_df, train_monthly, on=['date_block_num','shop_id','item_id'], how='left')
train_monthly.fillna(0, inplace=True)
時間軸からみた特徴も抽出してみましょう。
まずは、train_monthlyに’year’と’month’のカラムを追加します。
train_monthly['year'] = train_monthly['date_block_num'].apply(lambda x: ((x//12) + 2013))
train_monthly['month'] = train_monthly['date_block_num'].apply(lambda x: (x % 12))月間の平均と合計の売上推移を折れ線グラフで見てみます。
gp_month_mean = train_monthly.groupby(['month'], as_index=False)['item_cnt'].mean()
gp_month_sum = train_monthly.groupby(['month'], as_index=False)['item_cnt'].sum()
gp_category_mean = train_monthly.groupby(['item_category_id'], as_index=False)['item_cnt'].mean()
gp_category_sum = train_monthly.groupby(['item_category_id'], as_index=False)['item_cnt'].sum()
gp_shop_mean = train_monthly.groupby(['shop_id'], as_index=False)['item_cnt'].mean()
gp_shop_sum = train_monthly.groupby(['shop_id'], as_index=False)['item_cnt'].sum()
f,axes = plt.subplots(2,1,figsize=(22,10),sharex=True)
sns.lineplot(x='month',y='item_cnt',data=gp_month_mean,ax=axes[0]).set_title('Monthly mean')
sns.lineplot(x='month',y='item_cnt',data=gp_month_sum,ax=axes[1]).set_title('Monthly sum')
plt.show()
年末に向けて販売個数（平均値）が増加する傾向が見られます。
では、どのアイテムカテゴリーが売れているか見てみます。
f,axes = plt.subplots(2,1,figsize=(22,10),sharex=True)
sns.barplot(x='item_category_id',y='item_cnt',data=gp_category_mean,ax=axes[0],palette='rocket').set_title('Monthly mean')
sns.barplot(x='item_category_id',y='item_cnt',data=gp_category_sum,ax=axes[1],palette='rocket').set_title('Monthly sum')
plt.show()
一部のカテゴリーだけが、販売数の大部分を占めています。
それでは、どの店舗が売れているかみてましょう。
f,axes = plt.subplots(2,1,figsize=(22,10),sharex=True)
sns.barplot(x='shop_id',y='item_cnt',data=gp_shop_mean,ax=axes[0],palette='rocket').set_title('Monthly mean')
sns.barplot(x='shop_id',y='item_cnt',data=gp_shop_sum,ax=axes[1],palette='rocket').set_title('Monthly sum')
plt.show()
ほとんどの店が同じような販売率だが、なかでも3店はかなり高いです。
外れ値を確認してみましょう。
sns.jointplot(x='item_cnt',y='item_price',data=train_monthly,height=8)
plt.show()
sns.jointplot(x="item_cnt", y="transactions", data=train_monthly, height=8)
plt.show()
item_cntの特徴分布を見てみます。
plt.subplots(figsize=(22,8))
sns.boxplot(train_monthly['item_cnt'])
plt.show()
確認できた外れ値の除去をしましょう。
item_cnt” > 20 and < 0, “item_price” >= 400000は異常値として扱うので、削除する。
train_monthly = train_monthly.query('item_cnt >= 0 and item_cnt <= 20 and item_price < 400000')train_monthlyに‘item_cnt_month’を作成します。
train_monthly['item_cnt_month'] = train_monthly.sort_values('date_block_num').groupby(['shop_id', 'item_id'])['item_cnt'].shift(-1)
train_monthlyに‘item_price_unit’（単品価格）を作成します。
train_monthly['item_price_unit'] = train_monthly['item_price'] // train_monthly['item_cnt']
train_monthly['item_price_unit'].fillna(0, inplace=True)item_idに対して最小値、最大値のそれぞれのitem_priceを作成し、train_monthlyとマージします。
gp_item_price = train_monthly.sort_values('date_block_num').groupby(['item_id'],as_index=False).agg({'item_price':[np.min,np.max]})
gp_item_price.columns = ['item_id','hist_min_item_price','hist_max_item_price']
train_monthly = pd.merge(train_monthly,gp_item_price,how='left',on='item_id')各アイテムの価格が過去の（最安値/最高値）価格からどれだけ変化したかというカラムも作成します。
train_monthly['price_increase'] = train_monthly['item_price'] - train_monthly['hist_min_item_price']
train_monthly['price_decrease'] = train_monthly['hist_max_item_price'] - train_monthly['item_price']過去の値から統計的な値を算出するローリングウィンドウ法(window = 3 months)を使用します。
# Min value
f_min = lambda x: x.rolling(window=3, min_periods=1).min()
# Max value
f_max = lambda x: x.rolling(window=3, min_periods=1).max()
# Mean value
f_mean = lambda x: x.rolling(window=3, min_periods=1).mean()
# Standard deviation
f_std = lambda x: x.rolling(window=3, min_periods=1).std()
function_list = [f_min, f_max, f_mean, f_std]
function_name = ['min', 'max', 'mean', 'std']
for i in range(len(function_list)):
train_monthly[('item_cnt_%s' % function_name[i])] = train_monthly.sort_values('date_block_num').groupby(['shop_id', 'item_category_id', 'item_id'])['item_cnt'].apply(function_list[i])
# Fill the empty std features with 0
train_monthly['item_cnt_std'].fillna(0, inplace=True)ラグ特徴量をみます。
これは、時系列予測問題を教師あり学習問題に変換する古典的な方法です。
最も単純なアプローチは，現在の時間をｔとすると、前の時間(t-1)の値から次の時間(t+1)の値を予測することである。
時刻tの値は、時刻t-1の値に大きく影響されるのです。過去の値をラグといいますから、t-1がラグ1、t-2がラグ2というようになります。
lag_list = [1, 2, 3]
for lag in lag_list:
ft_name = ('item_cnt_shifted%s' % lag)
train_monthly[ft_name] = train_monthly.sort_values('date_block_num').groupby(['shop_id', 'item_category_id', 'item_id'])['item_cnt'].shift(lag)
# Fill the empty shifted features with 0
train_monthly[ft_name].fillna(0, inplace=True)アイテム販売数推移もカラムにいれましょう。
train_monthly['item_trend'] = train_monthly['item_cnt']
for lag in lag_list:
ft_name = ('item_cnt_shifted%s' % lag)
train_monthly['item_trend'] -= train_monthly[ft_name]
train_monthly['item_trend'] /= len(lag_list) + 1
train/validation/testを作成
trainセットは最初の3~27ブロック、validationは最後の5ブロック(28~32)、testは33ブロックを使用します。
最初の3ヶ月を省いているのは、3ヶ月のウィンドウを使って特徴量を生成しているからです。
train_set = train_monthly.query('date_block_num >= 3 and date_block_num < 28').copy()
validation_set = train_monthly.query('date_block_num >= 28 and date_block_num < 33').copy()
test_set = train_monthly.query('date_block_num == 33').copy()
train_set.dropna(subset=['item_cnt_month'], inplace=True)
validation_set.dropna(subset=['item_cnt_month'], inplace=True)
train_set.dropna(inplace=True)
validation_set.dropna(inplace=True)
print('Train set records:', train_set.shape[0])
print('Validation set records:', validation_set.shape[0])
print('Test set records:', test_set.shape[0])
print('Train set records: %s (%.f%% of complete data)' % (train_set.shape[0], ((train_set.shape[0]/train_monthly.shape[0])*100)))
print('Validation set records: %s (%.f%% of complete data)' % (validation_set.shape[0], ((validation_set.shape[0]/train_monthly.shape[0])*100)))
カテゴリ変数をそれぞれ、平均値でエンコーディングします。
# Shop 平均値でのエンコーディング
gp_shop_mean = train_set.groupby(['shop_id']).agg({'item_cnt_month': ['mean']})
gp_shop_mean.columns = ['shop_mean']
gp_shop_mean.reset_index(inplace=True)
# Item 平均値のエンコーディング
gp_item_mean = train_set.groupby(['item_id']).agg({'item_cnt_month': ['mean']})
gp_item_mean.columns = ['item_mean']
gp_item_mean.reset_index(inplace=True)
# Shop と item の平均値のエンコーディング
gp_shop_item_mean = train_set.groupby(['shop_id', 'item_id']).agg({'item_cnt_month': ['mean']})
gp_shop_item_mean.columns = ['shop_item_mean']
gp_shop_item_mean.reset_index(inplace=True)
# Year 平均値でのエンコーディング
gp_year_mean = train_set.groupby(['year']).agg({'item_cnt_month': ['mean']})
gp_year_mean.columns = ['year_mean']
gp_year_mean.reset_index(inplace=True)
# Month 平均値でのエンコーディング
gp_month_mean = train_set.groupby(['month']).agg({'item_cnt_month': ['mean']})
gp_month_mean.columns = ['month_mean']
gp_month_mean.reset_index(inplace=True)
# 上記をtrainセットに加える
train_set = pd.merge(train_set, gp_shop_mean, on=['shop_id'], how='left')
train_set = pd.merge(train_set, gp_item_mean, on=['item_id'], how='left')
train_set = pd.merge(train_set, gp_shop_item_mean, on=['shop_id', 'item_id'], how='left')
train_set = pd.merge(train_set, gp_year_mean, on=['year'], how='left')
train_set = pd.merge(train_set, gp_month_mean, on=['month'], how='left')
# 上記をvalidationセットに加える
validation_set = pd.merge(validation_set, gp_shop_mean, on=['shop_id'], how='left')
validation_set = pd.merge(validation_set, gp_item_mean, on=['item_id'], how='left')
validation_set = pd.merge(validation_set, gp_shop_item_mean, on=['shop_id', 'item_id'], how='left')
validation_set = pd.merge(validation_set, gp_year_mean, on=['year'], how='left')
validation_set = pd.merge(validation_set, gp_month_mean, on=['month'], how='left')
学習データ（X_train、X_validation）から‘item_cnt_month’, ‘date_block_num’を削除し、目的変数（Y_train、Y_validation）は‘item_cnt_month’とします。
X_train = train_set.drop(['item_cnt_month', 'date_block_num'], axis=1)
Y_train = train_set['item_cnt_month'].astype(int)
X_validation = validation_set.drop(['item_cnt_month', 'date_block_num'], axis=1)
Y_validation = validation_set['item_cnt_month'].astype(int)‘shop_id’, ‘item_id’, ‘year’, ‘month’をint32に変換し、X_train、X_validationに追加します。
int_features = ['shop_id', 'item_id', 'year', 'month']
X_train[int_features] = X_train[int_features].astype('int32')
X_validation[int_features] = X_validation[int_features].astype('int32')testセットを作成します。
latest_records = pd.concat([train_set, validation_set]).drop_duplicates(subset=['shop_id', 'item_id'], keep='last')
X_test = pd.merge(test, latest_records, on=['shop_id', 'item_id'], how='left', suffixes=['', '_'])
X_test['year'] = 2015
X_test['month'] = 9
X_test.drop('item_cnt_month', axis=1, inplace=True)
X_test[int_features] = X_test[int_features].astype('int32')
X_test = X_test[X_train.columns]
欠損値を各項目の中央値で置き換える。
sets = [X_train, X_validation, X_test]
for dataset in sets:
for shop_id in dataset['shop_id'].unique():
for column in dataset.columns:
shop_median = dataset[(dataset['shop_id'] == shop_id)][column].median()
dataset.loc[(dataset[column].isnull()) & (dataset['shop_id'] == shop_id), column] = shop_median
X_test.fillna(X_test.mean(), inplace=True)item_category_idを削除する。
X_train.drop(['item_category_id'], axis=1, inplace=True)
X_validation.drop(['item_category_id'], axis=1, inplace=True)
X_test.drop(['item_category_id'], axis=1, inplace=True)
モデルの構築
XGBoostを作成します。
XGBoostは、高効率、柔軟性、移植性に優れた最適化された分散型勾配ブースティング・ライブラリでKaggleで人気のモデルです。
今回は回帰の問題なので、XGBRegressor（）を用います。
# XGBoost
xgb_features = ['item_cnt','item_cnt_mean', 'item_cnt_std', 'item_cnt_shifted1',
'item_cnt_shifted2', 'item_cnt_shifted3', 'shop_mean',
'shop_item_mean', 'item_trend', 'mean_item_cnt']
xgb_train = X_train[xgb_features]
xgb_val = X_validation[xgb_features]
xgb_test = X_test[xgb_features]xgb_model = XGBRegressor(max_depth=8,
n_estimators=500,
min_child_weight=1000,
colsample_bytree=0.7,
subsample=0.7,
eta=0.3,
seed=0)
xgb_model.fit(xgb_train,
Y_train,
eval_metric='rmse',
eval_set=[(xgb_train,Y_train),(xgb_val,Y_validation)],
verbose=20,
early_stopping_rounds=20)XGBoostの特徴量の重要度をみてみる。
plt.rcParams['figure.figsize'] = (15,6)
plot_importance(xgb_model)
plt.show()
特徴量の重要度としては、‘shop_mean’がダントツ１位で‘shop_item_mean’が２位となっています。
まずは、これをもとに予測し、RMSEのスコアをみてみます。
xgb_train_pred = xgb_model.predict(xgb_train)
xgb_val_pred = xgb_model.predict(xgb_val)
xgb_test_pred = xgb_model.predict(xgb_test)print('Train rmse:', np.sqrt(mean_squared_error(Y_train, xgb_train_pred)))
print('Validation rmse:', np.sqrt(mean_squared_error(Y_validation, xgb_val_pred)))
次は、ランダムフォレスト（特徴量として用いるのは一部）を使用してみます。
rf_features = ['shop_id', 'item_id', 'item_cnt', 'transactions', 'year',
'item_cnt_mean', 'item_cnt_std', 'item_cnt_shifted1',
'shop_mean', 'item_mean', 'item_trend', 'mean_item_cnt']
rf_train = X_train[rf_features]
rf_val = X_validation[rf_features]
rf_test = X_test[rf_features]rf_model = RandomForestRegressor(n_estimators=50, max_depth=7, random_state=0, n_jobs=-1)
rf_model.fit(rf_train, Y_train)rf_train_pred = rf_model.predict(rf_train)
rf_val_pred = rf_model.predict(rf_val)
rf_test_pred = rf_model.predict(rf_test)print('Train rmse:', np.sqrt(mean_squared_error(Y_train, rf_train_pred)))
print('Validation rmse:', np.sqrt(mean_squared_error(Y_validation, rf_val_pred)))
lr_features = ['item_cnt', 'item_cnt_shifted1', 'item_trend', 'mean_item_cnt', 'shop_mean']
lr_train = X_train[lr_features]
lr_val = X_validation[lr_features]
lr_test = X_test[lr_features]そして、線形モデル（特徴量として用いるのは一部）も使用してみます。
線形モデルなので、特徴量をスケーリングします（単位や大きさを揃えないと特定の特徴量だけに反応してしまうモデルになってしますので）。
lr_scaler = MinMaxScaler()
lr_scaler.fit(lr_train)
lr_train = lr_scaler.transform(lr_train)
lr_val = lr_scaler.transform(lr_val)
lr_test = lr_scaler.transform(lr_test)lr_model = LinearRegression(n_jobs=-1)
lr_model.fit(lr_train, Y_train)lr_train_pred = lr_model.predict(lr_train)
lr_val_pred = lr_model.predict(lr_val)
lr_test_pred = lr_model.predict(lr_test)print('Train rmse:', np.sqrt(mean_squared_error(Y_train, lr_train_pred)))
print('Validation rmse:', np.sqrt(mean_squared_error(Y_validation, lr_val_pred)))
クラスタリングモデルとして、K近傍法を用います。
K近傍法とは、レコード間の距離をそれらの特徴量の値の差を用いて定義し、その距離が最も近いK個のレコードの目的変数から回帰または分類を行います。
例えば、グループA・グループBがあり、それぞれのグループに所属している人の特性がわかっているとします。そこに、新しい人がくるとその人の特性を調べて、近いほうのグループに割り当てるというようなものです。
knn_features = ['item_cnt', 'item_cnt_mean', 'item_cnt_std', 'item_cnt_shifted1',
'item_cnt_shifted2', 'shop_mean', 'shop_item_mean',
'item_trend', 'mean_item_cnt']
X_train_sampled = X_train[:100000]
Y_train_sampled = Y_train[:100000]
knn_train = X_train_sampled[knn_features]
knn_val = X_validation[knn_features]
knn_test = X_test[knn_features]こちらも特徴量のスケーリングを行います。
knn_scaler = MinMaxScaler()
knn_scaler.fit(knn_train)
knn_train = knn_scaler.transform(knn_train)
knn_val = knn_scaler.transform(knn_val)
knn_test = knn_scaler.transform(knn_test)knn_model = KNeighborsRegressor(n_neighbors=9, leaf_size=13, n_jobs=-1)
knn_model.fit(knn_train, Y_train_sampled)knn_train_pred = knn_model.predict(knn_train)
knn_val_pred = knn_model.predict(knn_val)
knn_test_pred = knn_model.predict(knn_test)print('Train rmse:', np.sqrt(mean_squared_error(Y_train_sampled, knn_train_pred)))
print('Validation rmse:', np.sqrt(mean_squared_error(Y_validation, knn_val_pred)))
スタッキング
第1レベルモデルからの予測値を用いて新しいデータセットを作成します。
ここでは簡単なアンサンブル手法を使います。
第1レベルモデルの予測を第2レベルモデルの入力として使います。この方法では第2レベルモデルは基本的に第1レベルモデルの予測を特徴として使い、どこに重きを置くべきかを学習します。
このテクニックを使うには、第1レベルモデルを使ってテストセットで予測を行い、第2レベルモデルでそれを使えるようにする必要があります。
また、第2レベルモデルに、余分な特徴（第1レベルモデルの予測）を含む完全な検証セットを渡して、解を見つけるための作業をもう少しさせることもできます。
first_level = pd.DataFrame(xgb_val_pred, columns=['xgbm'])
first_level['random_forest'] = rf_val_pred
first_level['linear_regression'] = lr_val_pred
first_level['knn'] = knn_val_pred
first_level['label'] = Y_validation.values
first_level.head()
first_level_test = pd.DataFrame(xgb_test_pred, columns=['xgbm'])
first_level_test['random_forest'] = rf_test_pred
first_level_test['linear_regression'] = lr_test_pred
first_level_test['knn'] = knn_test_pred
first_level_test.head()
第2レベルモデルとして線形回帰を用います。
これは、他のモデルを組み合わせて、全体としてより良い予測をすることを期待するモデルです。
meta_model = LinearRegression(n_jobs=-1)第1レベルモデルの予測値を特徴量として、validationセットで学習する。
first_level.drop('label', axis=1, inplace=True)
meta_model.fit(first_level, Y_validation)第1レベルモデルの予測値を特徴量として、testセットで予測を行う。
ensemble_pred = meta_model.predict(first_level)
final_predictions = meta_model.predict(first_level_test)print('Train rmse:', np.sqrt(mean_squared_error(ensemble_pred, Y_validation)))
ファイルの提出
prediction_df = pd.DataFrame(test['ID'], columns=['ID'])
prediction_df['item_cnt_month'] = final_predictions.clip(0., 20.)
prediction_df.to_csv('submission.csv', index=False)
prediction_df.head()
時系列分析でよく使うProphetについてもまとめましたので、こちらもご参考にしてください。
>>【完全版】Prophetー時系列分析の基本から外部変数追加、holidaysの追加、ハイパーパラメータの調整まで(python)
Kaggleで悩んだら
「Kaggle で勝つデータ分析の技術」
以下の書籍は、Kaggleを始める方には本当にオススメの書籍です。Kaggleでわからないことや悩んだことがあった方は、購入を検討してみください。
リンク
本だけでは物足りないという方は、動画のプラットフォームで学ぶこともオススメです。興味がございましたら、以下の無料のオンライン説明会に参加してみてはいかがでしょうか。
