はじめに
kaggle 超初心者のtomoです。
巷でKaggle上位ランカーはXgboostという学習モデルを利用している噂を耳に…
「自分もその波に乗っかりたい」
そんな安易な気持ちで初めました。
Xgboostとは
Xgboost=決定木×ブースティング
1:決定木(分類木&回帰木)
決定木単体では精度がそれほど高くない…
そこでアンサンブル学習のひとつであるブースティングを組み合わせると最強になるらしい。
2:ブースティング
- アンサンブル学習…簡単に言えば多数決をとる方法で、個々に別々の学習器として学習させたものを、融合させること。 ex)バギング,ブースティング,スタッキング
- ブースティング…アンサンブル学習のひとつ。学習データ&予測モデルを逐次的に生成&構築していく方法。
3:バギング(左)とブースティング(右)の違い
Xgboostの仲間紹介
- LightGBM…Xgboostの2年後に開発されたXgboostの改良モデル。計算負荷が軽く実装できる。
Xgboostの決定木はlevel（層）が成長していく。
LightGBMの決定木はleaf（葉）が成長していく。
CatBoost…データセットによってXgboostやLightGBMよりも精度が高くなる可能性を秘めている学習手法(特にカテゴリアルデータのとき有用)。
実装
今回は過去のkaggleコンペの
New York City Taxi Fare PredictionをXgboostを用いてやってみます。
やってみたい方はこちら
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()
考察
- key…訓練データにおいて不要なので、無視します。
- pickup-datetime…客が乗車した時刻
- pickup-longtitude…客を乗せた際の経度
- pickup_latitude…客を乗せた際の緯度
- dropoff_longitude…客が降車した際の経度
- dropoff_latitude…客が降車した際の緯度
- passenger_count…乗車客数
欠損値,型
#欠損値,型を確認
df.info()
考察
欠損値は含まれていないことが確認できる。(実際オリジナルデータの1/1000なので偶然この中に欠損値が含まれませんでした。)
データのクリーニング
各種統計量
df.describe()
考察
- 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()
#四分位数を使った方法(四分位範囲の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()
(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()
#乗客=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()
考察
上図から新たにカラムが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回連続で性能が上昇しない点で止まるか
グリッドリサーチでハイパーパラメータを最適化
#しらみ潰しに全部試す
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
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))
#仮に上の最適ハイパーパラメータは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()
#xgboostの便利な機能
#それぞれの特徴量が予測にどれだけの影響を与えているかを確認
fig, ax = plt.subplots(figsize=(12, 12))
xgb.plot_importance(model, max_num_features=12, height=0.8, ax=ax)
plt.show()
補足(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()
まとめ
Xgboostはモデル学習にかなりの計算時間を費やすため、あまりにこ大きなデータセットには向かないことが実験して分かりました。(今回の5500件なら問題なし)
上記の図に書かれている"影響度を表すグラフ"から特徴量エンジニアリングで作り出したものが上位にあることが分かりますね。
これだけ特徴量を決めるのが大事であることもXgboostを使って再認識させられました。
次回はLightGBMを使ってどれだけRMSEの値が減少するのか増加するのかを確認していきます。
