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

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

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

前回のは色々複雑すぎる...

不動産価格を機械学習で予測するKaggleのHouse Pricesですが、前回のスコアは0.12を切り、順位も上位20%に入ることができました。
引き換えに色々な試行錯誤、処理、機械学習をテストしたため、chain of thoughts, モデルが複雑になってしまいました。。

今回はできるだけシンプルなデータ処理と予測モデリングでより高い順位を目指したいと思います。
シンプルな予測モデルといえば、やはり線形回帰です。今回はラッソ回帰を中心に考えてみたいと思います。

import seaborn as sns
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

from scipy.stats import skew
from sklearn.linear_model import Lasso
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error

import lightgbm as lgb

今回は上記のライブラリ以外使いません。途中から新しくライブラリをインポートすることも行いません。
これで果たして上位に行けるのか、見てみましょう。今回は、できるだけシンプルにやります。

#データの読み込み
df_train = pd.read_csv('train.csv')
df_test = pd.read_csv('test.csv')

#外れ値の処理
df_train = df_train[~((df_train['GrLivArea'] > 4000) & (df_train['SalePrice'] < 300000))]

#Id及びSaleprice以外の全データを結合(左端=Id, 右端はSalePrice)
df_all = pd.concat((df_train.loc[:,'MSSubClass':'SaleCondition'], df_test.loc[:,'MSSubClass':'SaleCondition']))

#Feature Engineering - 前回の中で断トツでFeature Importanceが高かったTotalHousePorchSFを採用
df_all['TotalHousePorchSF'] = df_all['EnclosedPorch']+df_all['OpenPorchSF']+df_all['WoodDeckSF']+df_all['3SsnPorch']+df_all['ScreenPorch']
df_train['TotalHousePorchSF'] = df_train['EnclosedPorch']+df_train['OpenPorchSF']+df_train['WoodDeckSF']+df_train['3SsnPorch']+df_train['ScreenPorch']

外れ値も、前回わかっていた広大なエリアに対して破格の値段がついている(恐らくは何もない土地)の2点だけを除きました。

回帰分析の基礎に立ち返る

線形回帰分析を行う上で、残差の正規性や多重共線性の有無、説明変数ができるだけ正規分布している、線形の関係性があるなど気にすることがあります。
L1正則化項の存在で回帰係数が不安定になりづらいとはいえ、多重共線性の有無、説明変数・目的変数の分布の正規化ぐらいはやっておこうと思います。
尚、不動産価格に対して各種特徴量が概ね線形を保っているのはHouse Pricesについては自明なので、特に触れません。

多重共線性の確認

多重共線性を確認する方法はVIFや相関係数行列を目視するなどいくつかあります。
今回は相関係数行列を適当に眺めて相関が高い奴(0.8+)を省きます。カテゴリカル変数の共線性の評価はかなり複雑になるのでやりません...

#相関係数行列の可視化
fig, ax = plt.subplots(1, 1, figsize=(30, 30))
sns.heatmap(df_train.corr(), vmax=1, vmin=-1, center=0, annot=True, ax=ax)

f:id:prml2synth:20180610104857p:plain

相関係数が絶対値で0.8を超えているものはほぼ確実に多重共線性を引き起こすので、特徴量から削除します。4ペアあります
落とす際は、目的変数であるSalePriceに対しての相関係数が低いものから順に落としていきます。
この場合落とすのは1stFlrSF, TotRmsAbdGrd, GarageArea, GarageYrBltの4つを落とします。

df_all.drop(['1stFlrSF','GarageArea','TotRmsAbvGrd', 'GarageYrBlt'], axis=1, inplace=True)

説明変数・目的変数の対数変換

次に、SalePriceのヒストグラムを見てみます。

df_train["SalePrice"].hist(bins=30)

f:id:prml2synth:20180610104402p:plain

かなり歪度が高いですね。特にこのように高いほうに向かって裾野が広い分布は、対数変換をかけてやるとうまい具合に正規分布してくれます。
ちゃんとQ-Q Plotやシャピロウィルク検定などで正規性の確認をしても良いと思います。

df_train["SalePrice"] = np.log1p(df_train["SalePrice"])

#対数変換後の分布を確認
df_train["SalePrice"].hist(bins=30)

f:id:prml2synth:20180610104415p:plain

それっぽくなりました。尚、対数変換をかける場合は最後の予測の際にEXPを取ることを忘れないようにする必要があります。
同じことをすべてのノンカテゴリカルの特徴量に対して行います。scipyのskewで歪度が計算できるので、
今回は0.6より歪度が大きいものはすべて対数変換をかけることにします。やや雑ですが、シンプルにしておきたいので。

