Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
Help us understand the problem. What is going on with this article?

最強教師あり学習手法??Xgboostの正体を探る

はじめに

kaggle 超初心者のtomoです。
巷でKaggle上位ランカーXgboostという学習モデルを利用している噂を耳に…
「自分もその波に乗っかりたい」
そんな安易な気持ちで初めました。

Xgboostとは

Xgboost=決定木×ブースティング

1:決定木(分類木&回帰木)

分類木 回帰木

決定木単体では精度がそれほど高くない…

そこでアンサンブル学習のひとつであるブースティング組み合わせる最強になるらしい。

2:ブースティング

  • アンサンブル学習…簡単に言えば多数決をとる方法で、個々に別々の学習器として学習させたものを、融合させること。 ex)バギング,ブースティング,スタッキング
  • ブースティングアンサンブル学習のひとつ。学習データ&予測モデル逐次的に生成&構築していく方法。

3:バギング(左)とブースティング(右)の違い

バギング ブースティング

Xgboostの仲間紹介

  • LightGBM…Xgboostの2年後に開発されたXgboostの改良モデル計算負荷が軽く実装できる。

Xgboostの決定木level(層)が成長していく。
Xgboost

LightGBMの決定木leaf(葉)が成長していく。
LightGBM

CatBoost…データセットによってXgboostやLightGBMよりも精度が高くなる可能性を秘めている学習手法(特にカテゴリアルデータのとき有用)。

実装

今回は過去のkaggleコンペの
New York City Taxi Fare PredictionXgboostを用いてやってみます。
やってみたい方はこちら

image.png

https://www.kaggle.com/c/new-york-city-taxi-fare-prediction

データセット

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline 
plt.style.use('ggplot')
import seaborn as sns

#実際のデータは55000000個ある。でかすぎる…
#今回はその1/10000
df = pd.read_csv('train.csv', nrows=5500)
df.head()

image.png

考察

  • key…訓練データにおいて不要なので、無視します。
  • pickup-datetime…客が乗車した時刻
  • pickup-longtitude…客を乗せた際の経度
  • pickup_latitude…客を乗せた際の緯度
  • dropoff_longitude…客が降車した際の経度
  • dropoff_latitude…客が降車した際の緯度
  • passenger_count…乗車客数

欠損値,型

#欠損値,型を確認
df.info()

LightGBM

考察

欠損値は含まれていないことが確認できる。(実際オリジナルデータの1/1000なので偶然この中に欠損値が含まれませんでした。)

データのクリーニング

各種統計量

df.describe()

image.png

考察

  • fareのmin=-2.9になっているので処理を加える必要あり。
  • ニューヨークの緯度・経度は大体(-74°・41°)ということをおさえておく。(wikiで調べました。)
  • passsengerのmin=0になっているので処理を加える必要あり。
  • あくまで1/1000個しかないので注意。

(fare)運賃の前処理

#料金が正の値のみを抽出
df = df[df['fare_amount'] > 0]

#分布の形を確認
plt.figure(figsize=(12, 8))
plt.hist(df['fare_amount'], bins=100)
plt.show()

image.png

#四分位数を使った方法(四分位範囲の1.5倍以上はなれているものを除外)
#外れ値の処理の続き
def outlier_iqr(df, columns=None):
    if columns == None:
        columns = df.columns

    for col in columns:
        q1 = df[col].describe()['25%']
        q3 = df[col].describe()['75%']
        #四分位範囲
        iqr = q3 - q1 

        outlier_min = q1 - iqr * 1.5
        outlier_max = q3 + iqr * 1.5

        # 範囲から外れている値を除く
        df = df[(df[col] >= outlier_min) & (df[col] <= outlier_max)]

    return df

df = outlier_iqr(df, columns=['fare_amount'])
plt.figure(figsize=(12, 8))
plt.hist(df['fare_amount'], bins=50)
plt.show()

image.png

(pickup_longitude&pickup_latitude&dropoff_longitude&dropoff_latitude)緯度・経度の前処理

# 経度・緯度 +-1°で絞る
print('before', df.shape)
df = df[(df['pickup_longitude'] > -75) & (df['pickup_longitude'] < -73)]
df = df[(df['pickup_latitude'] > 40) & (df['pickup_latitude'] < 42)]
df = df[(df['dropoff_longitude'] > -75) & (df['dropoff_longitude'] < -73)]
df = df[(df['dropoff_latitude'] > 40) & (df['dropoff_latitude'] < 42)]

乗客数の前処理

df['passenger_count'].value_counts()

image.png

#乗客=0を取り除けばいいことがわかる
df = df[(df['passenger_count'] > 0)]

特徴量エンジニアリング

