Kaggleコンペ「HousePrices」~(1)まずは予測値を出すまでの一連の流れを実装してみる~
※統計・数学についてまだ知識が及んでいない部分があるため、一部間違った解釈をしている可能性も御座います。ご指摘・質問など御座いましたらコメントに書いて頂けると有り難いです。
今回取り組むテーマ
House Prices: Advanced Regression Techniques | Kaggle
駐車場の広さとか築年数などの説明変数を基に家がいくらで売れるか予測するという、スタンダードなタスクです。このタスクは親切にもtutorialが用意されており(Comprehensive data exploration with Python | Kaggle)、こちらでデータの概要や必要な前処理を把握することが出来ます。ざっと見た感じ、このタスクの難しそうなところは
- 変数選択(説明変数が78個あり、どれを使うか・どのように組み合わせて特徴量を作るかが肝になりそう)
- 手法選択(予測モデルを比較したり・モデルアンサンブルしたりしてモデルを決めていくことになりそう)
- 欠損値処理(訓練データ・提出データ共に欠損値が存在し、それをどのように扱うかも肝になりそう)
- データ数が少ない事に対する対策(訓練データが1460しかない)
が挙げられるのですが、まずは予測結果を出すまでの流れを作るのが大事だと思うので、細かいことはスルーして一連の流れを作っていこうと思います。
実装
データの読み込み・前処理
まずはモジュールとデータの読み込みを行っていきます。
#モジュールの読み込み import math import pandas as pd import matplotlib.pyplot as plt import seaborn as sns import numpy as np from scipy.stats import norm from sklearn.preprocessing import StandardScaler from scipy import stats import warnings from sklearn import tree# 決定木 from sklearn.tree import DecisionTreeRegressor from sklearn.model_selection import train_test_split #テストデータ分割 from sklearn.tree import export_graphviz #決定木の結果描画用 warnings.filterwarnings('ignore') %matplotlib inline #データの読み込み test = pd.read_csv("test.csv") train = pd.read_csv("train.csv") test_bp = test train_bp = train
次に前処理ですが、tutorialを参考に以下の事を実行しました。
- 変数選択
- 欠損値処理
- 異常値削除
- データが正規性を持つように対数変換
変数選択
今回は「目的変数と相関係数が高い9個に絞ったあと、説明変数同士で相関係数が0.8以上の組み合わせがあったら片方を削除する(多重共線性を回避するため)」という非常に荒い方法を取りました。この方法もゆくゆくは深掘っていきます。
#データの前処理 #変数選択(一旦相関係数が高い9個に絞る) k = 10 #number of variables for heatmap cols = corrmat.nlargest(k, 'SalePrice')['SalePrice'].index cm = np.corrcoef(train[cols].values.T) sns.set(font_scale=1) hm = sns.heatmap(cm, cbar=True, annot=True, square=True, fmt='.2f', annot_kws={'size': 10}, yticklabels=cols.values, xticklabels=cols.values) plt.show() train = train[cols] ##多重共線性を起こす可能性のある変数を削除 train = train.drop(["TotRmsAbvGrd","GarageArea","1stFlrSF"],axis = 1)
各変数同士の相関係数一覧です。TotRmsAbvGrdとGrLivArea,GarageAreaとGarageCars,1stFlrSFとTotalBsmtSFの相関係数が高いので、片方の変数を落として計6個の説明変数で予測を行っていきます。
欠損値処理
欠損値がある場合はリストワイズ削除するなり補完するなりする必要がありますが、今回は欠損値はなかったので何もしません。
#欠損値の確認(今回は欠損値無し) total = train.isnull().sum().sort_values(ascending=False) percent = (train.isnull().sum()/train.isnull().count()).sort_values(ascending=False) missing_data = pd.concat([total, percent], axis=1, keys=['Total', 'Percent']) missing_data.head(20)
異常値チェック
pairplotで各変数同士の散布図を確認し、挙動がおかしい2レコードを削除しました。
データが正規性を持つように対数変換
この部分は何をしているかというと、「正規分布になっていない変数に対して対数変換を行うことで正規分布に近づけている」という感じです。一般的に歪度が正である分布(右裾が長い)に対しては対数変換をすることで正規分布に近づく事が知られているため、対数変換を行っています。この処理を行う理由は正直完璧に理解できていないのですが、データに正規性を持たせることで「t検定など一般的な統計検定を行うことが出来る」「分散不均一性などの問題を回避できる(=時系列を考える時に関わってくる?)」というメリットがあるようです。
#正規性を持たせる前の分布 sns.distplot(train['SalePrice'], fit=norm);
このように目的変数(SalePrice)は右裾が長い分布になっているので、対数変換を行って正規性を持たせます。
#正規分布に近づけるために対数変換 train['SalePrice'] = np.log(train['SalePrice']) sns.distplot(train['SalePrice'], fit=norm);[f:id:avbio:20180506173538p:plain] #fig = plt.figure() res_saleprice = stats.probplot(train['SalePrice'], plot=plt)
他の説明変数に関しても変換が必要なものに対して行います。(TotalBsmtSFは0のレコードが存在し、0を対数変換することは出来ないので0の場合だけ避けるという少してくい感じの処理をしています。)
train['GrLivArea'] = np.log(train['GrLivArea']) sns.distplot(train['GrLivArea'], fit=norm); #fig = plt.figure() res_grLivarea = stats.probplot(train['GrLivArea'], plot=plt) train['HasBsmt'] = pd.Series(len(train['TotalBsmtSF']), index=train.index) train['HasBsmt'] = 0 train.loc[train['TotalBsmtSF']>0,'HasBsmt'] = 1 train.loc[train['HasBsmt']==1,'TotalBsmtSF'] = np.log(train['TotalBsmtSF']) sns.distplot(train[train['TotalBsmtSF']>0]['TotalBsmtSF'], fit=norm); #fig = plt.figure() res_hasbsmt = stats.probplot(train[train['TotalBsmtSF']>0]['TotalBsmtSF'], plot=plt) train = train.drop(['HasBsmt'],axis = 1)
モデル構築・検証・テストデータに対する予測結果算出
さて、前処理が完了したので予測するタスクに入っていきます。回帰の手法として今回は決定木を用いました。この段階では提出データ(train)には手を付けずに訓練データを分割して訓練データとテストデータに分けて汎化性能を確認するという方法で進めます。
#訓練データを用意 X = train.drop("SalePrice", axis=1) Y = train['SalePrice'] # 学習データとテストデータに分ける X_train, X_test, y_train, y_test = train_test_split(X,Y, random_state=50) # 決定木インスタンス（深さ5） tree_model = DecisionTreeRegressor(max_depth=5, random_state=50) tree_model.fit(X_train,y_train) print("train:",tree_model.__class__.__name__ ,tree_model.score(X_train,y_train)) print("test:",tree_model.__class__.__name__ , tree_model.score(X_test,y_test)) X_train.head()
train: DecisionTreeRegressor 0.833640464451 test: DecisionTreeRegressor 0.769666098077
こちらが予測結果です。訓練データ…83.4%,テストデータ…77.0%なので、過学習起こしているみたいですね…。原因としては、訓練データが少なすぎる事・枝分かれの判断基準や層の深さについて最適化されていない・特徴量に問題があるなどがありそうです。こちらも追って対策していきますが一旦突き進みます。
木の全体像も可視化しておきます(可視化するには別途Graphvizのインストールが必要になります。Graphviz をインストールするが参考になりました、ありがとうございます。
)
#決定木の可視化 temp_cols = cols[1:9] export_graphviz(tree_model, out_file='tree.dot', feature_names=X.columns)
上記プログラム実行後、「$ dot -Tpng tree.dot -o tree.png」を叩いて画像を出力することが出来ます。全体像だと細かくなってしまって分かりづらいので、全体像と一部抜粋した画像を載せます。
最後にこのモデルを提出データに適用させ、予測結果を算出します。
#提出データを使用している説明変数のみに絞ったあと、当てはめてみる temp_cols = train.columns temp_cols = temp_cols[1:7] test = test[temp_cols] #欠損値を平均値で置き換えた後、予測価格を出力する test = test.fillna(test.mean()) #正規性を持たせるために対数変換しているので、最後にネイピア数を掛けて価格に戻す predict = math.e ** tree_model.predict(test) f = open("output.csv","w") writer = csv.writer(f, lineterminator='\n') writer.writerow(predict) f.close()
あとは結果をkaggleのHPのSubmit Predictionsから提出すれば完了です!!
提出結果
Score = 0.370,1285位/15363位という結果でした。思ったより上位にいったみたいでびっくり。しかしながら最高スコアが289.94点なので、ここから上がっていくのが至難の技なのでしょう。こうやって予測結果の順位がすぐ見れるのは、モチベーションが上がりますね。
まとめ
今回は予測結果を出すための一連の流れを実装してみました。そこそこの順位には食い込めたのですが改良の余地は多いです。今後は下記の点を修正していく事になるかなと思います。