#カテゴリカルでない特徴量
non_categoricals = df_all.dtypes[df_all.dtypes != "object"].index

skewed_feats = df_train[non_categoricals].apply(lambda x: skew(x.dropna()))
skewed_feats = skewed_feats[skewed_feats > 0.6].index

#歪度0.7より大きい特徴量を対数変換
df_all[skewed_feats] = np.log1p(df_all[skewed_feats])

あとはOne Hot Encodingと、欠損値の補完です。最近気づいたんですがHouse Pricesについては欠損値は
全部平均で補完しても精度に影響ないようなので、今回はシンプルにすべて平均値で補完します。

#One Hot Encoding
df_all = pd.get_dummies(df_all)

#欠損値を平均値で補完
df_all = df_all.fillna(df_all.mean())

今回は、バイアス・バリアンスの観点からラッソ回帰と勾配ブースティングの二つのスタッキングを行います。

#学習データ、テストデータに分割
X = df_all[:df_train.shape[0]]
X_for_test = df_all[df_train.shape[0]:]
y = df_train.SalePrice

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=1111)
reg = Lasso(alpha=0.0004)
reg.fit(X_train, y_train)
Lasso(alpha=0.0004, copy_X=True, fit_intercept=True, max_iter=1000,
   normalize=False, positive=False, precompute=False, random_state=None,
   selection='cyclic', tol=0.0001, warm_start=False)
y_pred = reg.predict(X_test)
print("ラッソ回帰でのRMSE:",np.sqrt(mean_squared_error(y_pred, y_test)))
ラッソ回帰でのRMSE: 0.11467346011593593
lgb_train = lgb.Dataset(X_train,y_train)
params = {'task': 'train','boosting_type': 'gbdt','objective': 'regression','metric': {'l2'},'num_leaves': 256,
          'learning_rate': 0.01,'num_iterations':2000,'feature_fraction': 0.4,'bagging_fraction': 0.7,'bagging_freq': 5}

gbm = lgb.train(params, lgb_train, num_boost_round=1500)
 
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: 0.11716578498521343
print("LightGBM+Lassoのスタッキング時のRMSE:",np.sqrt(mean_squared_error(y_test, (y_pred*0.7+y_test_pred_lgb*0.3))))
LightGBM+Lassoのスタッキング時のRMSE: 0.11062571022201569

いくつかテストした結果これがもっとも良さそうなパラメーターでした。
グリッドサーチやCVも熱心には行っていませんが、これでスコアを見てみたいとおもいます。

#全データで学習
reg.fit(X, y)
lgb_train_full = lgb.Dataset(X,y)
gbm = lgb.train(params, lgb_train_full, num_boost_round=1500)
#ラッソ・LightGBMの予測及びスタッキング
pred = np.expm1(reg.predict(X_for_test))
pred2 = np.expm1(gbm.predict(X_for_test))
pred3 = (pred*0.7+pred2*0.3)
solution = pd.DataFrame({"id":df_test.Id, "SalePrice":pred3})
solution.to_csv("test_lasso.csv", index = False)

f:id:prml2synth:20180610104514p:plain

スコアは0.11482、順位は165位/5486位で上位3%にまで食い込みました。
そこまで複雑でない前処理・Feature Engineering、ラッソ回帰と勾配ブースティングのみ使ってシンプルなスタッキングを行うだけでもここまで上がれることを理解頂けたかと思います。
第3回となったHouse Pricesですが、この辺で一度終わりにして次の題材に進みたいと思います。

もし上記を元にスコアを更に上げたい!という場合はやはりIsolationForestなど多変量での外れ値処理、
深いFeature Engineering及びモデルの複雑化(アルゴリズムの追加と複層のスタッキング)などが必須になってくるでしょう。
NetFlixの例を考えても、実務的には上記ぐらいが運用の限界かなと思います。

不動産価格を機械学習で予測する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

不動産価格を機械学習で予測するKaggleに挑戦する [ベンチマーク編]

第一回ということで、Kaggleの不動産価格予測コンペ - House Prices: Advanced Regression Techniquesに挑戦してみたいと思います。
House Prices: Advanced Regression Techniques | Kaggle
タイタニックが分類問題のチュートリアルならば、House Pricesは回帰問題のチュートリアルという立ち位置でしょう。

個人的には分類問題の方が好きですが、早速本題に入ります。

私の場合はデータ分析の流れとして、まず最低限の処理(欠損値及びカテゴリカルデータの処理)だけを行ってベンチマークの精度とし、その後掘り下げるというアプローチをよくとります。
今回はベンチマーク編ということで、最低限の処理のみを行ってどのくらいのスコアになるかを確認したいと思います。

