水烏 - Kaggleに挑み続けるブログ

とあるインフラ企業のデータ分析担当がKaggleに挑んだその軌跡

不動産価格を機械学習で予測するKaggleに挑戦する [発展編1:上位20%]

移転しました。

#今回使うライブラリ
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline

from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error
from sklearn.linear_model import Lasso
from sklearn import svm
from sklearn.model_selection import GridSearchCV, train_test_split

前回は2690/5484チーム、スコアも1.5029とほぼど真ん中のスコアで、あまりパッとした結果になりませんでした。
内容的には欠損値の補完で大きな差が出るとも思えず、つまり改善すべきは下記の3点になると考えます。
- Feature Engineeringをせず与えられた特徴量のみ使用した
- 外れ値の確認を行っていない
- 機械学習アルゴリズム。線形と木とSVRでGrid Searchにホイッという手抜きぶり

発展編ではそれぞれ掘り下げて、スコアにどのような影響を与えるのか確認していきたいと思います。
前回の欠損値補完までのステップはスキップしています。get dummiesした後から始まります。

Feature Engineering

前回は与えられた特徴量以外に特に何も使用しませんでしたが、今回は色々作ってみたいと思います。
二乗しまくるとか、総当りで割りまくる、大きいグループでまとめるなど色々なアプローチが考えられますがまずは特徴量を増やしまくって、その後変数増加法で絞る作戦で行きたいと思います。

総当り作戦は面積系だけに留めます。

余談ですがdata_descriptionにもKaggleの説明にもベッドルームの数は'Bedroom'となっていますが、
実際のデータの特徴量名は'BedroomAbvGr'です。こういうの怖いなと思ったり。
しかもKitchenも同様。

あとFeature Engineeringって熱中すると軽く数時間飛んじゃうので、あらかじめ時間を決めた方がメリハリが出て時間の浪費も防げるかもしれません。
今回は1時間。
2乗3乗するのは、そうすることでデータ分布の形状が変わり、新たなパターンが見つかることを期待してのことです。

#日本人なら気になる、「で、何LDKなの?」
df_all['xLDK'] = df_all['BedroomAbvGr'] + df_all['KitchenAbvGr']

#何年に建ったかより、築年数の方が重要のはず
df_all['YearSinceBuilt'] = 2018 - df_all['YearBuilt']

#何年に売れたかより、売れてから経過している年数の方が重要のはず
df_all['YrSinceSold'] = 2018 - df_all['YrSold']

#ガレージの築年数
df_all['GarageSinceYrBlt'] = 2018 - df_all['GarageYrBlt']

#大きなくくりで集計してみる
df_all['TotalHousePorchSF'] = df_all['EnclosedPorch'] + df_all['OpenPorchSF'] + df_all['WoodDeckSF']+df_all['3SsnPorch']+df_all['ScreenPorch']
df_all['TotalHouseSF'] = df_all['1stFlrSF'] + df_all['2ndFlrSF'] + df_all['TotalBsmtSF']

#なんかLow Qualityな面積があるらしいので全体の面積から引く
df_all['TotalHouseSFHighQuality'] = df_all['TotalHouseSF'] - df_all['LowQualFinSF']

#割るぞ~
df_all['1stFlr_vs_2ndFlr'] = df_all['1stFlrSF'] / (df_all['2ndFlrSF'] + 0.001)
df_all['Bsmt_vs_1stAnd2ndFlr'] = df_all['TotalBsmtSF'] / (df_all['1stFlrSF'] + df_all['2ndFlrSF'] + 0.001)

#Overall Qualityが結構重要そうなので2乗と3乗してみる
df_all['OverallQual_2'] = df_all['OverallQual']**2
df_all['OverallQual_3'] = df_all['OverallQual']**3

高次元データにおける外れ値処理

前回は特に外れ値等を確認しませんでした。
79もの特徴量が存在し、家の価格という比較的身近な題材でもあり主観が入りやすくなると考えたためです。
他にも単変量でプロットして外れ値だと解釈して削除しても、実は別の特徴量も込みで考えれば外れ値ではなかったということも考えられます。

よって多変量で考える必要があり、多変量の外れ値検知手法としてはModified Stahel Donoho法などいくつかの方法が知られています。
今回はIsolation Forestと呼ばれる手法を使って外れ値の存在を確かめたいと思います。

Isolation Forestは名前が示す通り決定木ベースの外れ値検知法で、特徴量をランダムに選び、決定境界を適当に何度も引きます。
その後、あるサンプルを集団から分離させるのに必要だった分割回数(線を引いた回数)を利用して、外れ値を検出します。
外れ値ほど分離しやすく分割回数は減り、データが密集している付近の値ほど分離しにくく分割回数が増えるという具合です。

