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

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

#Xgboostとは
##Xgboost=決定木×ブースティング
###1:決定木(分類木&回帰木)
<img width="300" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/449715/143f3f1a-a737-2b66-3914-09d509d5097f.png" class="img" alt="分類木"> <img width="300" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/449715/bf01110b-f397-2af2-0f54-a7955a5ad438.png" class="img" alt="回帰木">

###決定木単体では精度がそれほど高くない…
そこで**アンサンブル学習のひとつであるブースティング**を**組み合わせる**と**最強**になるらしい。

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

###3:バギング(左)とブースティング(右)の違い
<img width="300" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/449715/ef41ebf2-d6d0-1258-6ab5-429b35f42937.png" class="img" alt="バギング">　<img width="240" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/449715/67a425fc-6615-3a0f-a729-f092eadb2457.png" class="img" alt="ブースティング">

#Xgboostの仲間紹介
- **LightGBM**…Xgboostの2年後に開発された**Xgboostの改良モデル**。**計算負荷が軽く**実装できる。

**Xgboostの決定木**は**level（層）**が成長していく。
<img width="300" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/449715/5d4705f0-d066-dc7c-26df-51f1212c10a7.png" class="img" alt="Xgboost">

**LightGBMの決定木**は**leaf（葉）**が成長していく。
<img width="300" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/449715/4a0cf689-7dd2-0002-3aa3-d3921ca81897.png" class="img" alt="LightGBM">

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

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

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/449715/bc085751-3123-bc46-88b4-c363384881f5.png)


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

#データセット
```py
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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/449715/2ff6dbce-e8dd-a6f9-eb9b-3604080e4de5.png)

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

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

<img width="500" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/449715/40fe4c81-19ed-ab99-fae4-db7dbfc58dac.png" class="img" alt="LightGBM">

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

#データのクリーニング
##各種統計量
```py
df.describe()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/449715/03595f59-332a-9e98-63a3-906a26bfe6f3.png)

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


##(fare)運賃の前処理
```py
#料金が正の値のみを抽出
df = df[df['fare_amount'] > 0]

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

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/449715/5ab0165c-634f-3552-3cfc-dfc7160f7bb0.png)


```py
#四分位数を使った方法(四分位範囲の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'])
```

```py
plt.figure(figsize=(12, 8))
plt.hist(df['fare_amount'], bins=50)
plt.show()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/449715/1db84712-264f-62a9-af12-f8025c1edb3e.png)


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

```py
# 経度・緯度 +-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)]
```

##乗客数の前処理
```py
df['passenger_count'].value_counts()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/449715/40498241-0719-c730-c1da-2b1519a9ccf5.png)

```py
#乗客=0を取り除けばいいことがわかる
df = df[(df['passenger_count'] > 0)]
```
#特徴量エンジニアリング
今回は新たに7つの特徴量を追加します。
- **「distance」**…走行距離
- **「azimuth」**…走行した方角
- **「hour」,「day」,「month」「dayofweek」「year」**
* **Series.dt.dayofweek**を使うと、**月曜日= 0、日曜日= 6**にして返してくれます。


```py
#距離の計算
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))
```

```py

#方角の計算
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)
```
```py
#距離のカラムを作成
df['distance'] = calculation_distance(df['pickup_longitude'],
                                      df['pickup_latitude'],
                                      df['dropoff_longitude'],
                                      df['dropoff_latitude'])
```

```py
#方角のカラムを作成
df['azimuth'] = calculation_azimuth(df['pickup_longitude'],
                                    df['pickup_latitude'],
                                    df['dropoff_longitude'],
                                    df['dropoff_latitude'])
```
```py
#日時のデータから特徴量を取得
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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/449715/db972992-f243-d24a-5a1a-109eed8c8425.png)

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

#モデル構築
```py
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)
```
```py
#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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/449715/356d29a3-ba07-c00b-cb06-2bbba900ce26.png)




#グリッドリサーチでハイパーパラメータを最適化
```py
#しらみ潰しに全部試す
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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/449715/636b59a2-0db5-c64e-e397-2ab8cc9016b1.png)

```py
cv_results = xgb.cv(
        params,
        dtrain,
        num_boost_round=10,
        seed=0,
        nfold=5,
        metrics={'rmse'},
        early_stopping_rounds=5
    )
cv_results
```

```py
#めっちゃ時間がかかるので、最初のデータを少なくした方がいいかも
#どのハイパーパラメータがいいのかを知る
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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/449715/6dc2c3d9-48ef-2f86-58e8-7eaae5a0b327.png)

```py
#仮に上の最適ハイパーパラメータは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')])
```

```py
#横軸=正解値,縦軸=予測値
#傾き=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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/449715/ce634143-72e7-b151-ceae-55c94cf53d74.png)

```py
#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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/449715/fbb87c09-14e6-b0e3-057d-5c21f59da690.png)

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

```py
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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/449715/0f3b6e84-40f0-2286-47dc-55728df738f2.png)





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

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