はじめに
今回SIGNATEのSOTA Challenge「オペレーション最適化に向けたシェアサイクルの利用予測」に参加して銅メダルを獲得することができたため、コンペの概要と自身の解法について書いていこうと思います。
コードについてはGitHubに公開しています。
コンペの概要
シェアサイクル事業を展開するA社がオペレーション改善のために、どのステーションで何台くらい自転車に空きがあるかを機械学習により予測したい。
⇒与えられたデータを用いて各々のステーションにおける空き台数を予測する
評価指標としてはRMSEが用いられている。RMSEの算出方法は以下の通り。
$$ RMSE = \sqrt{\frac{1}{n} \sum_{i=1}^n(y_i-\hat{y_i})^2} $$
データ概要
与えられたデータは以下の4つ。
1. 自転車の台数状況データ
各サイクルステーションで1時間ごとに記録された利用可能な自転車数（目的変数）の履歴データ
| カラム | データ型 | 説明 |
|---|---|---|
| id | int | インデックスとして使用 |
| year | int | 記録日時(年) |
| month | int | 記録日時(月) |
| day | int | 記録日時(日) |
| hour | int | 記録日時(時) |
| station_id | int | ステーションid |
| bikes_available | int | 利用可能な自転車数（台数） |
| predict | int | 予測対象フラグ（0=予測対象外, 1=予測対象 |
2. 利用者の移動履歴データ
利用者がシェアサイクルで移動した時間、起点駅、終点駅等を記録した移動履歴データ
| カラム | データ型 | 説明 |
|---|---|---|
| trip_id | int | 移動履歴ID |
| duration | int | 移動時間（秒) |
| start_date | char | 移動を開始した日時 |
| start_station_id | int | 移動を開始したステーションID |
| end_date | char | 移動を終了した日時 |
| end_station_id | int | 移動を終了したステーションID |
| bike_id | int | 利用された自転車ID |
| subscription_type | char | 利用者の登録種別（Customer=一時利用者, Subscriber=定額プラン利用者） |
3. ステーション情報
サイクルステーションの緯度・経度、ドック数（最大で停められる自転車数）、設置日のデータ
| カラム | データ型 | 説明 |
|---|---|---|
| station_id | int | ステーションID |
| lat | float | ステーション設置場所の緯度 |
| long | float | ステーション設置場所の経度 |
| dock_count | int | ドック数（最大で停められる自転車数） |
| city | char | ステーションが設置されている都市（具体的な都市名は非開示） |
| installation_date | char | ステーションが設置された日時 |
4. 気象情報
都市中心部における1日ごとの気象予報データ(0時時点)
各カラムの説明は割愛
方針
それぞれのデータセットに関してEDAを行いデータの理解を深め、LightGBMによるモデル構築を行う。
- 2013/9/1~2014/8/31のデータはすべて予測対象ではないため、学習データとしてモデル構築に用いる。
- 2014/9/1~2015/8/31のデータのうち予測対象でないものを検証用データとして用いる。
- 利用者の移動情報は時間帯ごと、都市ごとに集計して特徴量とする。
EDA
実施したデータ探索のうちいくつかを抜粋して記載する。
自転車の台数状況データ
時間帯別の空き台数をみると、夜間や早朝に多く、日中に少ない傾向にあることがわかる。
(ただしステーションによっては日中に多いところもある。)
自転車の台数状況データ
status = pd.read_csv('status.csv')
# 日付をdatetime型で作成
status['date'] = status['year'].astype(str) + status['month'].astype(str).str.zfill(2) + status['day'].astype(str).str.zfill(2)
status['date'] = pd.to_datetime(status['date'])
# 時間帯別に空き台数を確認
tmp_df = status.loc[status['date'] < '2014-09-01']
bikes_data_hour = tmp_df.groupby('hour')['bikes_available'].mean()
plt.figure(figsize=(10, 6))
plt.plot(bikes_data_hour.index, bikes_data_hour.values)
plt.title('時間帯別空き台数')
plt.show()
利用者の移動履歴データ
都市はcity1~city5までの5種類、ドック数は27,15,11,19,25,23の6パターンある。
ステーション情報
station = pd.read_csv('station.csv', parse_dates=['installation_date'])
station.groupby('city')['dock_count'].value_counts()
city dock_count
city1 15 10
19 4
11 1
27 1
city2 19 14
15 12
23 6
27 3
city3 15 6
25 1
city4 15 4
23 2
11 1
city5 11 2
15 2
23 1
Name: dock_count, dtype: int64
特徴量生成
すべてのデータセットを年月日時およびステーションIDをキーとして結合して新たにデータセットを作成し、predictが1のものを学習データ(train_df), 0のものをテストデータ(test_df)とする。
- 天候情報について
雨や霧の日を1とするダミー変数を作成する。 - 都市名について
Ligth GBMを使うため、ラベルエンコーディングを適用する。 - 時間帯別の返却数や貸出数について
集計結果を各ステーションのドック数で割る。(こうしたほうが目的変数との相関が高くなった) - 時刻について
9時から16時の間は1をとるダミー変数を作成。
特徴量生成
from sklearn.preprocessing import LabelEncoder
# 天候情報
train_df['events'] = train_df['events'].apply(lambda x: 1 if x is None else 0)
evaluate_df['events'] = evaluate_df['events'].apply(lambda x: 1 if x is None else 0)
test_df['events'] = test_df['events'].apply(lambda x: 1 if x is None else 0)
# 都市名をラベルエンコーディング
le = LabelEncoder()
encoded = le.fit_transform(train_df['city'].values)
train_df['city'] = encoded
encoded = le.fit_transform(evaluate_df['city'].values)
evaluate_df['city'] = encoded
encoded = le.fit_transform(test_df['city'].values)
test_df['city'] = encoded
# 返却数等をドック数で割る
train_df['rent_per_dock_1'] = train_df['rent_count_1'] / train_df['dock_count']
train_df['return_per_dock_1'] = train_df['return_count_1'] / train_df['dock_count']
train_df['start_subscriber_per_dock_1'] = train_df['start_subscriber_1'] / train_df['dock_count']
train_df['end_subscriber_per_dock_1'] = train_df['end_subscriber_1'] / train_df['dock_count']
evaluate_df['rent_per_dock_1'] = evaluate_df['rent_count_1'] / evaluate_df['dock_count']
evaluate_df['return_per_dock_1'] = evaluate_df['return_count_1'] / evaluate_df['dock_count']
evaluate_df['start_subscriber_per_dock_1'] = evaluate_df['start_subscriber_1'] / evaluate_df['dock_count']
evaluate_df['end_subscriber_per_dock_1'] = evaluate_df['end_subscriber_1'] / evaluate_df['dock_count']
test_df['rent_per_dock_1'] = test_df['rent_count_1'] / test_df['dock_count']
test_df['return_per_dock_1'] = test_df['return_count_1'] / test_df['dock_count']
test_df['start_subscriber_per_dock_1'] = test_df['start_subscriber_1'] / test_df['dock_count']
test_df['end_subscriber_per_dock_1'] = test_df['end_subscriber_1'] / test_df['dock_count']
train_df['rent_per_dock_7'] = train_df['rent_count_7'] / train_df['dock_count']
train_df['return_per_dock_7'] = train_df['return_count_7'] / train_df['dock_count']
train_df['start_subscriber_per_dock_7'] = train_df['start_subscriber_7'] / train_df['dock_count']
train_df['end_subscriber_per_dock_7'] = train_df['end_subscriber_7'] / train_df['dock_count']
evaluate_df['rent_per_dock_7'] = evaluate_df['rent_count_7'] / evaluate_df['dock_count']
evaluate_df['return_per_dock_7'] = evaluate_df['return_count_7'] / evaluate_df['dock_count']
evaluate_df['start_subscriber_per_dock_7'] = evaluate_df['start_subscriber_7'] / evaluate_df['dock_count']
evaluate_df['end_subscriber_per_dock_7'] = evaluate_df['end_subscriber_7'] / evaluate_df['dock_count']
test_df['rent_per_dock_7'] = test_df['rent_count_7'] / test_df['dock_count']
test_df['return_per_dock_7'] = test_df['return_count_7'] / test_df['dock_count']
test_df['start_subscriber_per_dock_7'] = test_df['start_subscriber_7'] / test_df['dock_count']
test_df['end_subscriber_per_dock_7'] = test_df['end_subscriber_7'] / test_df['dock_count']
# 時刻ダミー
train_df['daytime'] = train_df['hour'].apply(lambda x: 1 if x>=9 and x<=16 else 0)
evaluate_df['daytime'] = evaluate_df['hour'].apply(lambda x: 1 if x>=9 and x<=16 else 0)
test_df['daytime'] = test_df['hour'].apply(lambda x: 1 if x>=9 and x<=16 else 0)
モデリング
Light GBMで構築
KFoldで交差検証を行う
各Foldで評価データRMSEは3.55前後
light gbm
import lightgbm as lgb
from sklearn.model_selection import KFold, train_test_split
from sklearn.metrics import mean_squared_error
rmse_list = []
model_list = []
pred_list = []
kf = KFold(n_splits=5, shuffle=True, random_state=0)
for index, (train_index, test_index) in enumerate(kf.split(X, y)):
print(f'{index+1}回目')
X_train = X.iloc[train_index]
X_valid = X.iloc[test_index]
y_train = y.iloc[train_index]
y_valid = y.iloc[test_index]
gbm = lgb.LGBMRegressor(objective='mse', bagging_fraction=0.5, colsample_bytree=1, n_estimators=10000, learning_rate=0.01,
max_depth=5, min_child_samples=200, num_leaves=255, random_state=0, reg_alpha=1.5, reg_lambda=2)
eval_set = [(X_valid, y_valid)]
callbacks = []
callbacks.append(lgb.early_stopping(stopping_rounds=10))
callbacks.append(lgb.log_evaluation())
gbm.fit(X_train, y_train, eval_set=eval_set, callbacks=callbacks)
y_train_pred = gbm.predict(X_train)
y_valid_pred = gbm.predict(X_valid)
y_evaluate_pred = gbm.predict(X_evaluate)
y_test_pred = gbm.predict(X_test)
train_mse = mean_squared_error(y_train_pred, y_train)
valid_mse = mean_squared_error(y_valid_pred, y_valid)
evaluate_mse = mean_squared_error(y_evaluate_pred, y_evaluate)
train_rmse = np.sqrt(train_mse)
valid_rmse = np.sqrt(valid_mse)
evaluate_rmse = np.sqrt(evaluate_mse)
print(f'学習データRMSE: {train_rmse}')
print(f'検証データRMSE: {valid_rmse}')
rmse_list.append(evaluate_rmse)
pred_list.append(y_test_pred)
評価データRMSE
[3.5536094045907913, 3.533682412715398, 3.541065851811445, 3.561579978453426, 3.5489761957624895]
結果
予測値は各Foldで算出した予測値の平均とした。
提出した予測値のスコアが3.639で順位が70位(2022/09/23現在)
感想・改善点
- あまり込み入った特徴量エンジニアリングをしなくてもLightGBMで比較的高いスコア(メダルギリギリ圏内程度)が出たため、スコアの改善に苦労した。
- 特徴量として休日(土日)ダミーを用いたが、祝日に関しては考慮していなかったので祝日情報も特徴量として含めてもよかった。
- LightGBMなどのツリー系のモデルしか用いていないため、RNNなどの時系列予測のための深層学習モデルを構築する余地もあり。
- 記事作成にあたり改めてEDAを実施したところ、設置時期が遅いステーションがいくつかあり、それらに関しては学習データが足りないような気がする。