決定木ベースなので、ランダムフォレストなどと同じく正規化等なくてもうまい具合に動いてくれる上、データが多次元正規分布に従っているなどの仮定も必要ないのでとても便利に扱えます。
ツリーベースのアルゴリズムはこの安定・安心感が素晴らしいと思います。

尚、Isolation Forestはみんな大好きsklearnに実装されているので、簡単に呼べます。便利!
http://scikit-learn.org/stable/modules/generated/sklearn.ensemble.IsolationForest.html

#学習用データとテストデータを切り分ける
ntrain = df_train.shape[0]
train = df_all[:ntrain]
test = df_all[ntrain:]

y = y_train
#X = train.loc[:, train.columns != 'SalePrice']
#Isolation Forest
from sklearn.ensemble import IsolationForest
clf = IsolationForest(random_state=1234,n_estimators=500, contamination=0.01)
clf.fit(train)
IsolationForest(bootstrap=False, contamination=0.01, max_features=1.0,
        max_samples='auto', n_estimators=500, n_jobs=1, random_state=1234,
        verbose=0)
mv_outliers = pd.DataFrame(clf.predict(train))
mv_outliers.columns = ['OutlierFlag']
mv_outliers[mv_outliers == -1].count()
OutlierFlag    15
dtype: int64

Contaminationは外れ値の割合で、とりあえず1%にセットしています。
上記の15点について、まずどのようなデータが選ばれたのか確認したいと思います。
とりあえずSalePriceと1階の広さをPlotして確認します。

y_train = pd.DataFrame(y_train)
y_train.columns = ['SalePrice']
df_all_olcheck = pd.concat([mv_outliers,train, y_train], axis=1)
ax = sns.lmplot(x='1stFlrSF', y='SalePrice',fit_reg=False,
                data=df_all_olcheck, hue='OutlierFlag', size=7, aspect=1.5)

f:id:prml2synth:20180604100212p:plain

青く塗られているのが外れ値扱いされた点で、オレンジが正常値です。
もちろん多変量での外れ値検知のため、1階の広さと価格プロットのみでは見切れない何かが考慮されているのでしょう。

確かに明らかに外れている右のほうの2点も青く塗られているため、アルゴリズム自体は正しく働いている感触です。
ではこの外れ値を取り除いて、前回と同様にスコアを確認してみましょう。

df_all_r1 = df_all_olcheck[df_all_olcheck['OutlierFlag'] == 1]
df_all_r1.drop("OutlierFlag", axis=1, inplace=True)
y_train_r1 = df_all_r1["SalePrice"].values
df_all_r1.drop(["SalePrice"], axis=1, inplace=True)
#30%でチューニングを行う
X_train, X_test, y_train, y_test = train_test_split(df_all_r1, y_train_r1, test_size=0.3, random_state=1234)
#とりあえずランダムフォレスト、パラメーターは前回のグリッドサーチ結果と同様
rf = RandomForestRegressor(n_estimators = 500, max_depth=10, random_state=1234)
rf.fit(X_train,y_train)
RandomForestRegressor(bootstrap=True, criterion='mse', max_depth=10,
           max_features='auto', max_leaf_nodes=None,
           min_impurity_decrease=0.0, min_impurity_split=None,
           min_samples_leaf=1, min_samples_split=2,
           min_weight_fraction_leaf=0.0, n_estimators=500, n_jobs=1,
           oob_score=False, random_state=1234, verbose=0, warm_start=False)
#RMSE
y_pred2 = rf.predict(X_test)
print("ランダムフォレストでのRMSE:",np.sqrt(mean_squared_error(y_test, y_pred2)))
ランダムフォレストでのRMSE: 24143.2583576

前回よりもRMSEが大きくなってしまいました。
もちろんデータが違うので同一次元では語れません。Kaggleにアップロードして確認します。

y_pred_final = rf.predict(test)
submission = pd.concat((df_test_index, pd.DataFrame(y_pred_final)), axis=1)
submission.columns = ['Id', 'SalePrice']
submission.to_csv("submission.csv",sep=',',index=False)

結果、スコアは0.15051と前回とほぼ変わらず。
外れ値の処理を行っても、多少Featureを追加してもあまり変わらないようですね。

もちろん、無駄かどうかはまだ分かりません。別の機械学習アルゴリズムでは精度が上がるかもしれません。
ここまで確認できたら、後は機械学習のアルゴリズムを試しまくります!

様々な機械学習での予測

今のところランダムフォレストでしかスコアを試していません。
本当はL1/L2正則化回帰、ツリーベース、SVRなど色々と試すのがバイアス・バリアンスの観点から良さそうです。