1. データ分析の目的

まず今回のデータ分析の目的を確認します。
家を買うときに何LDKなのか、何坪なのかといったデータはもちろん不動産価格に影響するでしょう。

しかし今回はそれら以外にも多くの要因が不動産価格に影響することを確認します。
具体的には79の変数を元に、不動産価格を予測することになります。データはアメリカ合衆国アイオワ州エイムズの不動産です。

#今回使うライブラリ
import numpy as np
import pandas as pd

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

2. データセットについて

SalePriceが予測したい不動産価格で、それら以外に1階の広さや地下の天井の高さ、基礎の種類からプールの有無まで多数の特徴量が含まれています。
全部が全部不動産に影響することはないでしょうが、常識的には部屋の広さ、高さ、築年数などは大きく影響するはずです。
特徴量としては面積のように連続値のものもあれば、品質のようなカテゴリカル変数も存在しています。

ベンチマーク編では特徴量同士のデータ可視化や外れ値の処理などは行いません。
79もの特徴量があるため、いくつか重要と考えられる変数だけ抽出したとしてもそこには主観が入ってしまいます。
例えばSalePriceと1回の広さ1stFlrSFをプロットして、外れている点を削除するというのも単変量しか見ておらず、特徴量同士の相互作用を考慮できていません。

データはtrain.csvとtest.csvに分かれており、またKaggleにアップロードする際の形式としてsample_submission.csvが用意されています。
データの説明(Data_Description.txt)も用意されており、必ず読む必要があります。79変数あるのでここでは省きますが、欠損値の補完ですぐに必要になります。

#学習データとテストデータの読み込み
df_train = pd.read_csv("train.csv")
df_test = pd.read_csv("test.csv")

#予測に使わないIdを保持・データフレームから削除
df_train_index = df_train["Id"]
df_test_index = df_test["Id"]
df_train.drop(["Id"], axis=1, inplace=True)
df_test.drop(["Id"], axis=1, inplace=True)
#目的変数であるSalePriceを別に取っておく
y_train = df_train["SalePrice"].values
df_train.drop(["SalePrice"], axis=1, inplace=True)

#学習用データとテストデータを一度統合する
df_all = pd.concat((df_train,df_test)).reset_index(drop=True)
#データ数
print(df_train.shape)
print(df_test.shape)
(1460, 79)
(1459, 79)

データの数は学習用データ・テストデータともにほぼ同数。学習用データが1500レコード弱で比較的扱いやすいですね。

3. 欠損値の補完

データセットの中身を見てみると、結構データが存在しない(NaN)レコードが多いことに気が付きます。これは本当に欠損しているのか、あるいは単純に"該当なし"という意味なのか紐解きながら補完する必要があります。
いずれかのデータが欠損しているレコードを削除してしまうというのは一つの手ですが、可能な限り手元にあるデータを生かすためまずは補完します。

#欠損値の個数確認 / 2919データ
df_all.isnull().sum()[df_all.isnull().sum() != 0].sort_values(ascending=False)
PoolQC          2909
MiscFeature     2814
Alley           2721
Fence           2348
FireplaceQu     1420
LotFrontage      486
GarageFinish     159
GarageYrBlt      159
GarageQual       159
GarageCond       159
GarageType       157
BsmtExposure      82
BsmtCond          82
BsmtQual          81
BsmtFinType2      80
BsmtFinType1      79
MasVnrType        24
MasVnrArea        23
MSZoning           4
BsmtFullBath       2
BsmtHalfBath       2
Utilities          2
Functional         2
Exterior2nd        1
Exterior1st        1
SaleType           1
BsmtFinSF1         1
BsmtFinSF2         1
BsmtUnfSF          1
Electrical         1
KitchenQual        1
GarageCars         1
GarageArea         1
TotalBsmtSF        1
dtype: int64

こうしてみると結構な特徴量で欠損値があるようですが、例えばPoolQCはData Descriptionを読むとPoolQCがNAならプールなしという意味のようです。
ほかにもガレージ系でそれぞれ159個ずつ仲良く欠損していたりするので、値がないというより、"NA"のつもりが欠損値になっているのでしょう。
それでは早速上記リストとData Descriptionを見ながら欠損値を補完していきましょう。

#欠損値の補完

