CRUNKYおいしい

CRUNKY食べながらのんびり書いてる

Mercari Price Challenge -機械学習を使ったメルカリの価格予測 Ridge回帰 LightGBM

ポケモンデータ解析に続いて、またKaggleでのデータ解析ネタです。 今回の解析テーマはフリマアプリのメルカリです。提供されているのはアメリカで行われた商品毎の取引データです。データ内には出品者側が登録した商品名、商品コンディションランク、商品カテゴリ、ブランド名、価格、配送料負担サイド、商品概要が含まれています。自分もメルカリをたまに使用しますが、ユーザはカテゴリ欄や商品概要などは自由に選択したり書いているので、統率の取れているデータになっていない部分もあります。 現在賞金のかかったコンペが行われていて、一番精度の高い価格予測モデルを作成した人は10万ドル(約1,000万円)が貰えるそうです! まぁレベルが高すぎて、僕には関係ない話なのでコンペは無視してデータ解析の勉強をしていきたいと思います。 参加者の賢い人たちがコードを公開してくれているので、参考にしながら機械学習の方法を学んでいきます。 こちらの方を参考にします。(https://www.kaggle.com/iamprateek/submission-to-mercari-price-suggestion-challenge/notebook)
今回はけっこマニアックな内容になってしまいました。。出てくる主な学習ツールは、Ridge回帰、LightGBMです。

さてメルカリから提供されているデータは148万件もあります。今回148万件ものデータは処理するのが大変で、僕のPCでは処理計算に1時間ほどかかってしまうこともあり今回はサクサク進めていきたいので1万件のみを扱います。

まずはデータをざっくり俯瞰します。

import pandas as pd
train = pd.read_csv("train.csv", nrows=10000)
print ("train.shape: " + str(train.shape))
train.head()
train.shape: (10000, 8)

データにあるのはこんな情報ですね。以下写真の赤枠ですね。写真の取り方とか出品者のいいね数なども実際は影響しそうですが、今回はこれらの情報から価格が予想できるような学習モデルを導くのが目的です。

import matplotlib.pyplot as plt
%matplotlib inline
plt.figure(figsize=(10, 6))
plt.hist(train["price"])
plt.xlabel("price[$]")
plt.ylabel("count")
plt.title("Price Histogram")
plt.show()

ほとんどが200$以内の価格帯で取引されてますね。もうちょい拡大しますか。

plt.figure(figsize=(10, 6))
plt.hist(train["price"], bins=500)
plt.xlim(0, 200)
plt.xlabel("price[$]")
plt.ylabel("count")
plt.title("Price Histogram")
plt.show()

10~20$くらいが一番頻度の高い取引みたいです。1000~2000円程度ですので、実感的にも確かにそうだろうなって感じです。

次はカテゴリ別の価格分布を見てみたいです。カテゴリ名欄を分割して新カラムを作成し視覚化してみます。

import time
import seaborn as sns
def split_cat(text):
    try: return text.split("/")
    except: return ("No Label", "No Label", "No Label")
train["general_cat"], train["sub_cat1"], train["sub_cat2"] = \
    zip(*train['category_name'].apply(lambda x: split_cat(x)))
train.drop("category_name", axis=1)
start_time = time.time()
plt.figure(figsize=(16, 8))
ax = sns.violinplot(x="general_cat", y="price", data=train, inner=None)
ax = sns.swarmplot(x="general_cat", y="price", data=train, edgecolor="gray", hue="sub_cat1")
plt.xticks(rotation=30)
plt.ylim(0, 200)
plt.legend(loc=9, bbox_to_anchor=(0.5, -0.25), ncol=5)
plt.show()
print('Showing graph took {} secs.'.format(time.time() - start_time))

Showing graph took 1226.6749708652496 secs.

グラフを表示するだけに20分以上もかかってしまいました。。グラフをざっくりと考察してみると、キッズ用品や美容系は安値側にどっしりしたバイオリン図になっており比較的安価層に多く偏っていそうなことだったり、逆にElectronicsは安値から高値までバラついてそうなことがわかりまs。また安価帯の値段設定は, 10ドル、25ドル、40ドル、75ドルあたりなどに分布されてるものが多く、あまり連続的な分布ではなく階段状の分布がありそうですね。まぁ日本でも999円, 1999円, 9999円とか設定されやすい価格ってありますよね。そんな感じかなという印象です。

ちなみに一番高い価格でやり取りされた商品はどんな商品だったんでしょうか??見てみましょう。

train.sort_values(by="price", ascending=False).head(1)

Chanel Classic Flag Bag medium Caviar L 1506$
1500ドルですから約15万。まぁそこまで高くはなかったです。 シャネルのショルダーバッグみたいです。検索してみたところこれですかね

データ俯瞰はこれぐらいにして、本題である機械学習ツールを使った価格学習に入っていきます。 まずは機械学習のお決まり、学習器がちゃんと処理できるようにデータ整理です。

空欄になっているカラムについては"missing"という値を入れます。

def handle_missing_inplace(dataset):
    dataset['general_cat'].fillna(value='missing', inplace=True)
    dataset['sub_cat1'].fillna(value='missing', inplace=True)
    dataset['sub_cat2'].fillna(value='missing', inplace=True)
    dataset['brand_name'].fillna(value='missing', inplace=True)
    dataset['item_description'].fillna(value='missing', inplace=True)
handle_missing_inplace(train)
train['brand_name'].value_counts().head()
missing              4261
Nike                  350
PINK                  348
Victoria's Secret     322
LuLaRoe               201
Name: brand_name, dtype: int64

ブランド名のところは空白が多くて40%以上ものの欄にmissingが入りました。
次はブランド名やカテゴリ名なども出現頻度が極端に少ないものは学習しても仕方ないので、空欄と同じ扱いのmissing値に置換してしまいます。

def cutting(dataset):
    pop_brand = dataset['brand_name'].value_counts().loc[lambda x: x.index != 'missing'].index[:750]
    dataset.loc[~dataset['brand_name'].isin(pop_brand), 'brand_name'] = 'missing'
    pop_category1 = dataset['general_cat'].value_counts().loc[lambda x: x.index != 'missing'].index[:450]
    pop_category2 = dataset['sub_cat1'].value_counts().loc[lambda x: x.index != 'missing'].index[:450]
    pop_category3 = dataset['sub_cat2'].value_counts().loc[lambda x: x.index != 'missing'].index[:450]
    dataset.loc[~dataset['general_cat'].isin(pop_category1), 'general_cat'] = 'missing'
    dataset.loc[~dataset['sub_cat1'].isin(pop_category2), 'sub_cat1'] = 'missing'
    dataset.loc[~dataset['sub_cat2'].isin(pop_category3), 'sub_cat2'] = 'missing'
cutting(train)
train['brand_name'].value_counts().head()
missing              4272
Nike                  350
PINK                  348
Victoria's Secret     322
LuLaRoe               201
Name: brand_name, dtype: int64

ちょっとmissingが増えました。(参考コードではデータ量が多いのでこの作業の効力は大きいかもしれませんが、今回はデータ量はそれほど多くないので必要性はないかもしれません)
次はカテゴリ名と商品概要はテキストデータが入っており、次段階の解析で単語抽出を行いたいためデータ型をカテゴリ型に変換しておきます。

def to_categorical(dataset):
    dataset['general_cat'] = dataset['general_cat'].astype('category')
    dataset['sub_cat1'] = dataset['sub_cat1'].astype('category')
    dataset['sub_cat2'] = dataset['sub_cat2'].astype('category')
    dataset['item_condition_id'] = dataset['item_condition_id'].astype('category')
to_categorical(train)
train.dtypes
train_id                int64
name                   object
item_condition_id    category
category_name          object
brand_name             object
price                 float64
shipping                int64
item_description       object
general_cat          category
sub_cat1             category
sub_cat2             category
dtype: object

データの整理は終わりました。ここから少し難しくなります。
今回のデータ処理に対して回帰や分類をしようとする際に障壁となるのが、文字などの自然言語データです。これらを数学的に処理するのがひとつの壁になります。特に今回はユーザーがそれぞれの感覚で自由に書いているような文字列なので大変です。 このような自然言語データには単語の組み合わせやその傾向を見て、あぁこういうことを言ってるんだなとコンピュータに気づかせてあげるようにしなければなりません。そこでそういった処理に向いている便利な各種機械学習ライブラリを使っていきます。

まずは商品名に対してです。Scikit-learnという機械学習ライブラリにCountVectorizerというものがあり、これは単語の出現頻度を出すために使うようなクラスでこれを使ってみたいと思います。具体的には、fit関数に渡された文字列データに含まれる単語を単語毎にその出現回数をカウントして、それをベクトルとして表現してくれるクラスだそうです。 パラメータの min_dfは整数で設定し、学習データ全体で現れたサンプル数がmin_df未満の単語については特徴ベクトルからは除外する、というものです。

import numpy as np
from sklearn.feature_extraction.text import CountVectorizer
cv = CountVectorizer(min_df=5)
X_name = cv.fit_transform(train['name'])
print (X_name.shape)
occ = np.asarray(X_name.sum(axis=0)).ravel().tolist()
counts_cv = pd.DataFrame({'term': cv.get_feature_names(), 'occurrences': occ})
counts_cv.sort_values(by='occurrences', ascending=False).head(10)
(10000, 1508)

上記Outputは10,000個のテキストデータに対して1,508個のユニークな単語が含まれていたこと、出現数は"pink"という単語が563回で1番多く出てきたことを示します。5番目に多く出てきてる"lularoe"というのが聞きなれなかったので調べてみると女性用のファッションブランドぽいですね。このブランドの品が多く出品されているんだなとわかります。こういうものを使えば今流行っている服はなんだとかトレンドなんかも見えてきそうです。

同様にカテゴリ名に対しても同じ処理を行います。

cv = CountVectorizer(min_df=5)
combine_category = [train["general_cat"], train["sub_cat1"], train["sub_cat2"]]
X_category1 = cv.fit_transform(train['general_cat'])
print ("----general_cat----")
print (X_category1.shape)
occ = np.asarray(X_category1.sum(axis=0)).ravel().tolist()
counts_cv = pd.DataFrame({'term': cv.get_feature_names(), 'occurrences': occ})
print (counts_cv.sort_values(by='occurrences', ascending=False).head())
X_category2 = cv.fit_transform(train['sub_cat1'])
print ("----sub_cat1----")
print (X_category2.shape)
occ = np.asarray(X_category2.sum(axis=0)).ravel().tolist()
counts_cv = pd.DataFrame({'term': cv.get_feature_names(), 'occurrences': occ})
print (counts_cv.sort_values(by='occurrences', ascending=False).head())
X_category3 = cv.fit_transform(train['sub_cat2'])
print ("----sub_cat2----")
print (X_category3.shape)
occ = np.asarray(X_category3.sum(axis=0)).ravel().tolist()
counts_cv = pd.DataFrame({'term': cv.get_feature_names(), 'occurrences': occ})
print (counts_cv.sort_values(by='occurrences', ascending=False).head())
----general_cat----
(10000, 14)
    occurrences         term
13         4381        women
0          1508       beauty
5          1159         kids
2           816  electronics
7           642          men
----sub_cat1----
(10000, 109)
     occurrences         term
62           918       makeup
3            911  accessories
6            870      apparel
100          850         tops
9            849     athletic
----sub_cat2----
(10000, 317)
     occurrences         term
248          634       shirts
168          442     leggings
206          439        pants
2            438  accessories
288          396       tights

次は、商品概要が記入されているテキストデータの処理を行います。このデータは商品名やカテゴリ名とは違って多くの単語で構成される文章です。CountVectorizerは単語の出現数をカウントするのみでしたが、TfidfVectorizerというクラスを使えば各単語の出現頻度と希少性を計算して”重要そうな”単語を見つけ出すことができます。パラメータのmax_featuresは tfidf値を表示する単語数を指定します。tfidf値が高い順にとっていきます。ngramは複数単語の連続にも対応できます。

from sklearn.feature_extraction.text import TfidfVectorizer
tv = TfidfVectorizer(max_features=1000,
                         ngram_range=(1, 3),
                         stop_words='english')
X_description = tv.fit_transform(train['item_description'])
weights = np.asarray(X_description.mean(axis=0)).ravel().tolist()
weights_df = pd.DataFrame({'term': tv.get_feature_names(), 'weight': weights})
weights_df.sort_values(by='weight', ascending=False).head(10)

最もtfidf値が高いものが"description"となっていますが、これは"Not description yet"と"概要記載なし"のものが影響してそうです。あとは"新しい"だったり"中古"であったりという単語が重み付けされています。

次にブランド名の処理です。ブランド名はある程度限られたものしかないので各ブランドに該当するかどうかを0か1かでラベリングしていく方法が向いてそうです。例えばブランド名がNikeであれば1、Nikeでなければ0みたいな感じですね。このようなラベリングにはLabelBinarizerを使用します。

from sklearn.preprocessing import LabelBinarizer
lb = LabelBinarizer(sparse_output=True)
X_brand = lb.fit_transform(train['brand_name'])
occ = np.asarray(X_brand.sum(axis=0)).ravel().tolist()
counts_cv = pd.DataFrame({'term': lb.classes_, 'occurrences': occ})
counts_cv.sort_values(by='occurrences', ascending=False).head(10)

あとは商品コンディションと配送料負担ですがこれらはすでにラベリングがされているので、ダミー変数を使いましょう。ダミー変数については統計学用語ですのでわからない人は調べてください。またダミー変数を作成すると大量の疎行列(要素のほとんどが0の行列)が発生するためcsr_matrixという疎行列を効率よく処理できるクラスを使用します。(csr_matrixかcsc_matrixかどちらを使うべきかは、後工程で行に対する処理を行うか列に対する処理を行うかで使い分けるべきらしいです。今回は後工程で各行に対して回帰処理を行う予定なのでcsrを使います。)

from scipy.sparse import csr_matrix
X_dummies = csr_matrix(pd.get_dummies(train[['item_condition_id', 'shipping']],
                                          sparse=True).values)
X_dummies
<10000x6 sparse matrix of type '<class 'numpy.int64'>'
    with 14539 stored elements in Compressed Sparse Row format>

以上で10,000件のデータを行列群にすることができました。各10,000件の商品名、カテゴリ名を構成する単語出現数行列、商品概要を構成する単語の重み行列、ブランド名のラベリング、商品コンディション・配送料負担のダミー変数行列がありますのでこれらを学習器に入れていきましょう。それぞれ疎行列ですのでcsrで扱います。

from scipy.sparse import hstack
sparse_merge = hstack([X_dummies, X_description, X_brand, X_category1, X_category2, X_category3, X_name]).tocsr()
sparse_merge.shape
(10000, 3705)

3705個の因子を持ったデータフレームが出来ましたので、これらを回帰・分類学習器にかけることで各因子から予測して価格を算出してみます。ここでのミソとしては複数の学習器の結果を絶妙な割合で掛け合わせることで結果の精度を調整することです。今回使う学習器はRidge回帰とLightGBMという今流行りの勾配ブースティング学習器を使います。掛け合わせ割合は、Ridge: 0.3, LightGBM: 0.35 * 2個です。10,000件のデータは9,000個訓練用、1,000個検証用に分割します。

それぞれの学習器単体での対数平方平均二乗誤差(RMSLE: Root Mean Squared Logarithmic Error)を検証しながら進めていきます。

from sklearn.linear_model import Ridge
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error

X = sparse_merge
y = np.log1p(train["price"])
train_X, test_X, train_y, test_y = train_test_split(X, y, test_size = 0.1, random_state = 144) 

modelR = Ridge(alpha=.5, copy_X=True, fit_intercept=True, max_iter=100,
      normalize=False, random_state=101, solver='auto', tol=0.01)
modelR.fit(train_X, train_y)
predsR = modelR.predict(test_X)

def rmsle(y, y0):
     assert len(y) == len(y0)
     return np.sqrt(np.mean(np.power(np.log1p(y)-np.log1p(y0), 2)))

rmsleR = rmsle(predsR, test_y)
print ("Ridge Regression RMSLE = " + str(rmsleR))
Ridge Regression RMSLE = 0.150454850006
import lightgbm as lgb

train_XL1, valid_XL1, train_yL1, valid_yL1 = train_test_split(train_X, train_y, test_size = 0.1, random_state = 144) 
d_trainL1 = lgb.Dataset(train_XL1, label=train_yL1, max_bin=8192)
d_validL1 = lgb.Dataset(valid_XL1, label=valid_yL1, max_bin=8192)
watchlistL1 = [d_trainL1, d_validL1]
paramsL1 = {
        'learning_rate': 0.65,
        'application': 'regression',
        'max_depth': 3,
        'num_leaves': 60,
        'verbosity': -1,
        'metric': 'RMSE',
        'data_random_seed': 1,
        'bagging_fraction': 0.5,
        'nthread': 4
    }
modelL1 = lgb.train(paramsL1, train_set=d_trainL1, num_boost_round=8000, valid_sets=watchlistL1, \
early_stopping_rounds=5000, verbose_eval=500) 
predsL1 = modelL1.predict(test_X)
rmsleL1 = rmsle(predsL1, test_y)
print ("LightGBM1 RMSLE = " + str(rmsleL1))
Training until validation scores don't improve for 5000 rounds.
[500]   training's rmse: 0.411212   valid_1's rmse: 0.609769
[1000]  training's rmse: 0.35107    valid_1's rmse: 0.628343
[1500]  training's rmse: 0.31196    valid_1's rmse: 0.634885
[2000]  training's rmse: 0.283858   valid_1's rmse: 0.644041
[2500]  training's rmse: 0.260685   valid_1's rmse: 0.650578
[3000]  training's rmse: 0.243261   valid_1's rmse: 0.657922
[3500]  training's rmse: 0.229662   valid_1's rmse: 0.660836
[4000]  training's rmse: 0.217215   valid_1's rmse: 0.663572
[4500]  training's rmse: 0.206279   valid_1's rmse: 0.666011
[5000]  training's rmse: 0.195771   valid_1's rmse: 0.665814
Early stopping, best iteration is:
[96]    training's rmse: 0.528535   valid_1's rmse: 0.590428
LightGBM1 RMSLE = 0.153722873648
train_XL2, valid_XL2, train_yL2, valid_yL2 = train_test_split(train_X, train_y, test_size = 0.1, random_state = 101) 
d_trainL2 = lgb.Dataset(train_XL2, label=train_yL2, max_bin=8192)
d_validL2 = lgb.Dataset(valid_XL2, label=valid_yL2, max_bin=8192)
watchlistL2 = [d_trainL2, d_validL2]
paramsL2 = {
        'learning_rate': 0.85,
        'application': 'regression',
        'max_depth': 3,
        'num_leaves': 140,
        'verbosity': -1,
        'metric': 'RMSE',
        'data_random_seed': 2,
        'bagging_fraction': 1,
        'nthread': 4
    }
modelL2 = lgb.train(paramsL2, train_set=d_trainL2, num_boost_round=5500, valid_sets=watchlistL2, \
early_stopping_rounds=5000, verbose_eval=500) 
predsL2 = modelL2.predict(test_X)
rmsleL2 = rmsle(predsL2, test_y)
print ("LightGBM2 RMSLE = " + str(rmsleL2))
Training until validation scores don't improve for 5000 rounds.
[500]   training's rmse: 0.3948 valid_1's rmse: 0.660834
[1000]  training's rmse: 0.32773    valid_1's rmse: 0.677107
[1500]  training's rmse: 0.285268   valid_1's rmse: 0.687161
[2000]  training's rmse: 0.253586   valid_1's rmse: 0.693729
[2500]  training's rmse: 0.23149    valid_1's rmse: 0.702132
[3000]  training's rmse: 0.215441   valid_1's rmse: 0.707667
[3500]  training's rmse: 0.201827   valid_1's rmse: 0.710812
[4000]  training's rmse: 0.188848   valid_1's rmse: 0.713719
[4500]  training's rmse: 0.177358   valid_1's rmse: 0.715694
[5000]  training's rmse: 0.168306   valid_1's rmse: 0.71957
Early stopping, best iteration is:
[26]    training's rmse: 0.578606   valid_1's rmse: 0.62066
LightGBM2 RMSLE = 0.159046033321

各学習器単体での結果は、Ridge回帰: 0.1504、LightGBM1回目: 0.1524、LightGBM2回目: 0.1587でした。 これを組み合わせることで以下の様な結果が得られます。

preds = predsR*0.3 + predsL1*0.35 + predsL2*0.35
rmsle = rmsle(preds, test_y)
print ("Total RMSLE = " + str(rmsle))
Total RMSLE = 0.143722891253

掛け合わせることで、0.1435までスコアを上げることができました!ここからさらにスコアを上げていくには、各学習器のパラメータや学習器同士の割合をいろいろ調整していくことが必要になります。 では最後に実際の価格と予測価格の分布を見て終わります。

actual_price = np.expm1(test_y)
preds_price = np.expm1(preds)

plt.figure(figsize=(12,10))
cm = plt.cm.get_cmap('winter')
x_diff = np.clip(100 * ((preds_price - actual_price) / actual_price), -75, 75)
plt.scatter(x=actual_price, y=preds_price, c=x_diff, s=10, cmap=cm)
plt.colorbar()
plt.plot([0, 100], [0, 100], 'k--', lw=1)
plt.xlim(0, 100)
plt.ylim(0, 100)
plt.title('Actual vs. Predicted Prices')
plt.xlabel('Actual Prices [$]')
plt.ylabel('Predicted Prices [$]')
plt.show()

なんか今回はマニアックであまり面白くないものになりましたが、個人的には勉強になりました!以上です。