Deep learningとXGBoost - 「レンタル自転車の利用者予測」を題材に

# 目的
SIGNATEさんの練習問題として公開されている[「レンタル自転車の利用者予測」](https://signate.jp/competitions/114) に挑戦します。レンタル自転車の利用者数を予測します。[「重回帰分析」](https://qiita.com/Ambition_TY/private/34ee1b0eb932b6ec82f0)を試しましたが、あまり良い精度のモデルを得られませんでした。Deep leariningをやって、良い精度のモデルを作りたいです。

# 目的変数と説明変数を選択する
利用者数を予測するコンペですので、説明変数は今回も利用者数('cnt')です。

時間('hr')を説明変数で使うのをやめて、曜日('weekday')を代わりに使います。土日と平日で明確な特徴がありました。

![output_graph_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/688182/d658ed17-6df4-3f03-2209-9c18ff20ca8b.png)

その他の説明変数は、利用者数('cnt')と相関や関係がありそうな変数である温度('temp'),体感温度('atemp')、そして回帰分析でも使用した天気('weathersit')も追加します。天気('weathersit')はカテゴリーデータなのでダミー変数にして投入します。

目的変数・・・利用者('cnt')
説明変数・・・曜日('weekday')、温度('temp')、体感温度('atemp')、天気('weathersit')
でDeep learningを実施します。

# 必要なモジュールをインポートする

```python
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import tensorflow as tf
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Activation, Dense, Dropout, Input, BatchNormalization
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
```
tensorflowを使います。

# Deep learniningのコード
元データの読み込み込みは
```python
import pandas as pd
df = pd.read_csv("ファイルのパスをコピペします/train.csv")
```
で実行してください。左辺のdfは任意ですので、お好みの名称があればどうぞ。

# Deep learningの実装
いずれにしても目的変数と説明変数は必要です。説明変数は曜日、気温、体感温度、天気を使いますので、必要な目的変数(df1_y)と説明変数(df1_x)を用意します。必要な説明変数だけ取り出すために、locを使いました。

```python
df1_x = df.loc[:,['weekday','temp','atemp','weathersit']]
df1_y = df['cnt']
df1_x.head()
```
曜日の'weekday'を土日(0)と平日(1)の2値データにします。
```python
#新しい変数'weekday2'を作って、全て0(zero)を入れる
df1_x['weekday2'] = 0

# 平日を1にしたいので、1にするコードを用意する
df1_x.loc[df1_x['weekday'].isin(list(range(1,6))),'weekday2']=1

# 結果の件数を確認する
df1_x['weekday2'].value_counts()
```
土日が2495件、平日が6150件です。そしてこの2値にした曜日の変数('weekday2')と天気を示す変数('weathersit')をダミー変数に変換します。

```python
df1_x = pd.get_dummies(df1_x,columns=['weekday2','weathersit'])
```
元のweekdayは使わないのでdf1_xから落とします。念のため、df2_xで実行します。
```python
df2_x = df1_x.drop('weekday',axis=1)
```
df2_xを作ったので、yも合わせます。
```python
df2_y = df1_y
```
下準備を終えたので、実装に向かっていきます。

```python
# defaultはtrain dataが75%、test dataが25%に分割されます。
# また時系列のデータなので、shuffleをFalseにします。
train_X, test_X, train_y, test_y = train_test_split(df2_x, df2_y, random_state=42, shuffle=False)

# input_dimに説明変数の数を入れます。
# Denseで層を設定します。
model = Sequential()
model.add(Dense(256, input_dim=8, activation='relu'))
model.add(Dense(128, activation='relu'))
model.add(Dropout(0.5)) #128のNordを50%落とす
model.add(Dense(1))

# 損失関数にmse、最適化関数にadamを採用
model.compile(loss='mse', optimizer='adam')

# モデルを学習させます
# epochsは学習回数 データ数(6483) x 学習回数
# batch_size データ数(6483)を細かく分ける
train_model = model.fit(train_X, train_y, epochs= 60, batch_size= 24, verbose=1, validation_data=(test_X, test_y))
```
epochsが学習回数です。batch_sizeを24としたのは、データが1時間ごとのデータなので、24で1日分ということで、24を入れました。1週間であれば、168を入れます。

```python
# 予測値を出力します
y_pred = model.predict(test_X)

# 二乗誤差を出力します
mse= mean_squared_error(test_y, y_pred)
print("REG RMSE : %.2f" % (mse** 0.5))

# epoch毎の予測値の正解データとの誤差を表しています
# バリデーションデータのみ誤差が大きい場合、過学習を起こしています

loss=train_model.history['loss']
val_loss=train_model.history['val_loss']
epochs=len(loss)

plt.plot(range(epochs), loss, marker = '.', label = 'loss')
plt.plot(range(epochs), val_loss, marker = '.', label = 'val_loss')
plt.legend(loc = 'best')
plt.grid()
plt.xlabel('epoch')
plt.ylabel('loss')
plt.show()
```
学習結果です。RMSEが122.87でした
![output_deeplearining.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/688182/cbf0fd1b-11d9-cb1b-77da-4a3f81ef3693.png)

lossが14000あたりで停滞していますし、val_lossが飛んだり跳ねたりしています。あまり良いモデルではないようです。

test.csvもtrain dataと同じフォーマットにします。やる事は同じです。

```python
# Loading test data (test.csv)
import pandas as pd
df_test = pd.read_csv("ファイルのパスをコピペします/test.csv")
#df_test.head()

# df_xで使いたい説明変数を抜き出す (locの使い方に慣れましょう)
# 説明変数(df_x)と目的変数(df_y)を作る
df1_test_x = df_test.loc[:,['weekday','temp','atemp','weathersit']]
#df1_test_y = df['cnt'] #test.csvには'cnt'はありません。それを予想するのが今回の目的
#df1_test_x.head()

# 'weekday'を休日=0、平日=1の2値データにする
# isin以外での作成方法は「カテゴリー変換_6値を2値にする.ipynb」を参照

df1_test_x['weekday2'] = 0 #新しい変数'weekday2'を作って、全て0(zero)を入れる

# 平日を1にしたいので、1にするコードを用意する
df1_test_x.loc[df1_test_x['weekday'].isin(list(range(1,6))),'weekday2']=1
df1_test_x['weekday2'].value_counts()

# weathersitのカテゴリーデータをダミー変数(One-Hot Encoding)に置き換える
df1_test_x = pd.get_dummies(df1_test_x,columns=['weekday2','weathersit'])
df1_test_x

# 'weekday'が不要
df2_test_x = df1_test_x.loc[:,['temp','atemp','weekday2_0','weekday2_1','weathersit_1','weathersit_2','weathersit_3','weathersit_4']]
df2_test_x.count()

# test dataにモデルを適用する
X_test = df2_test_x[['atemp','weekday2_0','weekday2_1','weathersit_1','weathersit_2','weathersit_3','weathersit_4']] 
Y_test_pred = model.predict(X_test) #X_testデータを使って予測する
```
結果をcsvに保存します。
```python
sample = pd.read_csv('ファイルのパスをコピペします/sample_submit.csv', header=None)
sample[1] = Y_test_pred
sample.to_csv('ファイルのパスをコピペします/sample_submit.csv',index=None, header=None)
```
思ったほど精度が上がりません。そこでXGBoostを使います。

# XGBoost
決定木モデルのXGBoostでやってみます。XGBoostといつも使うモジュールをインポートします。

```python
# XGBoostをインストール
import xgboost as xgb
from sklearn.model_selection import GridSearchCV
from sklearn.metrics import mean_squared_error

# 使いそうなモジュールを読み込む
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import tensorflow as tf
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Activation, Dense, Dropout, Input, BatchNormalization
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
```

```python
train dataとtest dataを読み込みます。
df_train = pd.read_csv("ファイルのパスをコピペします/train.csv")
df_test = pd.read_csv("ファイルのパスをコピペします/test.csv")

# 説明変数、目的変数を格納
# 説明変数で使用しない'id','dteday','cnt'を落とします。
X_train = df_train.drop(['id','dteday','cnt'], axis=1)
# 目的変数は利用者数('cnt')なので、それだけ取り出します。
y_train = df_train['cnt']
X_test = df_test.drop(['id','dteday'], axis=1)
```
下準備が終わったので、XGBoostを実装します。

```python
# ハイパーパラメータ探索
# n_estimatorsで決定木の数を決めます。
reg_cv = GridSearchCV(reg, {'max_depth': list(range(2,10)), 'n_estimators': [50,100,200]}, verbose=1)
reg_cv.fit(X_train, y_train)

# 改めて最適パラメータで学習
reg = xgb.XGBRegressor(**reg_cv.best_params_)
reg.fit(X_train, y_train)

# 学習データを使って予測
pred_train = reg.predict(X_train)
```

出力結果は次の通りです。
![output_XGBoost_3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/688182/6e442acc-d32b-21e5-06da-fa04f4983334.png)
決定木の深さを示す'max_depth'が5、決定木の数を示す'n_estimators'が100、決定係数が0.823でした。それなりのモデルになっていると言って良いと思います。

二乗誤差を出力します
```python
mse= mean_squared_error(y_train, pred_train)
print("REG RMSE : %.2f" % (mse** 0.5))
```
RMSEは32.31まで向上しました。



予測値がtrain dataの利用者数とどのくらい一致しているでしょうか。
![output_XGBoost_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/688182/1f3deeff-3249-bc3c-26a1-3257590be06c.png)
良い感じで一致しています。しかし、ところどころ、予測値('pred')がマイナスになっていることが気になります。

重要な変数を見ることもできます。
![output_XGBoost_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/688182/91495994-4f45-f6c1-9923-5458ced0a348.png)
時間や曜日、気温が影響を与えているようです。

テストデータを使って予測します。
```python
# test.csvを読み込みます
df_test = pd.read_csv("ファイルのパスをコピペします/test.csv")

# 説明変数を格納（目的変数はtest dataで与えらえていないので不要です）
X_test = df_test.drop(['id','dteday'], axis=1)
```

モデルに当てはめます。
```python
pred_test = reg.predict(X_test)
```

先ほどと同じように結果をcsvファイルに変換して、SIGNATEのサイトに提出します。

結果は、回帰分析やDeep leariningよりも良くなり、暫定評価127まであがりました。

# まとめと考察
XGBoostを使う事で精度がよくなりました。Deep leariningやXGBoostのコードを書いてお仕舞いではなく、分析作業はそこからはじまるのであって、説明変数の吟味や層の数などいろいろパラメーターを変えて実行することを繰り返すことが必要だと思いました。