kaggle

【Kaggle】Predict Future Salesをやってみた

はじめに

以前記事にしていますが、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が伸びませんでした。

特徴量の生成の仕方などを見直す必要があるかもしれません。