🛰️
[solafune] 夜間光データから土地価格を予測 BaseLine
はじめに
こんにちは。mst8823という名前でいろんなコンペに参加している者です。
この記事は、コンペプラットフォームsolafuneで行われている夜間光コンペのベースラインとなります。(CV: 0.5353, LB: 0.4918)
directory structure
ディレクトリ構成は以下を想定しています。
- input: TrainDataSet.csv などの配布データ置き場
- notebooks: ベースラインノートブック置き場
- submission: 提出用ファイル置き場
├─input
├─notebooks
└─submission
import
最初に必要なライブラリをインポートします
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from lightgbm import LGBMModel
from sklearn import model_selection
from sklearn.metrics import mean_squared_error
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from tqdm import tqdm
config
INPUT = "../input/"
warnings.filterwarnings('ignore')
TARGET = "AverageLandPrice"
load data
使うデータたちを読み込みます
train = pd.read_csv(INPUT+"TrainDataSet.csv")
test = pd.read_csv(INPUT+"EvaluationData.csv")
sample_sub = pd.read_csv(INPUT+"UploadFileTemplate.csv")
feature engineering
前処理 & 特徴量作成をします
はじめに特徴量作成で使用する汎用的な関数を定義しておきます
# aggregationのagg_methodsで使用
def max_min(x):
return x.max()-x.min()
def q75_q25(x):
return x.quantile(0.75) - x.quantile(0.25)
# xfeat.aggregationとほぼ同じ
def aggregation(input_df, group_key, group_values, agg_methods):
new_df = []
for agg_method in agg_methods:
for col in group_values:
if callable(agg_method):
agg_method_name = agg_method.__name__
else:
agg_method_name = agg_method
new_col = f"agg_{agg_method_name}_{col}_grpby_{group_key}"
df_agg = (input_df[[col] + [group_key]].groupby(group_key)[[col]].agg(agg_method))
df_agg.columns = [new_col]
new_df.append(df_agg)
_df = pd.concat(new_df, axis=1).reset_index()
output_df = pd.merge(input_df[[group_key]], _df, on=group_key, how="left")
return output_df.drop(group_key, axis=1)
# group 内で diffをとる関数
def diff_aggregation(input_df, group_key, group_values, num_diffs):
dfs = []
for nd in num_diffs:
_df = input_df.groupby(group_key)[group_values].diff(nd)
_df.columns = [f'diff={nd}_{col}_grpby_{group_key}' for col in group_values]
dfs.append(_df)
output_df = pd.concat(dfs, axis=1)
return output_df
# group 内で shiftをとる関数
def shift_aggregation(input_df, group_key, group_values, num_shifts):
dfs = []
for ns in num_shifts:
_df = input_df.groupby(group_key)[group_values].shift(ns)
_df.columns = [f'shift={ns}_{col}_grpby_{group_key}' for col in group_values]
dfs.append(_df)
output_df = pd.concat(dfs, axis=1)
return output_df
つぎに実際に特徴量を作成する関数たちを定義します
# そのままの値の特徴量
def get_raw_features(input_df):
cols = [
"MeanLight",
"SumLight",
"Year"
]
return input_df[cols].copy()
# 面積(?)
def get_area_feature(input_df):
output_df = pd.DataFrame()
output_df["Area"] = input_df["SumLight"] / (input_df["MeanLight"]+1e-3)
return output_df
# PlaceIDをキーにした集約特徴量
def get_agg_place_id_features(input_df):
_input_df = pd.concat([input_df, get_area_feature(input_df)], axis=1)
group_key = "PlaceID"
group_values = ["MeanLight", "SumLight", "Area"]
agg_methods = ["min", "max", "median", "mean", "std", max_min, q75_q25]
output_df = aggregation(_input_df,
group_key=group_key,
group_values=group_values,
agg_methods=agg_methods)
return output_df
# Year をキーにした集約特徴量
def get_agg_year_features(input_df):
_input_df = pd.concat([input_df, get_area_feature(input_df)], axis=1)
group_key = "Year"
group_values = ["MeanLight", "SumLight", "Area"]
agg_methods = ["min", "max", "median", "mean", "std", max_min, q75_q25]
output_df = aggregation(_input_df,
group_key=group_key,
group_values=group_values,
agg_methods=agg_methods)
return output_df
# PlaceID をキーにしたグループ内差分
def get_diff_agg_place_id_features(input_df):
group_key = "PlaceID"
group_values = ["MeanLight", "SumLight"]
num_diffs = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21]
output_df = diff_aggregation(input_df,
group_key=group_key,
group_values=group_values,
num_diffs=num_diffs)
return output_df
# PlaceID をキーにしたグループ内シフト
def get_shift_agg_place_id_features(input_df):
group_key = "PlaceID"
group_values = ["MeanLight", "SumLight"]
num_shifts = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21]
output_df = shift_aggregation(input_df,
group_key=group_key,
group_values=group_values,
num_shifts=num_shifts)
return output_df
# pivot tabel を用いた特徴量
def get_place_id_vecs_features(input_df):
_input_df = pd.concat([input_df, get_area_feature(input_df)], axis=1)
# pivot table
area_df = pd.pivot_table(_input_df, index="PlaceID", columns="Year", values="Area").add_prefix("Area=")
mean_light_df = pd.pivot_table(_input_df, index="PlaceID", columns="Year", values="MeanLight").add_prefix("MeanLight=")
sum_light_df = pd.pivot_table(_input_df, index="PlaceID", columns="Year", values="SumLight").add_prefix("SumLight=")
all_df = pd.concat([area_df, mean_light_df, sum_light_df], axis=1)
# PCA all
sc_all_df = StandardScaler().fit_transform(all_df.fillna(0))
pca = PCA(n_components=64, random_state=2021)
pca_all_df = pd.DataFrame(pca.fit_transform(sc_all_df), index=all_df.index).rename(columns=lambda x: f"PlaceID_all_PCA={x:03}")
# PCA Area
sc_area_df = StandardScaler().fit_transform(area_df.fillna(0))
pca = PCA(n_components=16, random_state=2021)
pca_area_df = pd.DataFrame(pca.fit_transform(sc_area_df), index=all_df.index).rename(columns=lambda x: f"PlaceID_Area_PCA={x:03}")
# PCA MeanLight
sc_mean_light_df = StandardScaler().fit_transform(mean_light_df.fillna(0))
pca = PCA(n_components=16, random_state=2021)
pca_mean_light_df = pd.DataFrame(pca.fit_transform(sc_mean_light_df), index=all_df.index).rename(columns=lambda x: f"PlaceID_MeanLight_PCA={x:03}")
# PCA SumLight
sc_sum_light_df = StandardScaler().fit_transform(sum_light_df.fillna(0))
pca = PCA(n_components=16, random_state=2021)
pca_sum_light_df = pd.DataFrame(pca.fit_transform(sc_sum_light_df), index=all_df.index).rename(columns=lambda x: f"PlaceID_SumLight_PCA={x:03}")
df = pd.concat([all_df, pca_all_df, pca_area_df, pca_mean_light_df, pca_sum_light_df], axis=1)
output_df = pd.merge(_input_df[["PlaceID"]], df, left_on="PlaceID", right_index=True, how="left")
return output_df.drop("PlaceID", axis=1)
# PlaceIDをキーにしたグループ内相関係数
def get_corr_features(input_df):
_input_df = pd.concat([input_df, get_area_feature(input_df)], axis=1)
group_key = "PlaceID"
group_vlaues = [
["Year", "MeanLight"],
["Year", "SumLight"],
["Year", "Area"],
]
dfs = []
for gv in group_vlaues:
_df = _input_df.groupby(group_key)[gv].corr().unstack().iloc[:, 1].rename(f"Corr={gv[0]}-{gv[1]}")
dfs.append(pd.DataFrame(_df))
dfs = pd.concat(dfs, axis=1)
output_df = pd.merge(_input_df[[group_key]], dfs, left_on=group_key, right_index=True, how="left").drop(group_key, axis=1)
return output_df
# count 63
def get_count63_feature(input_df):
# 各地域でMeanLightが63をとった回数を特徴量にする
_mapping = input_df[input_df['MeanLight']==63].groupby('PlaceID').size()
output_df = pd.DataFrame()
output_df['count63'] = input_df['PlaceID'].map(_mapping).fillna(0)
return output_df
- 今回提供されているデータは、学習用データ・テストデータともにほぼ完全なパネルデータの構造をしています(学習用データには欠損がある場合がある)。
また、学習用データとテストデータには同じPlaceIDが存在しません。 - このことから、ある
PlaceID、つまり特定の場所の特性を特徴量として入れることが重要になりそうです。イメージとしてはPlaceIDごとの分散表現のようなものを与える感じでしょうか。
そのイメージで拵えた特徴量がplace_id_vecs_featuresです。 - この特徴量は
PlaceIDとYearで作成したピボットテーブルが主です。また、それをPCAで次元削減したものです。こうすることで、場所ごとで同一の特徴量を分散表現チックに獲得するとこができるのではないかという思惑です。 -
PlaceIDをキーとしたグループ内の相関係数を特徴量としているのも同様の考えです。
pivot tableで作った特徴量はこんな感じです。これをPlaceIDでマージします
Year MeanLight=1992 MeanLight=1993 MeanLight=1994 MeanLight=1995 MeanLight=1996
PlaceID
1 10.867647 8.985294 10.595589 10.794118 10.404411
2 48.815475 42.107143 50.172619 48.172619 48.398811
4 33.886139 29.480198 35.683167 35.009899 35.272278
8 32.555557 31.166666 35.027779 34.083332 33.861111
10 3.800000 5.400000 7.400000 5.400000 5.800000
最後に上記で定義した関数を実行します。
ここでは学習用データとテストデータを結合してから前処理していますが、気になる方は別々の処理に書き直してください<(_ _)>(PCA等の処理で少しleakyになります)
# 特徴量作成用の関数を実行する関数
def preprocess(train, test):
input_df = pd.concat([train, test]).reset_index(drop=True) # use concat data
# load process functions
process_blocks = [
get_raw_features,
get_area_feature,
get_agg_place_id_features,
get_agg_year_features,
get_diff_agg_place_id_features,
get_shift_agg_place_id_features,
get_place_id_vecs_features,
get_corr_features,
get_count63_feature
]
# preprocess
output_df = []
for func in tqdm(process_blocks):
_df = func(input_df)
output_df.append(_df)
output_df = pd.concat(output_df, axis=1)
# separate train and test
train_x = output_df.iloc[:len(train)]
test_x = output_df.iloc[len(train):].reset_index(drop=True)
train_y = train[TARGET]
return train_x, train_y, test_x
train_x, train_y, test_x = preprocess(train, test)
train_y = np.log1p(train_y) # 目的変数のlogをとる
train & predict
最初に使うモデルの定義とパラメータの設定を行います
今回はお馴染みlightGBMです
model = LGBMModel # 使うモデル
group = train["PlaceID"] # 今回は group kfold を用いる！
# model の parameter
model_params = {
"n_estimators": 10000,
"objective": 'rmse',
"learning_rate": 0.01,
"num_leaves": 31,
"random_state": 2021,
"n_jobs": -1,
"importance_type": "gain",
'colsample_bytree': .5,
"reg_lambda": 5
}
# fit 時の parameter
fit_params = {
"early_stopping_rounds": 100,
"verbose": False
}
次に cross validation 用の関数と評価関数を定義します
今回は、PlaceIDをキーとしたGroup K foldという方法を用います
class GroupKFold:
"""
GroupKFold with random shuffle with a sklearn-like structure
"""
def __init__(self, n_splits=4, shuffle=True, random_state=42):
self.n_splits = n_splits
self.shuffle = shuffle
self.random_state = random_state
def get_n_splits(self, X=None, y=None, group=None):
return self.n_splits
def split(self, X=None, y=None, group=None):
kf = model_selection.KFold(n_splits=self.n_splits, shuffle=self.shuffle, random_state=self.random_state)
unique_ids = group.unique()
for tr_group_idx, va_group_idx in kf.split(unique_ids):
# split group
tr_group, va_group = unique_ids[tr_group_idx], unique_ids[va_group_idx]
train_idx = np.where(group.isin(tr_group))[0]
val_idx = np.where(group.isin(va_group))[0]
yield train_idx, val_idx
# PlaceID をキーにした Group K fold
def make_gkf(X, y, n_splits=5, random_state=2020):
gkf = GroupKFold(n_splits=n_splits, random_state=random_state)
return list(gkf.split(X, y, group))
# 評価関数(log(目的変数)を使っているためrmseを用いる)
def root_mean_squared_error(y_true, y_pred):
return mean_squared_error(y_true, y_pred) ** .5
学習を始める前に feature importance を見てみます
# GBDTのfeature importanceを出力&可視化
def tree_importance(X, y, model, model_params, fit_params, cv, folds):
est = model(**model_params)
feature_importance_df = pd.DataFrame()
for i, (tr_idx, va_idx) in enumerate(cv(X, y, n_splits=folds)):
tr_x, va_x = X.values[tr_idx], X.values[va_idx]
tr_y, va_y = y.values[tr_idx], y.values[va_idx]
est.fit(tr_x, tr_y,
eval_set=[[va_x, va_y]],
**fit_params)
_df = pd.DataFrame()
_df['feature_importance'] = est.feature_importances_
_df['column'] = X.columns
_df['fold'] = i + 1
feature_importance_df = pd.concat([feature_importance_df, _df], axis=0, ignore_index=True)
order = feature_importance_df.groupby('column') \
.sum()[['feature_importance']] \
.sort_values('feature_importance', ascending=False).index[:50]
fig, ax = plt.subplots(figsize=(12, max(4, len(order) * .2)))
sns.boxenplot(data=feature_importance_df, y='column', x='feature_importance', order=order, ax=ax,
palette='viridis')
fig.tight_layout()
ax.grid()
ax.set_title('feature importance')
fig.tight_layout()
plt.show()
return fig, feature_importance_df
fig, _ = tree_importance(train_x, train_y, model, model_params, fit_params, cv=make_gkf, folds=5)
特徴量重要度をみるとplace_id_vecs_featuresで作った特徴量やPlaceIDをキーにした集約特徴量がけっこう効いているようです
学習の段階までたどり着けました。学習用の関数を定義してから、そこに上記で設定したモデルとパラメータをぶち込んで学習を回してみましょう
今回は 5fold, 3 seed average です
# cross validation で学習
def predict_cv(X, y, model, model_params, fit_params, cv, folds, seeds, metrics, name=""):
oof_seeds = []; scores_seeds = []; models = {}
for seed in seeds:
oof = []; va_idxes = []; scores = []
train_x, train_y = X.values, y.values
fold_idx = cv(train_x, train_y, n_splits=folds, random_state=seed)
if "random_state" in model_params:
model_params["random_state"] = seed
# train and predict by cv folds
for cv_num, (tr_idx, va_idx) in enumerate(fold_idx):
tr_x, va_x = train_x[tr_idx], train_x[va_idx]
tr_y, va_y = train_y[tr_idx], train_y[va_idx]
va_idxes.append(va_idx)
est = model(**model_params)
# fitting - train
est.fit(tr_x, tr_y,
eval_set=[[va_x, va_y]],
**fit_params)
model_name = f"{name}_SEED{seed}_FOLD{cv_num}_model"
models[model_name] = est # save model
# predict - validation
pred = est.predict(va_x)
oof.append(pred)
# validation score
score = metrics(va_y, pred)
scores.append(score)
print(f"SEED:{seed}, FOLD:{cv_num} >>>> val_score:{score:.4f}")
# sort as default
va_idxes = np.concatenate(va_idxes)
oof = np.concatenate(oof)
order = np.argsort(va_idxes)
oof = oof[order]
oof_seeds.append(oof)
scores_seeds.append(np.mean(scores))
oof = np.mean(oof_seeds, axis=0)
print(f"model:{name} score:{metrics(train_y, oof):.4f}\n")
return oof, models
# oof予測値と学習済みモデルをGet!
oof, models = predict_cv(train_x, train_y, model, # 説明変数, 目的変数, モデル
model_params, # モデルパラメータ
fit_params, # fit時のパラメータ
cv=make_gkf, # cv strategy
folds=5, # num folds
seeds=[0, 1, 2], # 3 seed average
metrics=root_mean_squared_error, # 評価関数
name="lgbm"
)
SEED:0, FOLD:1 >>>> val_score:0.5399
SEED:0, FOLD:2 >>>> val_score:0.5804
SEED:0, FOLD:3 >>>> val_score:0.5198
SEED:0, FOLD:4 >>>> val_score:0.5376
SEED:1, FOLD:0 >>>> val_score:0.5294
SEED:1, FOLD:1 >>>> val_score:0.5456
SEED:1, FOLD:2 >>>> val_score:0.5112
SEED:1, FOLD:3 >>>> val_score:0.5494
SEED:1, FOLD:4 >>>> val_score:0.5384
SEED:2, FOLD:0 >>>> val_score:0.5494
SEED:2, FOLD:1 >>>> val_score:0.5505
SEED:2, FOLD:2 >>>> val_score:0.5533
SEED:2, FOLD:3 >>>> val_score:0.5280
SEED:2, FOLD:4 >>>> val_score:0.5530
model:lgbm score:0.5353
結果は... CV:0.5353!!でした、どうなんでしょうか
次に学習済みモデルを用いてテストデータを予測してみます
# テストデータを学習済みモデルで予測
def predict_test(test_x, models):
preds = []
for name, est in models.items():
print(f"{name} >>>> ")
_pred = est.predict(test_x)
preds.append(_pred)
preds = np.mean(preds, axis=0)
return preds
# 提出用の予測値をGet!
preds = predict_test(test_x, models)
lgbm_SEED0_FOLD0_model >>>>
lgbm_SEED0_FOLD1_model >>>>
lgbm_SEED0_FOLD2_model >>>>
lgbm_SEED0_FOLD3_model >>>>
lgbm_SEED0_FOLD4_model >>>>
lgbm_SEED1_FOLD0_model >>>>
lgbm_SEED1_FOLD1_model >>>>
lgbm_SEED1_FOLD2_model >>>>
lgbm_SEED1_FOLD3_model >>>>
lgbm_SEED1_FOLD4_model >>>>
lgbm_SEED2_FOLD0_model >>>>
lgbm_SEED2_FOLD1_model >>>>
lgbm_SEED2_FOLD2_model >>>>
lgbm_SEED2_FOLD3_model >>>>
lgbm_SEED2_FOLD4_model >>>>
得られたoof予測値、テストデータ予測値、真の目的変数の値をプロットしてみます
def result_plot(train_y, oof, preds):
name = "result"
fig, ax = plt.subplots(figsize=(8, 8))
sns.distplot(train_y, label='trian_y')
sns.distplot(oof, label='oof')
sns.distplot(preds, label='preds')
ax.legend()
ax.grid()
ax.set_title(name)
fig.tight_layout()
plt.show()
result_plot(train_y, oof, preds)
真の値の分布とは少し違う感じですね
make submission
最後に提出用ファイルを拵えておしまいです。このときテストデータ予測値は対数をとった値になっているので直す必要があります
sample_sub["LandPrice"] = np.expm1(preds) # RMSE → RMSLE
sample_sub.to_csv('../submission/baseline001.csv', index=False)
投稿!!
summary
- 以上が今回の記事の内容です! 何か参考になるところなんかあれば嬉しいです。
- 大きく効いた特徴は、
PlaceIDの場所は一体どんな場所なのかを示せるような特徴量す。つまりPlaceID単位で同じ値をとる特徴量ですね。時間変化の重みを場所ごとで同一な値につけて変化を加える拡張や、モデル自体をパネルデータを扱えるもに変えるなどもできそうです。 - 面積はその土地を表す代表的な指標でもあるので、面積まわりで特徴量を考えてもいいのかもしれないです。
- 学習用データを、欠損値補完をし完全なパネルデータの形にすることも考えられます。
- クラスタリングで似たような地域をグルーピングすることで、新しい特徴量を作ることも可能だと思います。
- CV-LB相関については、自分の提出結果をみると、多少ブレがあるのももの相関はとれていそうです。
参考記事
この記事ではあまり触れていなかったことは以下の記事で詳しい説明があります。
私自身もめちゃくちゃ参考にしました！ありがとうございます！
Discussion