今回は新たに7つの特徴量を追加します。
- 「distance」…走行距離
- 「azimuth」…走行した方角
- 「hour」,「day」,「month」「dayofweek」「year」
* Series.dt.dayofweekを使うと、月曜日= 0、日曜日= 6にして返してくれます。

#距離の計算
def calculation_distance(x_1, y_1, x_2, y_2):
    # 赤道半径 (km)
    R = 6371
    # Radian角に変換
    _x1, _y1, _x2, _y2  = map(np.radians, [x_1, y_1, x_2, y_2])

    delta_x = _x2 - _x1
    delta_y = _y2 - _y1

    # 距離を計算
    a = np.sin(delta_y/2.0)**2 + np.cos(_y1) * np.cos(_y2) * np.sin(delta_x/2.0)**2
    return 2 * R * np.arcsin(np.sqrt(a))
#方角の計算
def calculation_azimuth(x_1, y_1, x_2, y_2):
    # Radian角に修正
    _x1, _y1, _x2, _y2 = map(np.radians, [x_1, y_1, x_2, y_2])

    delta_x = _x2 - _x1
    _y = np.sin(delta_x)
    _x = np.cos(_y1) * np.tan(_y2) - np.sin(_y1) * np.cos(delta_x)

    psi = np.rad2deg(np.arctan2(_y, _x))

    return np.where(psi < 0, 360 + psi, psi)
#距離のカラムを作成
df['distance'] = calculation_distance(df['pickup_longitude'],
                                      df['pickup_latitude'],
                                      df['dropoff_longitude'],
                                      df['dropoff_latitude'])
#方角のカラムを作成
df['azimuth'] = calculation_azimuth(df['pickup_longitude'],
                                    df['pickup_latitude'],
                                    df['dropoff_longitude'],
                                    df['dropoff_latitude'])
#日時のデータから特徴量を取得
df['pickup_datetime'] = pd.to_datetime(df['pickup_datetime'], format='%Y-%m-%d %H:%M:%S UTC')

df['hour'] = df['pickup_datetime'].dt.hour
df['day'] = df['pickup_datetime'].dt.day
df['month'] = df['pickup_datetime'].dt.month
df['dayofweek'] = df['pickup_datetime'].dt.dayofweek
df['year'] = df['pickup_datetime'].dt.year

df.head()

image.png

考察

上図から新たにカラムが7つ作られていることが確認できる。

モデル構築

from sklearn.model_selection import train_test_split
import xgboost as xgb

#不要カラム削除
X = df.drop(columns=['key', 'fare_amount', 'pickup_datetime']) #説明変数
y = df['fare_amount'] #目的変数

#訓練データを訓練データとテストデータに分割
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2,
                                                    random_state=0)
#xgboostはパラメータを決める必要がある…
params = {
    'silent': 1, 
    'max_depth': 6, #木の深さに上限
    'min_child_weight': 1, #葉の重みの下限
    'eta': 0.1, #学習率(どのくらいのスピードで学習するか)
    'tree_method': 'exact',
    'objective': 'reg:linear',#回帰分析だからreg:linear
    'eval_metric': 'rmse', #二乗誤差を評価の基準に
    'predictor': 'cpu_predictor' 
}

#xgboostをするには特殊な行列が必要…
dtrain = xgb.DMatrix(X_train, label=y_train)
dtest = xgb.DMatrix(X_test, label=y_test)
model = xgb.train(params=params,
                  dtrain=dtrain,
                  num_boost_round=1000, #1000個の決定木を重ねる
                  early_stopping_rounds=5, #5回連続で性能が上昇しなければそこで終了
                  evals=[(dtest, 'test')])

#つまり,1000回調整をするか、5回連続で性能が上昇しない点で止まるか

image.png

グリッドリサーチでハイパーパラメータを最適化

#しらみ潰しに全部試す
gridsearch_params = [
    (max_depth, eta)
    for max_depth in [6, 7, 8] #データによってrandom
    for eta in [0.1, 0.05, 0.01] #データによってrandom
]
gridsearch_params

image.png

cv_results = xgb.cv(
        params,
        dtrain,
        num_boost_round=10,
        seed=0,
        nfold=5,
        metrics={'rmse'},
        early_stopping_rounds=5
    )
cv_results
#めっちゃ時間がかかるので、最初のデータを少なくした方がいいかも
#どのハイパーパラメータがいいのかを知る
min_rmse = float('Inf') #初期値を無限大にする(atcorderっぽい)
best_param = []