#以下はNaN = NAかNoneの特徴量リスト。よって欠損値をそれぞれNAとNoneで補完する。
df_all["PoolQC"].fillna('NA', inplace=True)
df_all["MiscFeature"].fillna('None', inplace=True)
df_all["Alley"].fillna('NA', inplace=True)
df_all["Fence"].fillna('NA', inplace=True)
df_all["FireplaceQu"].fillna('NA', inplace=True)
df_all["GarageQual"].fillna('NA', inplace=True)
df_all["GarageFinish"].fillna('NA', inplace=True)
df_all["GarageCond"].fillna('NA', inplace=True)
df_all["GarageType"].fillna('NA', inplace=True)
df_all["BsmtCond"].fillna('NA', inplace=True)
df_all["BsmtExposure"].fillna('NA', inplace=True)
df_all["BsmtQual"].fillna('NA', inplace=True)
df_all["BsmtFinType2"].fillna('NA', inplace=True)
df_all["BsmtFinType1"].fillna('NA', inplace=True)
df_all["MasVnrType"].fillna('None', inplace=True)

#以下はNaN = 0の特徴量リスト。例えば地下なら、地下がないんだから0。みたいな。
df_all["GarageYrBlt"].fillna(0, inplace=True) # ガレージ築年数を0にするのも不思議な気はしますが、そもそもガレージがないので他に妥当な数字が思いつかず。
df_all["MasVnrArea"].fillna(0, inplace=True)
df_all["BsmtHalfBath"].fillna(0, inplace=True)
df_all["BsmtFullBath"].fillna(0, inplace=True)
df_all["TotalBsmtSF"].fillna(0, inplace=True)
df_all["BsmtUnfSF"].fillna(0, inplace=True)
df_all["BsmtFinSF2"].fillna(0, inplace=True)
df_all["BsmtFinSF1"].fillna(0, inplace=True)
df_all["GarageArea"].fillna(0, inplace=True)
df_all["GarageCars"].fillna(0, inplace=True)

#欠損レコード数が少なく、大半が一つの値をとっているためあまりに予測の役に立たなさそうな特徴量は単純に最頻値を代入
df_all["MSZoning"].fillna('RL', inplace=True)
df_all["Functional"].fillna('Typ', inplace=True)
df_all["Utilities"].fillna("AllPub", inplace=True)
df_all['SaleType'] = df_all['SaleType'].fillna(df_all['SaleType'].mode()[0])
df_all['Exterior2nd'] = df_all['Exterior2nd'].fillna(df_all['Exterior2nd'].mode()[0])
df_all['Exterior1st'] = df_all['Exterior1st'].fillna(df_all['Exterior1st'].mode()[0])
df_all['KitchenQual'] = df_all['KitchenQual'].fillna(df_all['KitchenQual'].mode()[0])
df_all['Electrical'] = df_all['Electrical'].fillna(df_all['Electrical'].mode()[0])

#LotFrontage - Linear feet of street connected to property
#これは補完方法が明らかかつ簡単で、近くのStreet名=Neighborhoodでグループし平均を取れば良い精度で補完できそう。
f = lambda x: x.fillna(x.mean())
df_all["LotFrontage"] = df_all.groupby("Neighborhood")["LotFrontage"].transform(f)
#欠損値がすべて補完されているか確認
df_all.isnull().sum()[df_all.isnull().sum() != 0].sort_values(ascending=False)
Series([], dtype: int64)

欠損レコードが消えていますね。

かなり地道な作業になりますが、特に欠損レコード数の多いものに対してどう補完するかで他のデータサイエンティストと差別化できます。
平均値が良いのか最頻値が良いのかなどはトライ&エラーしてみるのが一番簡単です。
もちろん離散値の場合に平均値を使うのは感覚的に変だとか、そういった判断をもとに決めるのも良いと思います

4.前処理

すぐに機械学習したいところですが、カテゴリー変数の処理を行います。

カテゴリー変数の扱い方はLabel EncodingやOne Hot Encodingなど様々な方法が存在します。
今回は大したレコード数ではないことと、特徴量の数も1000すらいかないと思われるのでOne Hot Encodingでダミー変数化します。

#One Hot Encoding
df_all = pd.get_dummies(df_all)
#特徴量数の確認
df_all.shape
(2919, 302)

学習データ・テストデータ合わせて3000レコード、300特徴量です。ややスパースなデータになっていそうですが、ベンチマーク編ではこれ以上の処理は行いません。

5. 機械学習で不動産価格を予測する

楽しいパートがやってきました。みんな大好き機械学習です。
ベンチマーク編なので、とりあえずラッソ回帰、ランダムフォレスト、SVRでやろうと思います。

チューニングはとりあえずざっくり辺りをつけてグリッドサーチします。

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

y = y_train
X = train.loc[:, train.columns != 'SalePrice']