今回は前述のランダムフォレスト、ラッソ回帰、に加えxgboost(勾配ブースティング)を試してみようと思います。

import xgboost as xgb
xgb_reg = xgb.XGBRegressor(max_depth=5, n_estimators=500, subsample=0.8)
xgb_reg.fit(X_train,y_train)
XGBRegressor(base_score=0.5, booster='gbtree', colsample_bylevel=1,
       colsample_bytree=1, gamma=0, learning_rate=0.1, max_delta_step=0,
       max_depth=5, min_child_weight=1, missing=None, n_estimators=500,
       n_jobs=1, nthread=None, objective='reg:linear', random_state=0,
       reg_alpha=0, reg_lambda=1, scale_pos_weight=1, seed=None,
       silent=True, subsample=0.8)
#RMSE
y_pred3 = xgb_reg.predict(X_test)
print("xgboostでのRMSE:",np.sqrt(mean_squared_error(y_test, y_pred3)))
xgboostでのRMSE: 21661.4138212

xgboostに切り替えた瞬間RMSEが2000以上下がりました。
これがブースティングの力なのか...

xgboostにした途端0.13564で2000位以内に入りました。騒がれるわけです。これはスゴイ。
しかもRandom Forestと同じように、特徴量の重要度を確認することもできます。
300を超える特徴量がありますが、ちょっと上位20個を確認してみましょう。

from xgboost import plot_importance
fig, ax = plt.subplots(1,1,figsize=(8,7))
plot_importance(xgb_reg, ax=ax, max_num_features=20)
plt.show()

f:id:prml2synth:20180604100219p:plain

Feature Engineeringで作成した"TotalHousePorchSF"が約330ある変数のうち重要度6位に!これは結構うれしいですね。
他にもTotalHouseSFやBsmt_vs_1stAnd2ndFlrがランクインしていて、Feature Engineeringやった甲斐があるというものです。
上記3つは家のことを考えながら作成した特徴量で、やはり適当に2乗3乗するより効くんだと実感。

にしても、LotFrontageとLotAreaがワンツーというのは面白いですね。
上位に入ってくれたTotalHousePorchSFも家のアプローチ全面積という意図ですし、家の中より外なのか。。。

ところで、私は大のSVM好きです。あの滑らかな決定境界には毎回ため息が漏れます。
ということでSVRもやりたいと思います。
SVRは結構Sensitiveなので、きちんと前処理を行う必要があります。

from sklearn.preprocessing import StandardScaler, RobustScaler, MinMaxScaler
scaler = MinMaxScaler()
X_train_svr = scaler.fit_transform(X_train)
X_test_svr = scaler.transform(X_test)
#SVR
svr = svm.SVR(C=1e6,epsilon=0.2,kernel='rbf')
svr.fit(X_train_svr, y_train)
#RMSE
y_pred = svr.predict(X_test_svr)
print("SVRでのRMSE:",np.sqrt(mean_squared_error(y_test, y_pred)))
SVRでのRMSE: 19752.2212481

・・・あれっ?SVRこんなに精度出るんだ!?
ちょっと気になったのでSVR単体でアップロードしてみます。

test_svr = scaler.transform(test)
y_pred_final = svr.predict(test_svr)
submission = pd.concat((df_test_index, pd.DataFrame(y_pred_final)), axis=1)
submission.columns = ['Id', 'SalePrice']
submission.to_csv("submission.csv",sep=',',index=False)

f:id:prml2synth:20180604101741p:plain なんとスコアが0.12611、順位は1303位まで上昇しました。上位20%は目前です。
あまりHouse PricesでSVRを使う人を見ない気がしますが、やっぱSVR強いですね。

スタッキング

更に、Kaggleで流行っているスタッキングという手法を使いたいと思います。

スタッキングはアンサンブル学習の一種とも考えられ、複数の学習器を組み合わせて予測を行うことです。

具体的には予測を段階にわけて、第一段階は各種機械学習のアルゴリズムを出力し、第二段階の入力としてその出力を別の機械学習アルゴリズムに流し込んで予測するという具合です。
第二段階は機械学習アルゴリズムを通しても良いし、単に平均、多数決、重み付けなどでも問題ないようです。

スタッキングで精度が上がる理由としては、それぞれの学習器とチューニングによるバイアスとバリアンスがうまいことならされる・・・ことかなと考えますが、アンサンブル学習の弱学習器が強い奴になった感じというざっくり理解で使います。

今回は、今まで試したxgboost, SVRとラッソ回帰の3つを使ってスタッキングしてみたいと思います。