for max_depth, eta in gridsearch_params:
    print('max_depth={}, eta={}'.format(max_depth, eta))

    params['max_depth'] = max_depth
    params['eta'] = eta

    #交差検証
    cv_results = xgb.cv(
        params,
        dtrain,
        num_boost_round=1000,
        seed=0,
        nfold=5, #5回繰り返す
        metrics={'rmse'},#今回のコンペの評価基準がrmseであるためこのように設定
        early_stopping_rounds=5
    )

    mean_rmse = cv_results['test-rmse-mean'].min()
    boost_rounds = cv_results['test-rmse-mean'].argmin() #何回目の学習で最終的な値を取り出せたか
    print('RMSE {} for {} rounds'.format(mean_rmse, boost_rounds))
    if mean_rmse < min_rmse:
        min_rmse = mean_rmse
        best_param = (max_depth, eta)

print('Best params {}, RMSE {}'.format(best_param, min_rmse))

image.png

#仮に上の最適ハイパーパラメータは6,0.05なら
params['max_depth'] = 6
params['eta'] = 0.05

model = xgb.train(params=params,
                  dtrain=dtrain,
                  num_boost_round=1000,
                  early_stopping_rounds=5,
                  evals=[(dtest, 'test')])
#横軸=正解値,縦軸=予測値
#傾き=1の直線が理想
prediction = model.predict(xgb.DMatrix(X_test), 
                           ntree_limit=model.best_ntree_limit)

plt.figure(figsize=(12, 12))
plt.scatter(y_test[:1000], prediction[:1000], alpha=0.2)
plt.show()

image.png

#xgboostの便利な機能
#それぞれの特徴量が予測にどれだけの影響を与えているかを確認
fig, ax = plt.subplots(figsize=(12, 12))
xgb.plot_importance(model, max_num_features=12, height=0.8, ax=ax)
plt.show()

image.png

補足(kaggle タイタニックでもざっくりやってみた)

import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import LabelEncoder
from sklearn.model_selection import train_test_split
import pandas as pd

train_df=pd.read_csv("train.csv")
test_df=pd.read_csv("test.csv")

all_df=pd.concat((train_df.loc[:,"Pclass":"Embarked"],test_df.loc[:,"Pclass":"Embarked"]))

all_df["Age"]=all_df["Age"].fillna(all_df["Age"].mean())
all_df["Fare"]=all_df["Fare"].fillna(all_df["Fare"].mean())
all_df["Embarked"]=all_df["Embarked"].fillna(all_df["Embarked"].mode()[0])

cat_features=["Sex","Embarked"]

for col in cat_features:
    lbl=LabelEncoder()
    all_df[col]=lbl.fit_transform(list(all_df[col].values))

all_df=all_df.drop(columns=["Name","Ticket","Cabin"])

train=all_df[:train_df.shape[0]]
test=all_df[train_df.shape[0]:]

y=train_df["Survived"]
ID=test_df["PassengerId"]

X_train,X_test,y_train,y_test=train_test_split(train,y,random_state=0)

import xgboost as xgb

params = {
    'silent': 1, 
    'max_depth': 6, #木の深さに上限
    'min_child_weight': 1, #葉の重みの下限
    'eta': 0.1, #学習率(どのくらいのスピードで学習するか)
    'tree_method': 'exact',
    'objective': "binary:logistic",#二値分類だからbinary:logistic
    'eval_metric': 'auc', #AUCを評価の基準に
    'predictor': 'cpu_predictor' 
}

#xgboostをするには特殊な行列が必要…
dtrain = xgb.DMatrix(X_train, label=y_train)
dtest = xgb.DMatrix(X_test, label=y_test)

model = xgb.train(params=params,
                  dtrain=dtrain,
                  num_boost_round=100, #100個の決定木を重ねる
                  early_stopping_rounds=10, #10回連続で性能が上昇しなければそこで終了
                  evals=[(dtest, 'test')]
                 )

prediction = model.predict(xgb.DMatrix(test), 
                           ntree_limit=model.best_ntree_limit)
prediction=np.where(prediction<0.5,0,1)

submisson=pd.DataFrame({
    "PassengerId":ID,
    "Survived":prediction
})
submisson.to_csv("submisson2.csv",index=False)

fig, ax = plt.subplots(figsize=(12, 12))
xgb.plot_importance(model, max_num_features=12, height=0.8, ax=ax)
plt.show()

image.png

まとめ

Xgboostはモデル学習にかなりの計算時間を費やすため、あまりにこ大きなデータセットには向かないことが実験して分かりました。(今回の5500件なら問題なし)

上記の図に書かれている"影響度を表すグラフ"から特徴量エンジニアリングで作り出したものが上位にあることが分かりますね。
これだけ特徴量を決めるのが大事であることもXgboostを使って再認識させられました。
次回はLightGBMを使ってどれだけRMSEの値が減少するのか増加するのかを確認していきます。

tomoxxx
こまったことを記事にするよ
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away