#30%でチューニングを行う
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=1234)
#モデルの呼び出し
lasso = Lasso()
rf = RandomForestRegressor()
svr = svm.SVR()
#グリッドサーチ用パラメータの設定
lasso_parameters = {'alpha':[0.1, 0.5, 1]}
rf_parameters= {'n_estimators':[100, 500, 2000], 'max_depth':[3, 5, 10]}
svr_parameters = {'C':[1e-1, 1e+1, 1e+3], 'epsilon':[0.05, 0.1, 0.3]} 

#グリッドサーチ
lasso_gs = GridSearchCV(lasso, lasso_parameters)
lasso_gs.fit(X_train,y_train)

rf_gs = GridSearchCV(rf, rf_parameters)
rf_gs.fit(X_train,y_train)

svr_gs = GridSearchCV(svr, svr_parameters)
svr_gs.fit(X_train,y_train)
GridSearchCV(cv=None, error_score='raise',
       estimator=SVR(C=1.0, cache_size=200, coef0=0.0, degree=3, epsilon=0.1, gamma='auto',
  kernel='rbf', max_iter=-1, shrinking=True, tol=0.001, verbose=False),
       fit_params=None, iid=True, n_jobs=1,
       param_grid={'C': [0.1, 10.0, 1000.0], 'epsilon': [0.05, 0.1, 0.3]},
       pre_dispatch='2*n_jobs', refit=True, return_train_score='warn',
       scoring=None, verbose=0)

それでは、予測してみたいと思います。

#ラッソ回帰
y_pred = lasso_gs.predict(X_test)
print("ラッソ回帰でのRMSE:",np.sqrt(mean_squared_error(y_test, y_pred)))

#ランダムフォレスト
y_pred2 = rf_gs.predict(X_test)
print("ランダムフォレストでのRMSE:",np.sqrt(mean_squared_error(y_test, y_pred2)))

#SVR
y_pred3 = svr_gs.predict(X_test)
print("SVRでのRMSE:",np.sqrt(mean_squared_error(y_test, y_pred3)))
ラッソ回帰でのRMSE: 38914.187417202746
ランダムフォレストでのRMSE: 24883.658315015313
SVRでのRMSE: 72294.23251731148

安定のランダムフォレストですね。今回のデータセットではやはり正規化などの処理を行わないとSVRは辛そうです。
では、実際にどのような値になっているか確認します。

df_visualize = pd.concat((pd.DataFrame(y_test), np.round(pd.DataFrame(y_pred))), axis=1)
df_visualize.columns = ['実際の値','予測値']
df_visualize.head(10)
実際の値 予測値
0 205000 203899.0
1 345000 355549.0
2 173900 183026.0
3 93500 45687.0
4 265900 223558.0
5 212000 228800.0
6 221000 257087.0
7 102000 108516.0
8 290000 328641.0
9 140000 156622.0

早速4つ目のデータは大外ししていますね。ほぼ半分。
これではスコアもあまり期待できなさそうですが、まずは確認したいと思います。
切り分けておいたテストデータで提出用のSalePrice表を作ります。

y_pred_final = rf_gs.predict(test)

次に提出用データのフォーマットを確認し、合わせます。

pd.read_csv('sample_submission.csv').head()
Id SalePrice
0 1461 169277.052498
1 1462 187758.393989
2 1463 183583.683570
3 1464 179317.477511
4 1465 150730.079977

予測値と、最初の方で分けておいたIdを合体させる。

submission = pd.concat((df_test_index, pd.DataFrame(y_pred_final)), axis=1)
submission.columns = ['Id', 'SalePrice']
#確認
submission.head()
Id SalePrice
0 1461 129225.399612
1 1462 159871.544036
2 1463 173464.563263
3 1464 183416.153801
4 1465 208079.458478

それでは、csvに吐き出してKaggleにアップロードします。
本来もっともワクワクする瞬間ですが、今回ベンチマーク編ではあまり期待できないのがわかっている為流れ作業的に...

submission.to_csv("submission.csv",sep=',',index=False)

では、一体スコアは何点だったでしょうか?
f:id:prml2synth:20180601225945p:plain


結果は2690位で、2018/6/1時点で参加者は5458チームのため、ちょうど真ん中あたりですね。

ベンチマーク編では予想通りスコアは残念な結果となってしまいましたが、次の発展編でデータ探査、精度評価及びチューニングを追求し、
どこまで順位を伸ばせるのか試してみたいと思います。

発展編1(上位20%)書きました。
www.mizugarasu.jp
発展編2(上位3%)やってみました。
www.mizugarasu.jp