#まだ予測値を取得していないラッソ回帰
lasso = Lasso(normalize=True, random_state=123, max_iter=10000)
lasso.fit(X_train,y_train)
Lasso(alpha=1.0, copy_X=True, fit_intercept=True, max_iter=10000,
   normalize=True, positive=False, precompute=False, random_state=123,
   selection='cyclic', tol=0.0001, warm_start=False)
#RMSE
y_pred4 = lasso.predict(X_test)
print("LassoでのRMSE:",np.sqrt(mean_squared_error(y_test, y_pred4)))
LassoでのRMSE: 21573.4611747

ラッソも悪くないですね。

まずは学習用データを使って、それぞれの結果をどれぐらいの割合で混ぜるのが良いか確認します。

#SVRはy_pred, xgboostはy_pred3, Lassoはy_pred4
y_stack = y_pred*(1/3) + y_pred3*(1/3) + y_pred4*(1/3)
print("StackingでのRMSE:",np.sqrt(mean_squared_error(y_test, y_stack)))
StackingでのRMSE: 18586.9432485

単純に平均しただけでもSVRを超えてしまいました。
この割合を使ってテストデータをスタッキングし、Kaggleに提出してみます。

y_pred_lasso = lasso.predict(test)
y_pred_xgb = xgb_reg.predict(test)
test_svr = scaler.transform(test)
y_pred_svr = svr.predict(test_svr)

y_stack = y_pred_svr*(1/3) + y_pred_xgb*(1/3) + y_pred_lasso*(1/3)

スコアは0.12373、順位は1152位で着実にスコアは上がっています。
ここからもう一歩スコアを伸ばしたいところです。

とりあえず、xgboostと並んで人気のlightBGMを使ってみます。

import lightgbm as lgb
lgb_train = lgb.Dataset(X_train,y_train)
lgb_eval = lgb.Dataset(X_test,y_test)
params = {
    'task': 'train',
    'boosting_type': 'gbdt',
    'objective': 'regression',
    'metric': {'l2'},
    'num_leaves': 256,
    'learning_rate': 0.005,
    'num_iterations':3000,
    'feature_fraction': 0.4,
    'bagging_fraction': 0.7,
    'bagging_freq': 5,
}
gbm = lgb.train(params,lgb_train,num_boost_round=1500,valid_sets=lgb_eval)
 
y_test_pred_lgb = gbm.predict(np.array(X_test), num_iteration=gbm.best_iteration)

#RMSE
print("LightGBMでのRMSE:",np.sqrt(mean_squared_error(y_test, y_test_pred_lgb)))
LightGBMでのRMSE: 20927.1470379
y_stack = y_pred3*0.2+y_pred*0.55+y_test_pred_lgb*0.25
print("StackingでのRMSE:",np.sqrt(mean_squared_error(y_test, y_stack)))
StackingでのRMSE: 18524.2845651

わずかですが、LightBGMなしの頃よりもRMSEが低下しています。
では、これでKaggleにSubmitしてみましょう。

y_pred_lasso = lasso.predict(test)
y_pred_xgb = xgb_reg.predict(test)
test_svr = scaler.transform(test)
y_pred_svr = svr.predict(test_svr)
y_pred_lgb = gbm.predict(np.array(test), num_iteration=gbm.best_iteration)
y_stack = y_pred_xgb*(0.2)+y_pred_svr*(0.55)+y_pred_lgb*(0.25)
submission = pd.concat((df_test_index, pd.DataFrame(y_stack)), axis=1)
submission.columns = ['Id', 'SalePrice']
submission.to_csv("submission.csv",sep=',',index=False)

f:id:prml2synth:20180604102014p:plain ついにスコアは0.11968、順位は844位まで浮上しました。
上位20%を突破し15%に迫る勢いです。

今回は外れ値の処理や正規化処理(Lasso)、Feature Engineeringにスタッキングと一通りのテクニックを試しました。
スタッキングは、個別の機械学習では限界を感じている時に試してみると精度が上がることが多いです。
特に決定境界が異なる(例えば軸直行しか引けないRandom Forest v.s. SVM)ものを組み合わせると、汎化性能の観点で好ましいと思われます。

しかし、オッカムの剃刀のようにモデルがどんどん複雑化することを好まないとする意見も根強くあります。 例えばNetflix Prizeで1位になったモデルが数十個のモデルを複数層にスタッキングしたもので、その運用コストから結局使われなかった例があります。
こちらはコストの話なので本筋とは少しずれるかも知れませんが、やはり可能な限りシンプルかつ説明可能なデータ処理・予測モデル構築で済ませるのがベストでしょう。

次回は単独の機械学習手法を使って更にスコアを伸ばすことができるか挑戦したいと思います。

上位3%以上を狙いたい方は以下を参考ください。

www.mizugarasu.jp